/*C * Original project: Lars Arge, Jeff Chase, Pat Halpin, Laura Toma, Dean * Urban, Jeff Vitter, Rajiv Wickremesinghe 1999 * * GRASS Implementation: Lars Arge, Helena Mitasova, Laura Toma 2002 * * Copyright (c) 2002 Duke University -- Laura Toma * * Copyright (c) 1999-2001 Duke University -- * Laura Toma and Rajiv Wickremesinghe * * Redistribution and use in source and binary forms, with or without * modification, are permitted provided that the following conditions * are met: * 1. Redistributions of source code must retain the above copyright * notice, this list of conditions and the following disclaimer. * 2. Redistributions in binary form must reproduce the above copyright * notice, this list of conditions and the following disclaimer in the * documentation and/or other materials provided with the distribution. * 3. All advertising materials mentioning features or use of this software * must display the following acknowledgement: * This product includes software developed by Duke University * 4. Neither the name of the University nor the names of its contributors * may be used to endorse or promote products derived from this software * without specific prior written permission. * * THIS SOFTWARE IS PROVIDED BY THE TRUSTEES AND CONTRIBUTORS ``AS IS'' AND * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE * ARE DISCLAIMED. IN NO EVENT SHALL THE TRUSTEES OR CONTRIBUTORS BE LIABLE * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF * SUCH DAMAGE. *C*/ #include #include #include #include "fill.h" #include "option.h" #include "common.h" #include "water.h" #include "sortutils.h" #include "streamutils.h" #include "filldepr.h" #include "grid.h" /* globals in common.H extern statsRecorder *stats; stats file extern userOptions *opt; command-line options extern struct Cell_head *region; header of the region extern dimension_type nrows, ncols; */ #define FILL_DEBUG if(0) #define FILL_SAVEALL if(0) /* defined in this module */ void recordWatersheds(AMI_STREAM *labeledWater); long assignDirections(AMI_STREAM *statstr, AMI_STREAM *platstr, AMI_STREAM *waterstr); void assignFinalDirections(AMI_STREAM *statstr, AMI_STREAM *platstr, AMI_STREAM *waterstr); template void mergeStreamGridGrid(AMI_STREAM *grid1, AMI_STREAM *grid2, dimension_type rows, dimension_type cols, AMI_STREAM *str, FUN fo, AMI_STREAM *outStream); void merge2waterBase(AMI_STREAM *unsortedWaterStr, AMI_STREAM *dirStr, AMI_STREAM *elStr, AMI_STREAM *merge); AMI_STREAM* merge2waterGrid(AMI_STREAM *unsortedWaterStr, AMI_STREAM *dirStr, AMI_STREAM *elStr); AMI_STREAM * findBoundariesMain(AMI_STREAM *labeledWater); /* ********************************************************************** */ /* some helper classes */ /* ********************************************************************** */ class printElevation { public: char *operator()(const elevation_type &p) { static char buf[20]; sprintf(buf, "%.1f", (float)p); return buf; } }; class printDirection { public: char *operator()(const direction_type &p) { static char buf[20]; sprintf(buf, "%3d", p); return buf; } char *operator()(const waterWindowBaseType &p) { static char buf[20]; sprintf(buf, "%3d", p.dir); return buf; } #if(0) char *operator()(const waterWindowBaseType &p) { static char buf[3]; buf[0] = directionSymbol(p.dir); buf[1] = '\0'; return buf; } #endif }; class printLabel { public: char *operator()(const labelElevType&p) { static char buf[8]; sprintf(buf, CCLABEL_FMT, p.getLabel()); return buf; } char *operator()(const waterGridType &p) { static char buf[8]; sprintf(buf, CCLABEL_FMT, p.getLabel()); return buf; } char *operator()(const waterType &p) { static char buf[8]; sprintf(buf, CCLABEL_FMT, p.getLabel()); return buf; } }; class printDepth { public: char *operator()(const waterGridType &p) { static char buf[3]; sprintf(buf, "%1u", p.depth); return buf; } }; char * verbosedir(char *s) { static char buf[BUFSIZ]; sprintf(buf, "dump/%s", s); return buf; } /* ---------------------------------------------------------------------- */ /* fill the terrain (if necessary) and compute flow direction stream; elstr must exist and contain grid data before call, filledstr and dirstr are created; elstr is deleted and replaced with the classified elstr, which has boundary nodata distinguished from inner nodata */ AMI_STREAM* computeFlowDirections(AMI_STREAM*& elstr, AMI_STREAM*& filledstr, AMI_STREAM*& dirstr, AMI_STREAM *& labeledWater) { Rtimer rt, rtTotal; AMI_STREAM *elstr_reclass=NULL; AMI_STREAM *winstr=NULL; AMI_STREAM *statstr=NULL; AMI_STREAM *platstr=NULL; AMI_STREAM *waterstr=NULL; AMI_STREAM *mergedWaterStr=NULL; AMI_STREAM *boundaryStr=NULL; AMI_STREAM *waterWindows=NULL; rt_start(rtTotal); assert(elstr && filledstr == NULL && dirstr == NULL && labeledWater == NULL); stats->comment("------------------------------"); stats->comment("COMPUTING FLOW DIRECTIONS"); /* adjust nodata -- boundary nodata distinguished from inner nodata */ stats->comment("classifying nodata (inner & boundary)"); elstr_reclass = classifyNodata(elstr); delete elstr; elstr = elstr_reclass; /* ---------------------------------------------------------------------- */ /* find the plateaus. */ /* ---------------------------------------------------------------------- */ stats->comment("----------", opt->verbose); stats->comment("assigning preliminary directions"); rt_start(rt); dirstr = new AMI_STREAM; winstr = new AMI_STREAM(); statstr = new AMI_STREAM; platstr = findPlateaus(elstr, nrows, ncols, nodataType::ELEVATION_NODATA, winstr, dirstr, statstr); delete winstr; /* not used; not made */ rt_stop(rt); stats->recordTime("findingPlateaus", rt); stats->recordLength("plateaus", platstr); stats->recordLength("plateau stats", statstr); FILL_SAVEALL { /* printStream(*stats, statstr); */ FILL_DEBUG cout << "sort plateauStr (by ij): "; AMI_STREAM *tmp = sort(platstr, ijCmpPlateauType()); printStream2Grid(tmp, nrows, ncols, verbosedir("label1.asc"), plateauType::printLabel); delete tmp; } /* ---------------------------------------------------------------------- */ /* assign labels and directions & BFS ordering. depressions have labels, but no direction information. */ /* ---------------------------------------------------------------------- */ rt_start(rt); waterstr = new AMI_STREAM(); assignDirections(statstr, platstr, waterstr); delete platstr; delete statstr; rt_stop(rt); stats->recordTime("assigning directions", rt); *stats << "maxWatershedCount=" << labelFactory::getLabelCount() << endl; rt_start(rt); mergedWaterStr = merge2waterGrid(waterstr, dirstr, elstr); delete dirstr; delete waterstr; rt_stop(rt); stats->recordTime("merging", rt); stats->recordLength("mergedWaterStr", mergedWaterStr); /* ---------------------------------------------------------------------- */ /* watershed analysis */ /* IN: mergedWaterStr, labelFactory::... */ /* ---------------------------------------------------------------------- */ stats->comment("--------------", opt->verbose); stats->comment("generating watersheds and watershed graph"); rt_start(rt); waterWindows = new AMI_STREAM(); createWaterWindows(mergedWaterStr, nrows, ncols, waterWindows); delete mergedWaterStr; rt_stop(rt); stats->recordTime("creating windows", rt); stats->recordLength("waterWindows", waterWindows); /* ---------------------------------------------------------------------- */ rt_start(rt); labeledWater = new AMI_STREAM(); boundaryStr = new AMI_STREAM(); generateWatersheds(&waterWindows, nrows, ncols, labeledWater, boundaryStr); /* do we need to make boundaries here?? */ delete waterWindows; /* cout << "bogus boundary length = " << boundaryStr->stream_len() << endl;*/ assert(boundaryStr->stream_len() == 0); delete boundaryStr; assert(labeledWater->stream_len() == nrows * ncols); rt_stop(rt); stats->recordTime("generating watersheds", rt); /* ---------------------------------------------------------------------- */ /* find boundaries */ /* ---------------------------------------------------------------------- */ FILL_DEBUG cerr << "sort labeledWater (by ij):"; sort(&labeledWater, ijCmpLabelElevType()); #ifdef SAVE_ASCII cerr << "saving WATERSHED grid as watershed_grid\n"; printStream2Grid(labeledWater, nrows, ncols, "watershed.asc", labelElevType::printLabel); #endif boundaryStr = findBoundariesMain(labeledWater); /* ---------------------------------------------------------------------- */ /* filling */ /* valid streams are: boundaryStr, labeledWater */ /* ---------------------------------------------------------------------- */ rt_start(rt); elevation_type *raise; /*find the raise elevations */ FILL_DEBUG cerr << "sort boundaryStr (by elev): "; sort(&boundaryStr, elevCmpBoundaryType()); raise = fill_depression(boundaryStr, labelFactory::getLabelCount()); delete boundaryStr; rt_stop(rt); stats->recordTime("filling depressions", rt); /*fill the terrain*/ rt_start(rt); filledstr = new AMI_STREAM(); commit_fill(labeledWater, raise, labelFactory::getLabelCount(), filledstr); assert(filledstr->stream_len() == nrows * ncols); delete [] raise; rt_stop(rt); stats->recordTime("updating filled grid", rt); stats->recordLength("filledstr", filledstr); /* ---------------------------------------------------------------------- */ /* find plateaus again and reassign directions */ /* ---------------------------------------------------------------------- */ stats->comment("------------------------------"); stats->comment("REASSIGNING DIRECTIONS"); rt_start(rt); winstr = NULL; dirstr = new AMI_STREAM(); statstr = new AMI_STREAM(); platstr = findPlateaus(filledstr, nrows, ncols, nodataType::ELEVATION_NODATA, winstr, dirstr, statstr); rt_stop(rt); stats->recordTime("findingPlateaus2", rt); stats->recordLength("final plateaus", platstr); stats->recordLength("final plateau stats", statstr); FILL_SAVEALL { FILL_DEBUG cout << "sort plateauStr (by ij): "; AMI_STREAM *tmp = sort(platstr, ijCmpPlateauType()); printStream2Grid(tmp, nrows, ncols, verbosedir("plateaus.asc"), plateauType::printLabel); delete tmp; } /* assign final directions */ rt_start(rt); waterstr = new AMI_STREAM(); long dc = assignDirections(statstr, platstr, waterstr); if(dc) { *stats << "WARNING: " << dc << " depressions (islands) detected\n"; } delete platstr; delete statstr; rt_stop(rt); stats->recordTime("final directions", rt); /* merge */ rt_start(rt); AMI_STREAM *flowStream; /*STREAM_TO_OPTION(flowStream, "flowStream");*/ char path[BUFSIZ]; char* base_dir = getenv(STREAM_TMPDIR); assert(base_dir); sprintf(path, "%s/flowStream", base_dir); flowStream = new AMI_STREAM(path); /*flowStream->persist(PERSIST_PERSISTENT); */ stats->comment("creating flowStream: ", flowStream->sprint()); merge2waterBase(waterstr, dirstr, filledstr, flowStream); delete waterstr; rt_stop(rt); stats->recordTime("merge water,dir,elev to flow", rt); rt_stop(rtTotal); #ifdef SAVE_ASCII /*write grids as ascii file */ printGridStream(filledstr, nrows, ncols, "filled_elev.asc", printElevation()); printGridStream(flowStream, nrows, ncols, "direction.asc", printDirection()); #endif stats->recordTime("Total compute flow direction running time", rtTotal); stats->comment("compute flow directions done."); return flowStream; } /* ---------------------------------------------------------------------- */ void recordWatersheds(AMI_STREAM *labeledWater) { AMI_err ae; long labelCount = 0; AMI_STREAM *tmp; *stats << "--- watershed stats" << endl; FILL_DEBUG cout << "sort labeledWater (by wat label): "; tmp = sort(labeledWater, labelCmpLabelElevType()); labelElevType *p; cclabel_type prev(LABEL_UNDEF); tmp->seek(0); while((ae = tmp->read_item(&p)) == AMI_ERROR_NO_ERROR) { if(p->getLabel() != prev) { labelCount++; prev = p->getLabel(); } } assert(ae == AMI_ERROR_END_OF_STREAM); *stats << "watershed count = " << labelCount << endl; *stats << "---" << endl; stats->flush(); delete tmp; } /* ********************************************************************** */ /* assign directions to plateaus that have sinks; * reassign labels to depressions (don't drain out). * all plateaus are written out to the waterstr. */ long assignDirections(AMI_STREAM *statstr, AMI_STREAM *platstr, AMI_STREAM *waterstr) { size_t fmem; AMI_err ae; plateauStats *ps; stats->comment("----------", opt->verbose); stats->comment("assigning directions on plateaus"); labelFactory::reset(); /* we are relabeling now */ statstr->seek(0); platstr->seek(0); fmem = getAvailableMemory(); long depressionCount=0; long spillCount=0; while((ae = statstr->read_item(&ps)) == AMI_ERROR_NO_ERROR) { if(ps->size*sizeof(gridElement) > fmem) { cerr << "WARNING: grid larger than memory (ignored)" << endl; } assert(ps->label != LABEL_NODATA); if(ps->hasSpill) { spillCount++; grid *platGrid = new grid(ps->iMin, ps->jMin, ps->iMax, ps->jMax, ps->size, ps->label); platGrid->load(*platstr); platGrid->assignDirections(opt->d8 ? 1 : 0); platGrid->save(*waterstr); /* this doesn't save labels */ delete platGrid; } else { /* depression - just give contiguous labels only */ depressionCount++; cclabel_type label = labelFactory::getNewLabel(); for(int i=0; isize; i++) { plateauType *pt; platstr->read_item(&pt); pt->cclabel = label; /* assign new label */ waterType wt(*pt); /* write it out */ ae = waterstr->write_item(wt); assert(ae == AMI_ERROR_NO_ERROR); } } } *stats << "depression count = " << depressionCount << endl; *stats << "spill count = " << spillCount << endl; return depressionCount; } /* ********************************************************************** */ /* assign directions to plateaus that have sinks; * check that there are no depressions. */ void assignFinalDirections(AMI_STREAM *statstr, AMI_STREAM *platstr, AMI_STREAM *waterstr) { AMI_err ae; plateauStats *ps; stats->comment("assigning final directions"); statstr->seek(0); platstr->seek(0); while((ae = statstr->read_item(&ps)) == AMI_ERROR_NO_ERROR) { if(ps->hasSpill) { grid *platGrid = new grid(ps->iMin, ps->jMin, ps->iMax, ps->jMax, ps->size, ps->label); platGrid->load(*platstr); platGrid->assignDirections(opt->d8 ? 1 : 0); platGrid->save(*waterstr); /* this doesn't save labels */ delete platGrid; } else { /* could be legitimate */ cerr << "WARNING: depression detected: " << *ps << endl; continue; } } }; /* ********************************************************************** */ class directionElevationMerger { public: waterGridType operator()(elevation_type el, direction_type dir, waterType p) { /* check that no (boundary) nodata values got in here */ assert(el != nodataType::ELEVATION_BOUNDARY); assert(!is_nodata(el)); /* p should be a valid grid cell */ return waterGridType(el, p.getDirection(), p.getLabel(), p.getDepth()); } waterGridType operator()(elevation_type el, direction_type dir) { waterGridType wg(el, dir); if(el == nodataType::ELEVATION_BOUNDARY) { /* hack XXX (approved) */ wg.setLabel(LABEL_BOUNDARY); } /* nodata => boundary or undef */ assert(!is_nodata(el) || (wg.getLabel()==LABEL_BOUNDARY || wg.getLabel()==LABEL_UNDEF)); return wg; } }; /* ---------------------------------------------------------------------- */ AMI_STREAM * merge2waterGrid(AMI_STREAM *unsortedWaterStr, AMI_STREAM *dirStr, AMI_STREAM *elStr) { AMI_STREAM *waterStr; FILL_DEBUG cout << "sort waterStr (by ij): "; waterStr = sort(unsortedWaterStr, ijCmpWaterType()); FILL_SAVEALL printStream2Grid(waterStr, nrows, ncols, verbosedir("platlabels.asc"), printLabel()); AMI_STREAM *mergedWaterStr = new AMI_STREAM(); mergeStreamGridGrid(elStr, dirStr, nrows, ncols, waterStr, directionElevationMerger(), mergedWaterStr); delete waterStr; FILL_SAVEALL printGridStream(mergedWaterStr, nrows, ncols, verbosedir("mergedlabels.asc"), printLabel()); assert(mergedWaterStr->stream_len()); return mergedWaterStr; } /* ---------------------------------------------------------------------- */ void merge2waterBase(AMI_STREAM *unsortedWaterStr, AMI_STREAM *dirStr, AMI_STREAM *elStr, AMI_STREAM *merge) { AMI_STREAM *waterStr; FILL_DEBUG cout << "sort waterStr (by ij): "; waterStr = sort(unsortedWaterStr, ijCmpWaterType()); mergeStreamGridGrid(elStr, dirStr, nrows, ncols, waterStr, directionElevationMerger(), merge); delete waterStr; } /* ---------------------------------------------------------------------- */ /* * merge 2 grids and a stream together to form a new grid. * str should be sorted in ij order */ template void mergeStreamGridGrid(AMI_STREAM *grid1, AMI_STREAM *grid2, dimension_type rows, dimension_type cols, AMI_STREAM *str, FUN fo, AMI_STREAM *outStream) { T1 *t1p; T2 *t2p; T3 *t3p; AMI_err aeUpd, ae; grid1->seek(0); grid2->seek(0); str->seek(0); aeUpd = str->read_item(&t3p); assert(aeUpd == AMI_ERROR_NO_ERROR || aeUpd == AMI_ERROR_END_OF_STREAM); for(dimension_type row = 0; row < rows; row++) { for(dimension_type col = 0; col < cols; col++) { ae = grid1->read_item(&t1p); assert(ae == AMI_ERROR_NO_ERROR); ae = grid2->read_item(&t2p); assert(ae == AMI_ERROR_NO_ERROR); T4 t4; if(aeUpd == AMI_ERROR_NO_ERROR && t3p->i == row && t3p->j == col) { /* cell present in stream */ t4 = fo(*t1p, *t2p, *t3p); aeUpd = str->read_item(&t3p); assert(aeUpd == AMI_ERROR_NO_ERROR || aeUpd == AMI_ERROR_END_OF_STREAM); } else { /* not in stream */ t4 = fo(*t1p, *t2p); } ae = outStream->write_item(t4); assert(ae == AMI_ERROR_NO_ERROR); } /*assert(outStream->stream_len() == (row+1) * cols); */ } assert(outStream->stream_len() == rows * cols); return; } /* ---------------------------------------------------------------------- */ /* make boundaryStr from labeledWater */ AMI_STREAM * findBoundariesMain(AMI_STREAM *labeledWater) { AMI_STREAM *boundaryStr; Rtimer rt; rt_start(rt); boundaryStr = new AMI_STREAM(); FILL_SAVEALL printGridStream(labeledWater, nrows, ncols, verbosedir("labels.asc"), printLabel()); findBoundaries(labeledWater, nrows, ncols, boundaryStr); stats->recordLength("all boundaries", boundaryStr); FILL_SAVEALL { FILL_DEBUG cout << "sort boundaryStr (by ij): "; sort(&boundaryStr, ijCmpBoundaryType()); removeDuplicatesEx(&boundaryStr, ijCmpBoundaryType()); printStream2Grid(boundaryStr, nrows, ncols, verbosedir("boundary.asc"), boundaryType::print); } FILL_DEBUG cout << "sort boundaryStr (by wat label): "; sort(&boundaryStr, waterCmpBoundaryType()); removeDuplicatesEx(&boundaryStr, boundaryCmpBoundaryType()); rt_stop(rt); stats->recordTime("generating boundaries", rt); stats->recordLength("boundary stream", boundaryStr); /*if(GETOPT("veryfillverbose")) printStream(cout, boundaryStr); SAVE_ON_OPTION(boundaryStr, "saveBoundaryStream"); */ return boundaryStr; }