modify gagewatershed function for RWD from NHDPlus data

Author: nazmussazib

src/gagewatershed.cpp

@@ -53,7 +53,7 @@ using namespace std;
 
 
 
-int gagewatershed( char *pfile, char *swfile,char *wfile, char* datasrc,char* lyrname,int uselyrname,int lyrno, char *idfile,int writeid,int useswg,int writeupid, char *upidfile) 
+int gagewatershed( char *pfile,char *wfile, char* datasrc,char* lyrname,int uselyrname,int lyrno, char *idfile,int writeid, int writeupid,char *upidfile) 
 {//1
 
 	MPI_Init(NULL,NULL);{
@@ -68,7 +68,7 @@ int gagewatershed( char *pfile, char *swfile,char *wfile, char* datasrc,char* ly
 	int *ids=NULL;
 	int *dsids=NULL;
 	int idnodata;
-	
+
 	double begint = MPI_Wtime();
 	tiffIO p(pfile,SHORT_TYPE);
 	long totalX = p.getTotalX();
@@ -143,20 +143,6 @@ int gagewatershed( char *pfile, char *swfile,char *wfile, char* datasrc,char* ly
 	p.read(xstart, ystart, ny, nx, flowData->getGridPointer());
 	//printf("Pfile read");  fflush(stdout);
 
-	//if using subwatershed grid for rapid watershed delineation, get information from file
-	tdpartition *swData;
-	if( useswg == 1){
-		tiffIO sw(swfile,LONG_TYPE);
-		if(!p.compareTiff(sw)){
-			printf("File sizes do not match\n%s\n",swfile);
-			MPI_Abort(MCW,5);
-			return 1;  
-		} 
-		swData = CreateNewPartition(sw.getDatatype(), totalX, totalY, dxA, dyA, sw.getNodata()); 
-		swData->localToGlobal(0, 0, xstart, ystart);
-		sw.read(xstart, ystart, swData->getny(), swData->getnx(), swData->getGridPointer());
-	}
-
 	//Begin timer
 	double readt = MPI_Wtime();
 	//printf("Read time %lf\n",readt);
@@ -174,26 +160,22 @@ int gagewatershed( char *pfile, char *swfile,char *wfile, char* datasrc,char* ly
 	node temp;
 	queue<node> que;
 
-	
-
 	for( int i=0; i<numOutlets; i++)
 	{
 		p.geoToGlobalXY(x[i], y[i], outletsX[i], outletsY[i]);
 		int xlocal, ylocal;
 		wshed->globalToLocal(outletsX[i], outletsY[i],xlocal,ylocal);
-	
 		if(wshed->isInPartition(xlocal,ylocal)){   //xOutlets[i], yOutlets[i])){
 			wshed->setData(xlocal,ylocal,(int32_t)ids[i]);  //xOutlets[i], yOutlets[i], (long)ids[i]);
 							//Push outlet cells on to que
 						temp.x = xlocal;
 						temp.y = ylocal;
 						que.push(temp);	
-					
-
 		}
 		dsids[i]= idnodata;  
 	}
-	
+
+
 	// con is used to check for contamination at the edges
 	long i,j;
 	short k;
@@ -223,24 +205,20 @@ int gagewatershed( char *pfile, char *swfile,char *wfile, char* datasrc,char* ly
 	wshed->share();
 	// open upid file to write and update upstream watershed id file ( need to confirm)
 	FILE *fp1;
-
-	          fp1 = fopen(upidfile,"a");
-			  if (!rank) fprintf(fp1,"Upstream Contributing Watershed ID:\n"); 
+    fp1 = fopen(upidfile,"a");
+	if (!rank) fprintf(fp1,"Upstream Contributing coordinates:\n"); 
 
 	//flowData->getData(500,600,tempShort);
 	finished = false;
 	//Ring terminating while loop
 	while(!finished) {
 		while(!que.empty()){
 			//Takes next node with no contributing neighbors
-			int idex=0;
 			temp = que.front();
 			que.pop();
 			i = temp.x;
 			j = temp.y;	
-			
 			//  if not labeled, label with downstream result
-			
 			if(wshed->isNodata(i,j))
 			{
 				flowData->getData(i,j,k);
@@ -264,34 +242,16 @@ int gagewatershed( char *pfile, char *swfile,char *wfile, char* datasrc,char* ly
 						p.globalXYToGeo(gx,gy,wi,wj);
 						//printf("%f, %f", wi,wj);
 						fprintf (fp1,"%f, %f\n",wi,wj);
-						if(useswg==1);// {long sss=swData->getData(in,jn,tempLong);
-						//fprintf (fp1,"%ld\n",sss);}
+				
 				}
 				if(sdir > 0) 
 				{
-					int idex=0;
 					if((sdir-k==4 || sdir-k==-4))
 					{
 						//  add to Q
 						if(flowData->hasAccess(in,jn) && !flowData->isNodata(in,jn) && wshed->isNodata(in,jn)){
 							//Decrement the number of contributing neighbors in neighbor
 							neighbor->addToData(in,jn,(short)-1);
-							// for rapid watershed delineation
-							if(useswg==1){    
-							        // short swidup=swData->getData(in,jn,tempSubwgrid); // get upstream watershed ID
-							        // short swidme=swData->getData(i,j,tempSubwgrid); // get outlet point watershed ID
-									//  if upstream ID is not equal to outlet point watershed ID then write into upidfile 
-									//  swidup may be nodata value therefore also check is it greater than or equal to zero
-							       //  if((swidup!=swidme) & (swidup>=0)){
-								      //     write into upidfile
-										// printf("Upstream edge is reached\n");
-								       //  fprintf (fp1,"%d\n",swidup);
-									   //   } 
-									
-							}
-							
-
-							
 
 							//Check if neighbor needs to be added to que
 							if(flowData->isInPartition(in,jn) && neighbor->getData(in, jn, tempShort) == 0 ){
@@ -300,12 +260,10 @@ int gagewatershed( char *pfile, char *swfile,char *wfile, char* datasrc,char* ly
 								que.push(temp);
 							}
 						}
-	
 						if(!wshed->isNodata(in,jn))
 						{
 							int32_t idup=wshed->getData(in,jn,tempLong);
 							int32_t idme=wshed->getData(i,j,tempLong);
-							;
 							//  Find the array index for idup
 							int iidex;
 							for(iidex=0; iidex<numOutlets; iidex++)
@@ -321,7 +279,7 @@ int gagewatershed( char *pfile, char *swfile,char *wfile, char* datasrc,char* ly
 		//Pass information
 		neighbor->addBorders();
 		wshed->share();
-	  
+
 		//If this created a cell with no contributing neighbors, put it on the queue
 		for(i=0; i<nx; i++){
 			if(neighbor->getData(i, -1, tempShort)!=0 && neighbor->getData(i, 0, tempShort)==0){
@@ -342,18 +300,13 @@ int gagewatershed( char *pfile, char *swfile,char *wfile, char* datasrc,char* ly
 		finished = que.empty();
 		finished = wshed->ringTerm(finished);
 	}
-
 	//  Reduce all values to the 0 process
 	int *dsidsr;	
     dsidsr=new int[numOutlets];
 	MPI_Reduce(dsids, dsidsr,numOutlets,MPI_INT, MPI_MAX,0, MCW);
-//	MPI_Reduce(upids, dsidsr,numOutlets,MPI_INT, MPI_MAX,0, MCW);
 
 	//Stop timer
 	double computet = MPI_Wtime();
-	
-	
-	
 
 //  Write id.txt with from and to links
 	if(writeid == 1)
@@ -370,7 +323,6 @@ int gagewatershed( char *pfile, char *swfile,char *wfile, char* datasrc,char* ly
 		}
 	}
 
-
 	//Create and write TIFF file
 	long lNodata = MISSINGLONG;
 	tiffIO wshedIO(wfile, LONG_TYPE, &lNodata, p);
@@ -390,4 +342,3 @@ int gagewatershed( char *pfile, char *swfile,char *wfile, char* datasrc,char* ly
 
 
 
-

src/gagewatershedmn.cpp

@@ -141,25 +141,15 @@ int main(int argc,char **argv)
 			}
 			else goto errexit;
 		}
-		 else if(strcmp(argv[i],"-swg")==0)
-		{
-			i++;
-			if(argc > i)
-			{
-				strcpy(swfile,argv[i]);
-				useswg=1;
-				i++;
-			}
-			else goto errexit;
-	 }
+	
 
 		else 
 		{
 			goto errexit;
 		}
 	}
 
-    if( (err=gagewatershed(pfile,swfile,wfile,datasrc,lyrname,uselyrname,lyrno,idfile,writeid,useswg,writeupid,upidfile)) != 0)
+    if( (err=gagewatershed(pfile,wfile,datasrc,lyrname,uselyrname,lyrno,idfile,writeid,writeupid,upidfile)) != 0)
         printf("Gage watershed error %d\n",err);
 
 	return 0;