/**************************************************************************** * * MODULE: r.terraflow * * COPYRIGHT (C) 2007, 2010 Laura Toma * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by * the Free Software Foundation; either version 2 of the License, or * (at your option) any later version. * * This program is distributed in the hope that it will be useful, * but WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * GNU General Public License for more details. * *****************************************************************************/ #include #include #include #include #include #include #include #ifdef HAVE_STATVFS_H #include #endif extern "C" { #include #include #include } #include "option.h" #include "common.h" /* declares the globals */ #include "fill.h" #include "flow.h" #include "nodata.h" #include "grass2str.h" #include "water.h" #include "sortutils.h" /* globals: in common.H extern statsRecorder *stats; extern userOptions* opt; extern struct Cell_head region; */ /* #define JUMP2FLOW */ /* define it only if you want to skip the flow direction computation and jump directly to computing flow accumulation; the flowstream must exist in /STREAM_DIR/flowStream */ /* ---------------------------------------------------------------------- */ void parse_args(int argc, char *argv[]) { /* input elevation grid */ struct Option *input_elev; input_elev = G_define_standard_option(G_OPT_R_ELEV); /* output filled elevation grid */ struct Option *output_elev; output_elev = G_define_standard_option(G_OPT_R_OUTPUT); output_elev->key = "filled"; output_elev->description= _("Name for output filled (flooded) elevation raster map"); /* output direction grid */ struct Option *output_dir; output_dir = G_define_standard_option(G_OPT_R_OUTPUT); output_dir->key = "direction"; output_dir->description= _("Name for output flow direction raster map"); /* output sinkwatershed grid */ struct Option *output_watershed; output_watershed = G_define_standard_option(G_OPT_R_OUTPUT); output_watershed->key = "swatershed"; output_watershed->description= _("Name for output sink-watershed raster map"); /* output flow accumulation grid */ struct Option *output_accu; output_accu = G_define_standard_option(G_OPT_R_OUTPUT); output_accu->key = "accumulation"; output_accu->description= _("Name for output flow accumulation raster map"); #ifdef OUTPUT_TCI struct Option *output_tci; output_tci = G_define_standard_option(G_OPT_R_OUTPUT); output_tci->key = "tci"; output_tci->description= _("Name for output topographic convergence index (tci) raster map"); #endif /* MFD/SFD flag */ struct Flag *sfd_flag; sfd_flag = G_define_flag() ; sfd_flag->key = 's'; sfd_flag->description= _("SFD (D8) flow (default is MFD)"); /* D8CUT value*/ struct Option *d8cut; d8cut = G_define_option(); d8cut->key = "d8cut"; d8cut->type = TYPE_DOUBLE; d8cut->required = NO; d8cut->label = _("Routing using SFD (D8) direction"); d8cut->description = _("If flow accumulation is larger than this value it is routed using " "SFD (D8) direction (meaningfull only for MFD flow). " "If no answer is given it defaults to infinity."); /* main memory */ struct Option *mem; mem = G_define_option() ; mem->key = "memory"; mem->type = TYPE_INTEGER; mem->required = NO; mem->answer = "300"; /* 300MB default value */ mem->description = _("Maximum runtime memory size (in MB)"); /* temporary STREAM path */ struct Option *streamdir; streamdir = G_define_option() ; streamdir->key = "stream_dir"; streamdir->type = TYPE_STRING; streamdir->required = NO; //streamdir->answer = ""; streamdir->description= _("Directory to hold temporary files (they can be large)"); /* stats file */ struct Option *stats_opt; stats_opt = G_define_option() ; stats_opt->key = "stats"; stats_opt->type = TYPE_STRING; stats_opt->required = NO; stats_opt->description= _("Name of file containing runtime statistics"); if (G_parser(argc, argv)) { exit (EXIT_FAILURE); } /* ************************* */ assert(opt); opt->elev_grid = input_elev->answer; opt->filled_grid = output_elev->answer; opt->dir_grid = output_dir->answer; opt->watershed_grid = output_watershed->answer; opt->flowaccu_grid = output_accu->answer; #ifdef OUTPUT_TCI opt->tci_grid = output_tci->answer; #endif opt->d8 = sfd_flag->answer; if (!d8cut->answer) { opt->d8cut = MAX_ACCU; } else { opt->d8cut = atof(d8cut->answer); } opt->mem = atoi(mem->answer); if (!streamdir->answer) { const char *tmpdir = G_tempfile(); if (G_mkdir(tmpdir) == -1) G_fatal_error(_("Unable to create temp dir")); opt->streamdir = G_store(tmpdir); } else opt->streamdir = streamdir->answer; opt->verbose = G_verbose() == G_verbose_max(); opt->stats = stats_opt->answer; /* somebody should delete the options */ } /* ---------------------------------------------------------------------- */ /* check compatibility of map header and region header */ void check_header(char* cellname) { const char *mapset; mapset = G_find_raster(cellname, ""); if (mapset == NULL) { G_fatal_error(_("Raster map <%s> not found"), cellname); } /* read cell header */ struct Cell_head cell_hd; Rast_get_cellhd (cellname, mapset, &cell_hd); /* check compatibility with module region */ if (!((region->ew_res == cell_hd.ew_res) && (region->ns_res == cell_hd.ns_res))) { G_fatal_error(_("cell file %s resolution differs from current region"), cellname); } else { if (opt->verbose) { G_message(_("cell %s header compatible with region header"), cellname); fflush(stderr); } } /* check type of input elevation raster and check if precision is lost */ RASTER_MAP_TYPE data_type; data_type = Rast_map_type(opt->elev_grid, mapset); #ifdef ELEV_SHORT G_verbose_message(_("Elevation stored as SHORT (%dB)"), sizeof(elevation_type)); if (data_type == FCELL_TYPE) { G_warning(_("raster %s is of type FCELL_TYPE " "--precision may be lost."), opt->elev_grid); } if (data_type == DCELL_TYPE) { G_warning(_("raster %s is of type DCELL_TYPE " "--precision may be lost."), opt->elev_grid); } #endif #ifdef ELEV_FLOAT G_verbose_message( _("Elevation stored as FLOAT (%dB)"), sizeof(elevation_type)); if (data_type == CELL_TYPE) { G_warning(_("raster %s is of type CELL_TYPE " "--you should use r.terraflow.short"), opt->elev_grid); } if (data_type == DCELL_TYPE) { G_warning(_("raster %s is of type DCELL_TYPE " "--precision may be lost."), opt->elev_grid); } #endif } /* ---------------------------------------------------------------------- */ void check_args() { /* check if filled elevation grid name is valid */ if (G_legal_filename (opt->filled_grid) < 0) { G_fatal_error(_("<%s> is an illegal file name"), opt->filled_grid); } /* check if output grid names are valid */ if (G_legal_filename (opt->dir_grid) < 0) { G_fatal_error(_("<%s> is an illegal file name"), opt->dir_grid); } if (G_legal_filename (opt->filled_grid) < 0) { G_fatal_error(_("<%s> is an illegal file name"), opt->filled_grid); } if (G_legal_filename (opt->flowaccu_grid) < 0) { G_fatal_error(_("<%s> is an illegal file name"), opt->flowaccu_grid); } if (G_legal_filename (opt->watershed_grid) < 0) { G_fatal_error(_("<%s> is an illegal file name"), opt->watershed_grid); } #ifdef OUTPU_TCI if (G_legal_filename (opt->tci_grid) < 0) { G_fatal_error(_("<%s> is an illegal file name"), opt->tci_grid); } #endif /* check compatibility with region */ check_header(opt->elev_grid); /* what else ? */ } /* ---------------------------------------------------------------------- */ void record_args(int argc, char **argv) { time_t t = time(NULL); char buf[BUFSIZ]; if(t == (time_t)-1) { perror("time"); exit(1); } #ifdef __MINGW32__ strcpy(buf, ctime(&t)); #else ctime_r(&t, buf); buf[24] = '\0'; #endif stats->timestamp(buf); *stats << "Command Line: " << endl; for(int i=0; ielev_grid << "\n"; *stats << "output (flooded) elevations grid: " << opt->filled_grid << "\n"; *stats << "output directions grid: " << opt->dir_grid << "\n"; *stats << "output sinkwatershed grid: " << opt->watershed_grid << "\n"; *stats << "output accumulation grid: " << opt->flowaccu_grid << "\n"; #ifdef OUTPUT_TCI *stats << "output tci grid: " << opt->tci_grid << "\n"; #endif if (opt->d8) { stats ->comment("SFD (D8) flow direction"); } else { stats->comment("MFD flow direction"); } sprintf(buf, "D8CUT=%f", opt->d8cut); stats->comment(buf); size_t mm_size = (size_t) opt->mem << 20; /* (in bytes) */ char tmp[100]; formatNumber(tmp, mm_size); sprintf(buf, "Memory size: %s bytes", tmp); stats->comment(buf); } /* ---------------------------------------------------------------------- */ void setFlowAccuColorTable(char* cellname) { struct Colors colors; const char *mapset; struct Range r; mapset = G_find_raster(cellname, ""); if (mapset == NULL) { G_fatal_error (_("Raster map <%s> not found"), cellname); } if (Rast_read_range(cellname, mapset, &r) == -1) { G_fatal_error(_("cannot read range")); } /*fprintf(stderr, "%s range is: min=%d, max=%d\n", cellname, r.min, r.max);*/ int v[6]; v[0] = r.min; v[1] = 5; v[2] = 30; v[3] = 100; v[4] = 1000; v[5] = r.max; Rast_init_colors(&colors); Rast_add_c_color_rule(&v[0], 255,255,255, &v[1], 255,255,0, &colors); Rast_add_c_color_rule(&v[1], 255,255,0, &v[2], 0,255,255, &colors); Rast_add_c_color_rule(&v[2], 0,255,255, &v[3], 0,127,255, &colors); Rast_add_c_color_rule(&v[3], 0,127,255, &v[4], 0,0,255, &colors); Rast_add_c_color_rule(&v[4], 0,0,255, &v[5], 0,0,0, &colors); Rast_write_colors(cellname, mapset, &colors); Rast_free_colors(&colors); } /* ---------------------------------------------------------------------- */ void setSinkWatershedColorTable(char* cellname) { struct Colors colors; const char *mapset; struct Range r; mapset = G_find_raster(cellname, ""); if (mapset == NULL) { G_fatal_error (_("Raster map <%s> not found"), cellname); } if (Rast_read_range(cellname, mapset, &r) == -1) { G_fatal_error(_("cannot read range")); } Rast_init_colors(&colors); Rast_make_random_colors(&colors, 1, r.max); Rast_write_colors(cellname, mapset, &colors); Rast_free_colors(&colors); } /* print the largest interm file that will be generated during r.terraflow */ void printMaxSortSize(long nodata_count) { char buf[BUFSIZ]; long long fillmaxsize = (long long)nrows*ncols*sizeof(waterWindowType); long long flowmaxsize = (long long)(nrows*ncols - nodata_count)*sizeof(sweepItem); long long maxneed = (fillmaxsize > flowmaxsize) ? fillmaxsize: flowmaxsize; maxneed = 2*maxneed; /* need 2*N to sort */ G_message( "total elements=%ld, nodata elements=%ld", (long)nrows*ncols, nodata_count); G_message( "largest temporary files: "); G_message( "\t\t FILL: %s [%d elements, %dB each]", formatNumber(buf, fillmaxsize), nrows * ncols, sizeof(waterWindowType)); G_message( "\t\t FLOW: %s [%ld elements, %dB each]", formatNumber(buf, flowmaxsize), (long)(nrows * ncols - nodata_count), sizeof(sweepItem)); G_message( "Will need at least %s space available in %s", formatNumber(buf, maxneed), /* need 2*N to sort */ getenv(STREAM_TMPDIR)); #ifdef HAVE_STATVFS_H fprintf(stderr, "Checking current space in %s: ", getenv(STREAM_TMPDIR)); struct statvfs statbuf; statvfs(getenv(STREAM_TMPDIR), &statbuf); float avail = statbuf.f_bsize*statbuf.f_bavail; fprintf(stderr, "available %ld blocks x %ldB = %.0fB", (long)statbuf.f_bavail, statbuf.f_bsize, avail); if (avail > maxneed) { fprintf(stderr, ". OK.\n"); } else { fprintf(stderr, ". Not enough space available.\n"); exit(EXIT_FAILURE); } #endif } /* ---------------------------------------------------------------------- */ int main(int argc, char *argv[]) { struct GModule *module; Rtimer rtTotal; char buf[BUFSIZ]; /* initialize GIS library */ G_gisinit(argv[0]); module = G_define_module(); #ifdef ELEV_SHORT module->description = _("Flow computation for massive grids (integer version)."); #endif #ifdef ELEV_FLOAT module->description = _("Flow computation for massive grids (float version)."); #endif G_add_keyword(_("raster")); G_add_keyword(_("hydrology")); /* read user options; fill in global */ opt = (userOptions*)malloc(sizeof(userOptions)); assert(opt); region = (struct Cell_head*)malloc(sizeof(struct Cell_head)); assert(region); parse_args(argc, argv); /* get the current region and dimensions */ G_get_set_window(region); check_args(); int nr = Rast_window_rows(); int nc = Rast_window_cols(); if ((nr > dimension_type_max) || (nc > dimension_type_max)) { G_fatal_error(_("[nrows=%d, ncols=%d] dimension_type overflow -- " "change dimension_type and recompile"), nr, nc); } else { nrows = (dimension_type)nr; ncols = (dimension_type)nc; } G_verbose_message( _("Region size is %d x %d"), nrows, ncols); /* check STREAM path (the place where intermediate STREAMs are placed) */ sprintf(buf, "%s=%s",STREAM_TMPDIR, opt->streamdir); /* don't pass an automatic variable; putenv() isn't guaranteed to make a copy */ putenv(G_store(buf)); if (getenv(STREAM_TMPDIR) == NULL) { fprintf(stderr, "%s:", STREAM_TMPDIR); G_fatal_error("not set"); } else { fprintf(stderr, "STREAM temporary files in %s ", getenv(STREAM_TMPDIR)); fprintf(stderr, "(THESE INTERMEDIATE STREAMS WILL NOT BE DELETED IN CASE OF ABNORMAL TERMINATION OF THE PROGRAM. TO SAVE SPACE PLEASE DELETE THESE FILES MANUALLY!)\n"); } if (opt->stats) { /* open the stats file */ stats = new statsRecorder(opt->stats); record_args(argc, argv); { char buf[BUFSIZ]; long grid_size = nrows * ncols; *stats << "region size = " << formatNumber(buf, grid_size) << " elts " << "(" << nrows << " rows x " << ncols << " cols)\n"; stats->flush(); } } /* set up STREAM memory manager */ size_t mm_size = (size_t) opt->mem << 20; /* opt->mem is in MB */ MM_manager.set_memory_limit(mm_size); if (opt->verbose) { MM_manager.warn_memory_limit(); } else { MM_manager.ignore_memory_limit(); } MM_manager.print_limit_mode(); /* initialize nodata */ nodataType::init(); if (stats) *stats << "internal nodata value: " << nodataType::ELEVATION_NODATA << endl; /* start timing -- after parse_args, which are interactive */ rt_start(rtTotal); #ifndef JUMP2FLOW /* read elevation into a stream */ AMI_STREAM *elstr=NULL; long nodata_count; elstr = cell2stream(opt->elev_grid, elevation_type_max, &nodata_count); /* print the largest interm file that will be generated */ printMaxSortSize(nodata_count); /* -------------------------------------------------- */ /* compute flow direction and filled elevation (and watersheds) */ AMI_STREAM *dirstr=NULL; AMI_STREAM *filledstr=NULL; AMI_STREAM *flowStream=NULL; AMI_STREAM *labeledWater = NULL; flowStream=computeFlowDirections(elstr, filledstr, dirstr, labeledWater); delete elstr; /* write streams to GRASS raster maps */ stream2_CELL(dirstr, nrows, ncols, opt->dir_grid); delete dirstr; #ifdef ELEV_SHORT stream2_CELL(filledstr, nrows, ncols, opt->filled_grid); #else stream2_CELL(filledstr, nrows, ncols, opt->filled_grid,true); #endif delete filledstr; stream2_CELL(labeledWater, nrows, ncols, labelElevTypePrintLabel(), opt->watershed_grid); setSinkWatershedColorTable(opt->watershed_grid); delete labeledWater; #else AMI_STREAM *flowStream; char path[GPATH_MAX]; sprintf(path, "%s/flowStream", streamdir->answer); flowStream = new AMI_STREAM(path); fprintf(stderr, "flowStream opened: len=%d\n", flowStream->stream_len()); fprintf(stderr, "jumping to flow accumulation computation\n"); #endif /* -------------------------------------------------- */ /* compute flow accumulation (and tci) */ AMI_STREAM *outstr=NULL; computeFlowAccumulation(flowStream, outstr); /* delete flowStream -- deleted inside */ /* write output stream to GRASS raster maps */ #ifdef OUTPUT_TCI stream2_FCELL(outstr, nrows, ncols, printAccumulation(), printTci(), opt->flowaccu_grid, opt->tci_grid); #else stream2_FCELL(outstr, nrows, ncols, printAccumulation(), opt->flowaccu_grid); #endif setFlowAccuColorTable(opt->flowaccu_grid); delete outstr; rt_stop(rtTotal); if (stats) { stats->recordTime("Total running time: ", rtTotal); stats->timestamp("end"); } G_done_msg(" "); /* free the globals */ free(region); free(opt); if (stats) delete stats; return 0; }