Skip to content

Commit 0a959ca

Browse files
dtarbdtarb
committed
Add minimum area capability to catchment outlets
1 parent 87e37d2 commit 0a959ca

File tree

3 files changed

+44
-26
lines changed

3 files changed

+44
-26
lines changed

‎src/CatchOutlets.cpp

+31-24
Original file line numberDiff line numberDiff line change
@@ -75,6 +75,7 @@ struct netlink {
7575
int32_t uslinkno1;
7676
int32_t uslinkno2;
7777
double doutend;
78+
double uscontarea;
7879
double x, y;
7980
};
8081
OGRLayerH hLayer, hLayerpt;
@@ -83,7 +84,7 @@ OGRGeometryH hGeometry, hGeometrypt;
8384

8485
netlink* allinks;
8586

86-
void CheckPoint(int thelink, double downout, double mindist)
87+
void CheckPoint(int thelink, double downout, double mindist, double minarea)
8788
{
8889
int istream=0;
8990
while (linknos[istream] != thelink && istream < nstreams) {
@@ -92,33 +93,37 @@ void CheckPoint(int thelink, double downout, double mindist)
9293
if (istream < nstreams) // Here a stream was found. This traps the passing over streams with none upstream
9394
{
9495
if (allinks[istream].doutend - downout > mindist) {
95-
//printf("%d, %g, %g\n", linknos[istream], allinks[istream].x, allinks[istream].y);
96-
hGeometrypt = OGR_G_CreateGeometry(wkbPoint);// create geometry
97-
OGRFeatureH hFeaturept = OGR_F_Create(OGR_L_GetLayerDefn(hLayerpt)); // create new feature with null fields and no geometry
98-
OGR_F_SetFieldInteger(hFeaturept, 0, linknos[istream]);
99-
OGR_F_SetFieldInteger(hFeaturept, 1, allinks[istream].dslinkno);
100-
OGR_F_SetFieldInteger(hFeaturept, 2, allinks[istream].uslinkno1);
101-
OGR_F_SetFieldInteger(hFeaturept, 3, allinks[istream].uslinkno2);
102-
OGR_F_SetFieldDouble(hFeaturept, 4, allinks[istream].doutend);
103-
OGR_F_SetFieldInteger(hFeaturept, 5, gwcount);
104-
gwcount = gwcount + 1;
105-
OGR_G_SetPoint_2D(hGeometrypt, 0, allinks[istream].x, allinks[istream].y);
106-
OGR_F_SetGeometry(hFeaturept, hGeometrypt);
107-
OGR_G_DestroyGeometry(hGeometrypt);
108-
OGR_L_CreateFeature(hLayerpt, hFeaturept);
109-
OGR_F_Destroy(hFeaturept);
110-
111-
downout = allinks[istream].doutend;
96+
if (allinks[istream].uscontarea > minarea) { // Only write if min area is exceeded
97+
//printf("%d, %g, %g\n", linknos[istream], allinks[istream].x, allinks[istream].y);
98+
hGeometrypt = OGR_G_CreateGeometry(wkbPoint);// create geometry
99+
OGRFeatureH hFeaturept = OGR_F_Create(OGR_L_GetLayerDefn(hLayerpt)); // create new feature with null fields and no geometry
100+
OGR_F_SetFieldInteger(hFeaturept, 0, linknos[istream]);
101+
OGR_F_SetFieldInteger(hFeaturept, 1, allinks[istream].dslinkno);
102+
OGR_F_SetFieldInteger(hFeaturept, 2, allinks[istream].uslinkno1);
103+
OGR_F_SetFieldInteger(hFeaturept, 3, allinks[istream].uslinkno2);
104+
OGR_F_SetFieldDouble(hFeaturept, 4, allinks[istream].doutend);
105+
OGR_F_SetFieldInteger(hFeaturept, 5, gwcount);
106+
gwcount = gwcount + 1;
107+
OGR_G_SetPoint_2D(hGeometrypt, 0, allinks[istream].x, allinks[istream].y);
108+
OGR_F_SetGeometry(hFeaturept, hGeometrypt);
109+
OGR_G_DestroyGeometry(hGeometrypt);
110+
OGR_L_CreateFeature(hLayerpt, hFeaturept);
111+
OGR_F_Destroy(hFeaturept);
112+
113+
downout = allinks[istream].doutend;
114+
}
115+
}
116+
if (allinks[istream].uscontarea > minarea) { // Only recurse up if minarea is satisified
117+
if (allinks[istream].uslinkno1 >= 0) CheckPoint(allinks[istream].uslinkno1, downout, mindist, minarea);
118+
if (allinks[istream].uslinkno2 >= 0) CheckPoint(allinks[istream].uslinkno2, downout, mindist, minarea);
112119
}
113-
if (allinks[istream].uslinkno1 >= 0) CheckPoint(allinks[istream].uslinkno1, downout, mindist);
114-
if (allinks[istream].uslinkno2 >= 0) CheckPoint(allinks[istream].uslinkno2, downout, mindist);
115120
}
116121
}
117122

118123

119124
// int catchoutlets(char *pfile, char *streamnetsrc, char *streamnetlyr, char *outletsdatasrc, char *outletslayer, int lyrno, float maxdist)
120125

121-
int catchoutlets(char *pfile, char *streamnetsrc, char *outletsdatasrc, double mindist, int gwstartno)
126+
int catchoutlets(char *pfile, char *streamnetsrc, char *outletsdatasrc, double mindist, int gwstartno, double minarea)
122127
{
123128

124129
MPI_Init(NULL,NULL);{
@@ -257,8 +262,9 @@ int catchoutlets(char *pfile, char *streamnetsrc, char *outletsdatasrc, double m
257262
OGR_Fld_SetWidth(hFieldDefn, 10); // set field width
258263
OGR_L_CreateField(hLayerpt, hFieldDefn, 0);
259264

265+
// Fields we use but do not add to new feature class
260266
int32_t doutstartfield = OGR_FD_GetFieldIndex(hFDefn, "DOUTSTART");
261-
// Dont need to create new field for this
267+
int32_t uscontareafield = OGR_FD_GetFieldIndex(hFDefn, "USContArea");
262268

263269
int32_t istream = 0;
264270
gwcount = gwstartno; // Initialize count
@@ -293,6 +299,7 @@ int catchoutlets(char *pfile, char *streamnetsrc, char *outletsdatasrc, double m
293299
OGR_G_GetPoint(hGeometry, num_points-1, &X1, &Y1, &Z1);
294300
allinks[istream].doutend = OGR_F_GetFieldAsDouble(hFeature, doutstartfield); // Note here start field
295301
}
302+
allinks[istream].uscontarea = OGR_F_GetFieldAsDouble(hFeature, uscontareafield);
296303
allinks[istream].x = X1;
297304
allinks[istream].y = Y1;
298305
//printf("%d, %lf, %lf,%lf\n", istream, X1, Y1, Z1);
@@ -342,8 +349,8 @@ int catchoutlets(char *pfile, char *streamnetsrc, char *outletsdatasrc, double m
342349
OGR_F_Destroy(hFeaturept);
343350
double downout = 0.0;
344351
// Recursive calls to traverse up
345-
if (allinks[istream].uslinkno1 >= 0) CheckPoint(allinks[istream].uslinkno1, downout,mindist);
346-
if (allinks[istream].uslinkno2 >= 0) CheckPoint(allinks[istream].uslinkno2, downout,mindist);
352+
if (allinks[istream].uslinkno1 >= 0) CheckPoint(allinks[istream].uslinkno1, downout,mindist,minarea);
353+
if (allinks[istream].uslinkno2 >= 0) CheckPoint(allinks[istream].uslinkno2, downout,mindist,minarea);
347354
}
348355
}
349356
}

‎src/CatchOutlets.h

+1-1
Original file line numberDiff line numberDiff line change
@@ -2,5 +2,5 @@
22

33

44
//int catchoutlets(char *pfile, char *streamnetsrc, char *streamnetlyr, char *outletsdatasrc, char *outletslayer, int lyrno, float maxdist);
5-
int catchoutlets(char *pfile, char *streamnetsrc,char *outletsdatasrc, double mindist, int gwstartno);
5+
int catchoutlets(char *pfile, char *streamnetsrc,char *outletsdatasrc, double mindist, int gwstartno, double minarea);
66
// Incomplete implementation no capability for layers within data sources yet

‎src/CatchOutletsmn.cpp

+12-1
Original file line numberDiff line numberDiff line change
@@ -58,6 +58,7 @@ int main(int argc,char **argv)
5858
{
5959
char pfile[MAXLN],streamnetsrc[MAXLN], streamnetlyr[MAXLN]="",outletsdatasrc[MAXLN],outletslayer[MAXLN]="";
6060
double mindist = 0.0;
61+
double minarea = 0.0;
6162
int err,i,uselyrname=0,lyrno=0,gwstartno=1;
6263

6364
if(argc < 7)
@@ -151,12 +152,22 @@ int main(int argc,char **argv)
151152
}
152153
else goto errexit;
153154
}
155+
else if (strcmp(argv[i], "-minarea") == 0)
156+
{
157+
i++;
158+
if (argc > i)
159+
{
160+
sscanf(argv[i], "%lf", &minarea);
161+
i++;
162+
}
163+
else goto errexit;
164+
}
154165
else
155166
{
156167
goto errexit;
157168
}
158169
}
159-
if(err=catchoutlets(pfile,streamnetsrc,outletsdatasrc,mindist,gwstartno) != 0)
170+
if(err=catchoutlets(pfile,streamnetsrc,outletsdatasrc,mindist,gwstartno,minarea) != 0)
160171
printf("Catchment Outlets error %d\n",err);
161172

162173
return 0;

0 commit comments

Comments
 (0)