[GRASS-SVN] r60475 - in grass-addons/grass7/raster: . r.findtheriver
svn_grass at osgeo.org
svn_grass at osgeo.org
Sun May 25 06:12:57 PDT 2014
Author: hcho
Date: 2014-05-25 06:12:57 -0700 (Sun, 25 May 2014)
New Revision: 60475
Added:
grass-addons/grass7/raster/r.findtheriver/
grass-addons/grass7/raster/r.findtheriver/Makefile
grass-addons/grass7/raster/r.findtheriver/README
grass-addons/grass7/raster/r.findtheriver/README.md
grass-addons/grass7/raster/r.findtheriver/ReadMe-OSX.rtf
grass-addons/grass7/raster/r.findtheriver/main.c
grass-addons/grass7/raster/r.findtheriver/point_list.c
grass-addons/grass7/raster/r.findtheriver/point_list.h
grass-addons/grass7/raster/r.findtheriver/r.findtheriver.html
Log:
r.findtheriver: Updated to GRASS 7 and added to grass7 addon
Added: grass-addons/grass7/raster/r.findtheriver/Makefile
===================================================================
--- grass-addons/grass7/raster/r.findtheriver/Makefile (rev 0)
+++ grass-addons/grass7/raster/r.findtheriver/Makefile 2014-05-25 13:12:57 UTC (rev 60475)
@@ -0,0 +1,10 @@
+MODULE_TOPDIR = ../..
+
+PGM = r.findtheriver
+
+LIBES = $(RASTERLIB) $(GISLIB)
+DEPENDENCIES = $(RASTERDEP) $(GISDEP)
+
+include $(MODULE_TOPDIR)/include/Make/Module.make
+
+default: cmd
Added: grass-addons/grass7/raster/r.findtheriver/README
===================================================================
--- grass-addons/grass7/raster/r.findtheriver/README (rev 0)
+++ grass-addons/grass7/raster/r.findtheriver/README 2014-05-25 13:12:57 UTC (rev 60475)
@@ -0,0 +1,73 @@
+COPYRIGHT
+---------
+(C) 2013 by the University of North Carolina at Chapel Hill
+
+This program is free software under the GNU General Public License
+(>=v2). Read the file COPYING that comes with GRASS for details.
+
+
+AUTHOR
+------
+Brian Miles - brian_miles at unc.edu
+
+
+DESCRIPTION
+-----------
+r.findtheriver finds the nearest stream pixel to a coordinate pair
+using an upstream accumulating area (UAA) raster map. This is
+necessary because the coordinates for streamflow gages are often not
+perfectly registered to the topography represented by a digital
+elevation model (DEM) map. This presents a problem when trying to
+derive a watershed contributing area using r.water.outlet; if the
+streamflow gage does not fall on the stream as represented in the
+DEM, r.water.outlet can fail to derive the watershed area.
+
+The basic assumption is that the UAA for "stream" pixels will be much
+higher than for adjacent "non-stream" pixels.
+r.findtheriver attempts to "snap" the coordinates of the
+streamflow gage to the "true" stream location by first identifying
+stream pixels within a search window, and then selecting the stream
+pixel that is closest (cartesian distance) to the input gage
+coordinates. Stream pixels are identified by searching the UAA
+raster window for pixels that exceed a threshold. This is done by
+computing the log10 of the UAA value for the pixel corresponding to
+the gage coordinates and subtracting from it the log10 of each pixel
+in the window; for a given pixel if this difference is greater than
+the threshold, the pixel is deemed to be a stream pixel.
+
+r.findtheriver will automatically compute the window and threshold if
+they are not supplied by the user. The window is determined based on
+a THRESHOLD_DISTANCE / cell resolution of the UAA map. The threshold
+is determined by subtracting the log10 of the UAA value at the input
+gage coordinate from the log10 of the maximum UAA value of the map,
+and then rounding down to the nearest integer, in other words:
+threshold = floor( log(maxUAA) - log(gageUAA) ).
+
+The closest stream pixel is printed to standard output. If no stream
+pixels were found nothing is printed.
+
+
+BUILDING
+--------
+OS X
+
+First read William Kyngesburye's "Building Addon Modules for Mac OS X
+GRASS.app (see ReadMe-OSX.rtf). To install, you will generally do:
+
+make GRASS_HOME=./ GRASS_APP=/Applications/GRASS-6.4.app
+sudo make GRASS_HOME=./ GRASS_APP=/Applications/GRASS-6.4.app deploy
+
+Assuming you have unpacked r.findtheriver into your local copy of
+modbuild/module (see ReadMe-OSX.rtf) and the current directory is
+r.findtheriver.
+
+
+Linux/Unix
+
+First read: http://grasswiki.osgeo.org/wiki/Compile_and_Install#Addons
+
+To install, you will generally do:
+
+sudo make MODULE_TOPDIR=/usr/lib/grass64/ INST_DIR=/usr/lib/grass64/
+
+Note: you will have to change /usr/lib/grass64 to the location of your GRASS headers.
Added: grass-addons/grass7/raster/r.findtheriver/README.md
===================================================================
--- grass-addons/grass7/raster/r.findtheriver/README.md (rev 0)
+++ grass-addons/grass7/raster/r.findtheriver/README.md 2014-05-25 13:12:57 UTC (rev 60475)
@@ -0,0 +1,73 @@
+COPYRIGHT
+---------
+(C) 2013 by the University of North Carolina at Chapel Hill
+
+This program is free software under the GNU General Public License
+(>=v2). Read the file COPYING that comes with GRASS for details.
+
+
+AUTHOR
+------
+Brian Miles - brian_miles at unc.edu
+
+
+DESCRIPTION
+-----------
+r.findtheriver finds the nearest stream pixel to a coordinate pair
+using an upstream accumulating area (UAA) raster map. This is
+necessary because the coordinates for streamflow gages are often not
+perfectly registered to the topography represented by a digital
+elevation model (DEM) map. This presents a problem when trying to
+derive a watershed contributing area using r.water.outlet; if the
+streamflow gage does not fall on the stream as represented in the
+DEM, r.water.outlet can fail to derive the watershed area.
+
+The basic assumption is that the UAA for "stream" pixels will be much
+higher than for adjacent "non-stream" pixels.
+r.findtheriver attempts to "snap" the coordinates of the
+streamflow gage to the "true" stream location by first identifying
+stream pixels within a search window, and then selecting the stream
+pixel that is closest (cartesian distance) to the input gage
+coordinates. Stream pixels are identified by searching the UAA
+raster window for pixels that exceed a threshold. This is done by
+computing the log10 of the UAA value for the pixel corresponding to
+the gage coordinates and subtracting from it the log10 of each pixel
+in the window; for a given pixel if this difference is greater than
+the threshold, the pixel is deemed to be a stream pixel.
+
+r.findtheriver will automatically compute the window and threshold if
+they are not supplied by the user. The window is determined based on
+a THRESHOLD_DISTANCE / cell resolution of the UAA map. The threshold
+is determined by subtracting the log10 of the UAA value at the input
+gage coordinate from the log10 of the maximum UAA value of the map,
+and then rounding down to the nearest integer, in other words:
+threshold = floor( log(maxUAA) - log(gageUAA) ).
+
+The closest stream pixel is printed to standard output. If no stream
+pixels were found nothing is printed.
+
+
+BUILDING
+--------
+OS X
+
+First read William Kyngesburye's "Building Addon Modules for Mac OS X
+GRASS.app (see ReadMe-OSX.rtf). To install, you will generally do:
+
+make GRASS_HOME=./ GRASS_APP=/Applications/GRASS-6.4.app
+sudo make GRASS_HOME=./ GRASS_APP=/Applications/GRASS-6.4.app deploy
+
+Assuming you have unpacked r.findtheriver into your local copy of
+modbuild/module (see ReadMe-OSX.rtf) and the current directory is
+r.findtheriver.
+
+
+Linux/Unix
+
+First read: http://grasswiki.osgeo.org/wiki/Compile_and_Install#Addons
+
+To install, you will generally do:
+
+sudo make MODULE_TOPDIR=/usr/lib/grass64/ INST_DIR=/usr/lib/grass64/
+
+Note: you will have to change /usr/lib/grass64 to the location of your GRASS headers.
Added: grass-addons/grass7/raster/r.findtheriver/ReadMe-OSX.rtf
===================================================================
--- grass-addons/grass7/raster/r.findtheriver/ReadMe-OSX.rtf (rev 0)
+++ grass-addons/grass7/raster/r.findtheriver/ReadMe-OSX.rtf 2014-05-25 13:12:57 UTC (rev 60475)
@@ -0,0 +1,168 @@
+{\rtf1\ansi\ansicpg1252\cocoartf1187\cocoasubrtf370
+\cocoascreenfonts1{\fonttbl\f0\fswiss\fcharset0 Helvetica;\f1\fnil\fcharset0 Monaco;}
+{\colortbl;\red255\green255\blue255;}
+\margl1440\margr1440\vieww11220\viewh11220\viewkind0
+\pard\tx720\tx1440\tx2160\tx2880\tx3600\tx4320\tx5040\tx5760\tx6480\tx7200\tx7920\tx8640\pardirnatural\qc
+
+\f0\b\fs36 \cf0 Building Addon Modules for\
+Mac OS X GRASS.app\
+
+\b0\fs32 (without GEM)\
+\pard\tx720\tx1440\tx2160\tx2880\tx3600\tx4320\tx5040\tx5760\tx6480\tx7200\tx7920\tx8640\pardirnatural
+
+\fs28 \cf0 \
+GEM is a nice idea, but there are a few issues. Not many modules are setup for GEM (there are a couple in the new SVN repository, and some by the GEM author). It puts the resulting files in the GRASS binary folder and alters the GUI menus in the application package, and this is Not a Good Thing to do with Mac OS X applications.\
+\
+The method is similar to GEM in that it uses a minimal template of the GRASS source and uses the needed headers and libraries from GRASS.app. But there is less setup needed (for developers), and you don't need admin privileges to install and use the module, except in some more complex cases. Though a user does need to do a little setup, unlike with GEM.\
+\
+It's a bit rough right now. It doesn't have a module registry like GEM. It doesn't add modules to the GUI menu, though the menu is dynamically maintained by the GRASS.app startup.\
+\
+I haven't tried it with a module that also includes its own library.\
+\
+
+\b Note:
+\b0 if there is no makefile in the module source, it's probably just a script that you can simply copy to the bin folder (see \ul Use\ulnone section below).\
+\
+\
+
+\b Requirements\
+
+\b0 \
+- Xcode - at least version 2.3 (scripts shouldn't need this).\
+\
+- GRASS.app (you don't need the full GRASS source).\
+\
+- Module source. It should be compatible with the GRASS.app version you use.\
+\
+- This modbuild template.\
+\
+\
+
+\b Setup\
+
+\b0 \
+
+\f1\fs24 $VERSION
+\f0\fs28 is the major.minor GRASS version.\
+\
+- copy the modbuild folder from
+\f1\fs24 /Library/GRASS/$VERSION/
+\f0\fs28 to a location you have read/write access to. You could use your home
+\f1\fs24 /Library/GRASS/$VERSION/
+\f0\fs28 , for example.\
+\
+- Unpack the module source and put that source folder in the '
+\f1\fs24 module
+\f0\fs28 ' folder in the modbuild folder. If the module source is setup for GEM, you'll have to dig into the src/[section] subfolder to find the module's source folder (ie src/raster/module.source.folder). You should end up with your template folder like so:\
+\
+
+\f1\fs24 modbuild\
+ dist.some-platform-version\
+ include\
+ License.rtf\
+ module\
+ module_source_folder\
+ description.html (this might have the module name instead)\
+ main.c\
+ makefile\
+ (other_src.c)\
+ (other_header.h)\
+ ReadMe.rtf\
+ tools\
+
+\f0\fs28 \
+
+\b NOTE:
+\b0 make sure that there are no spaces in the path to the modbuild folder.\
+\
+- If the module uses extra data file(s) that are stored in the
+\f1\fs24 etc/
+\f0\fs28 folder (v.in.dwg does this), make sure it uses the new
+\f1\fs24 G_find_etc()
+\f0\fs28 function (for C modules) or the
+\f1\fs24 g.findetc
+\f0\fs28 module (for script modules) instead of the hardwired
+\f1\fs24 $GISBASE/etc/
+\f0\fs28 . This will mean either asking the author to make this change or attempting it yourself. See the developer documentation for usage of G_find_etc(), and the GRASS Help for usage of g.findetc.\
+\
+- In a Terminal,
+\f1\fs24 cd
+\f0\fs28 to the module source folder, then:\
+\
+
+\f1\fs24 $ make GRASS_HOME='/path/to/grass/modbuild' GRASS_APP='/path/to/GRASS.app'\
+
+\f0\fs28 \
+The modbuild path is the path to the
+\f1\fs24 modbuild
+\f0\fs28 folder, not the module source folder. Make sure to include the correct GRASS.app name in the GRASS_APP path (ie GRASS-6.3.app).\
+\
+\
+
+\b Use
+\b0 \
+\
+The module files will be in the
+\f1\fs24 dist
+\f0\fs28 folder in this template. A C source module will be in
+\f1\fs24 bin
+\f0\fs28 , script modules will be in
+\f1\fs24 scripts
+\f0\fs28 , and documentation will be in
+\f1\fs24 docs/html
+\f0\fs28 . GRASS.app has a couple default locations for addon modules that use the same configuration as the GRASS_ADDON_PATH environment variable:
+\f1\fs24 /Library/GRASS/$VERSION/Modules
+\f0\fs28 and
+\f1\fs24 ~/Library/GRASS/$VERSION/Modules
+\f0\fs28 . The version is the Major.minor GRASS version. Modules that should be available to all users should be installed in
+\f1\fs24 /Library
+\f0\fs28 (requires admin privileges), and any modules in
+\f1\fs24 ~/Library
+\f0\fs28 will only be available to the user that installs them. The user GRASS folder will be created the first time GRASS is run.\
+\
+
+\i Installation is not automatic.
+\i0 You must manually copy the files because the install targets in the GRASS make system only exist in the GRASS source.\
+\
+- Both
+\b scripts
+\b0 and
+\b binary
+\b0 executables can be run from the same location, so install them in the
+\f1\fs24 bin/
+\f0\fs28 folder.\
+- Any
+\b libraries
+\b0 should go in the
+\f1\fs24 lib/
+\f0\fs28 folder.\
+-
+\b Documentation
+\b0 goes in the
+\f1\fs24 docs/html/
+\f0\fs28 folder.\
+-
+\b Data
+\b0 files should go into
+\f1\fs24 etc/
+\f0\fs28 .\
+\pard\tx720\tx1440\tx2160\tx2880\tx3600\tx4320\tx5040\tx5760\tx6480\tx7200\tx7920\tx8640
+\cf0 \
+\
+\pard\tx720\tx1440\tx2160\tx2880\tx3600\tx4320\tx5040\tx5760\tx6480\tx7200\tx7920\tx8640
+
+\fs24 \cf0 \'a9 2006-2007 by the GRASS Development Team\
+This program is free software under the GNU General Public License (>=v2).\
+\pard\tx720\tx1440\tx2160\tx2880\tx3600\tx4320\tx5040\tx5760\tx6480\tx7200\tx7920\tx8640
+
+\fs28 \cf0 \
+\
+- William Kyngesburye\
+\pard\tx720\tx1440\tx2160\tx2880\tx3600\tx4320\tx5040\tx5760\tx6480\tx7200\tx7920\tx8640
+{\field{\*\fldinst{HYPERLINK "mailto:kyngchaos at kyngchaos.com"}}{\fldrslt \cf0 kyngchaos at kyngchaos.com}}
+\fs24 \
+\pard\tx720\tx1440\tx2160\tx2880\tx3600\tx4320\tx5040\tx5760\tx6480\tx7200\tx7920\tx8640
+{\field{\*\fldinst{HYPERLINK "http://www.kyngchaos.com"}}{\fldrslt
+\fs28 \cf0 http://www.kyngchaos.com/}}
+\fs28 \
+}
\ No newline at end of file
Added: grass-addons/grass7/raster/r.findtheriver/main.c
===================================================================
--- grass-addons/grass7/raster/r.findtheriver/main.c (rev 0)
+++ grass-addons/grass7/raster/r.findtheriver/main.c 2014-05-25 13:12:57 UTC (rev 60475)
@@ -0,0 +1,480 @@
+
+/****************************************************************************
+ *
+ * MODULE: r.findtheriver
+ *
+ * AUTHOR(S): Brian Miles <brian_miles unc.edu> (May 20, 2013)
+ * Update to GRASS 7 and returning the input cell as needed
+ * by Huidae Cho <grass4u gmail.com>
+ *
+ * PURPOSE: Finds the nearest stream pixel to a coordinate pair using an
+ * upstream accumulating area map.
+ *
+ * COPYRIGHT: (C) 2013 by the University of North Carolina at Chapel Hill
+ *
+ * This program is free software under the GNU General Public
+ * License (>=v2). Read the file COPYING that comes with GRASS
+ * for details.
+ *
+ *****************************************************************************/
+
+#include <stdlib.h>
+#include <string.h>
+#include <assert.h>
+#include <math.h>
+#include <grass/raster.h>
+#include <grass/glocale.h>
+#include "point_list.h"
+
+#define THRESHOLD_DISTANCE 100
+#define DEFAULT_COORD_SEP ' '
+
+/*
+ * global function declaration
+ */
+
+ /**
+\brief Find stream pixels in a given window of the UAA raster
+
+\param[in] fd UAA raster file descriptor
+\param[in] name Name of UAA raster map
+\param[in] mapset Name of the GRASS mapset from which UAA raster should be read
+\param[in] dataType Data type of UAA raster
+\param[in] windowSize Size of the search window, must be an odd integer >= 3
+\param[in] threshold to use for distinguishing stream pixels by comparing log(UAA) values
+\param[in] nrows_less_one Number of rows in the current window, minus 1
+\param[in] ncols_less_one Number of columns in the current window, minus 1
+\param[in] currRow, Row of pixel that will be compared to pixels in the window
+\param[in] currCol, Column of pixel that will be compared to pixels in the window
+
+\return Pointer to PointList_t representing a list of stream pixels found in the window
+*/
+PointList_t *find_stream_pixels_in_window(int fd, char *name,
+ const char *mapset,
+ RASTER_MAP_TYPE dataType,
+ int windowSize, double threshold,
+ int nrows_less_one,
+ int ncols_less_one, int currRow,
+ int currCol);
+
+int main(int argc, char *argv[])
+{
+ struct Cell_head cellhd; /* it stores region information,
+ and header information of rasters */
+ char name[GNAME_MAX]; /* input raster name */
+
+ const char *mapset; /* mapset name */
+
+ int nrows, ncols;
+
+ int rowIdx, colIdx, nrows_less_one, ncols_less_one;
+
+ int infd; /* file descriptor */
+
+ RASTER_MAP_TYPE data_type; /* type of the map (CELL/DCELL/...) */
+
+ struct GModule *module; /* GRASS module for parsing arguments */
+
+ struct Option *optInput, *optWindow, *optThreshold, *optE, *optN, *optSep;
+
+ double E, N;
+
+ struct Cell_head window;
+
+ int windowSize;
+
+ double threshold;
+
+ size_t SEP_SIZE = 2;
+
+ char *sep = (char *)G_malloc(SEP_SIZE * sizeof(char));
+
+ /* initialize GIS environment */
+ G_gisinit(argv[0]); /* reads grass env, stores program name to G_program_name() */
+
+ /* initialize module */
+ module = G_define_module();
+ G_add_keyword(_("raster"));
+ G_add_keyword(_("hydrology"));
+ module->description =
+ _("Find the stream pixel nearest the input coordinate");
+
+ /* Define command options */
+ optInput = G_define_option();
+ optInput->key = "accumulation";
+ optInput->type = TYPE_STRING;
+ optInput->required = YES;
+ optInput->gisprompt = "old,cell,raster";
+ optInput->description =
+ _("Name of input upstream accumulation area raster map");
+
+ optWindow = G_define_option();
+ optWindow->key = "window";
+ optWindow->type = TYPE_INTEGER;
+ optWindow->key_desc = "x";
+ optWindow->multiple = NO;
+ optWindow->required = NO;
+ optWindow->description =
+ _
+ ("The size of the window, in pixels, to search in for stream pixels. Must be an odd integer. If not supplied, window will be inferred based on raster resolution.");
+
+ optThreshold = G_define_option();
+ optThreshold->key = "threshold";
+ optThreshold->type = TYPE_DOUBLE;
+ optThreshold->key_desc = "x";
+ optThreshold->multiple = NO;
+ optThreshold->required = NO;
+ optThreshold->description =
+ _
+ ("The threshold for distinguishing log(UAA) values of stream and non-stream pixels. If not supplied, threshold will be inferred from minimum and maximum raster values.");
+
+ optE = G_define_option();
+ optE->key = "easting";
+ optE->type = TYPE_STRING;
+ optE->key_desc = "x";
+ optE->multiple = NO;
+ optE->required = YES;
+ optE->description = _("The map E grid coordinates");
+
+ optN = G_define_option();
+ optN->key = "northing";
+ optN->type = TYPE_STRING;
+ optN->key_desc = "y";
+ optN->multiple = NO;
+ optN->required = YES;
+ optN->description = _("The map N grid coordinates");
+
+ optSep = G_define_option();
+ optSep->key = "separator";
+ optSep->type = TYPE_STRING;
+ optSep->key_desc = "y";
+ optSep->multiple = NO;
+ optSep->required = NO;
+ optSep->description =
+ _
+ ("Coordinate separator. Defaults to ' '. Must be 1 character in length.");
+
+ /* options and flags parser */
+ if (G_parser(argc, argv))
+ exit(EXIT_FAILURE);
+
+ G_get_window(&window);
+
+ /* stores options and flags to variables */
+ strncpy(name, optInput->answer, GNAME_MAX);
+
+ /* returns NULL if the map was not found in any mapset,
+ * mapset name otherwise */
+ mapset = G_find_raster(name, "");
+ if (mapset == NULL) {
+ G_fatal_error(_("Raster map <%s> not found"), name);
+ }
+
+ /* Get raster metadata */
+ Rast_get_cellhd(name, mapset, &cellhd);
+
+ if (NULL != optWindow->answer) {
+ windowSize = atoi(optWindow->answer);
+ }
+ else {
+ /* Determine window size */
+ double cellRes = (cellhd.ew_res + cellhd.ns_res) / 2;
+
+ windowSize = THRESHOLD_DISTANCE / cellRes;
+ if (!(windowSize & 1))
+ windowSize++;
+ }
+ if ((windowSize < 2) || !(windowSize & 1)) {
+ G_warning(_
+ ("Invalid window size %s. Window size must be an odd integer >= 3\n"),
+ optWindow->answer);
+ G_usage();
+ exit(EXIT_FAILURE);
+ }
+ G_verbose_message(_("Stream search window size %d\n"), windowSize);
+
+ if (NULL != optThreshold->answer) {
+ threshold = atof(optThreshold->answer);
+ }
+ else {
+ /* Automatically determine the threshold */
+ threshold = -1.0;
+ }
+ if (threshold != -1.0 && threshold <= 0.0) {
+ G_warning(_("Invalid threshold %s. Threshold must be > 0.0\n"),
+ optThreshold->answer);
+ G_usage();
+ exit(EXIT_FAILURE);
+ }
+
+ if (!G_scan_easting(*optE->answers, &E, G_projection())) {
+ G_warning(_("Illegal east coordinate <%s>\n"), optE->answer);
+ G_usage();
+ exit(EXIT_FAILURE);
+ }
+ if (!G_scan_northing(*optN->answers, &N, G_projection())) {
+ G_warning(_("Illegal north coordinate <%s>\n"), optN->answer);
+ G_usage();
+ exit(EXIT_FAILURE);
+ }
+ G_verbose_message(_("Input coordinates, easting %f, northing %f\n"), E, N);
+
+ if (NULL == optSep->answer) {
+ sep[0] = DEFAULT_COORD_SEP;
+ }
+ else {
+ if (strlen(optSep->answer) != 1) {
+ G_warning(_("Separator must be 1 character in length\n"));
+ G_usage();
+ exit(EXIT_FAILURE);
+ }
+ strncpy(sep, optSep->answer, SEP_SIZE);
+ }
+
+ /* Determine the inputmap type (CELL/FCELL/DCELL) */
+ data_type = Rast_map_type(name, mapset);
+
+ /* Open the raster - returns file descriptor (>0) */
+ if ((infd = Rast_open_old(name, mapset)) < 0)
+ G_fatal_error(_("Unable to open raster map <%s>"), name);
+
+ G_get_set_window(&window);
+ nrows = Rast_window_rows();
+ ncols = Rast_window_cols();
+ nrows_less_one = nrows - 1;
+ ncols_less_one = ncols - 1;
+
+ rowIdx = (int)Rast_northing_to_row(N, &window);
+ colIdx = (int)Rast_easting_to_col(E, &window);
+
+ //double currNearestE, prevNearestE, currNearestN, prevNearestN;
+ PointList_t *streamPixels =
+ find_stream_pixels_in_window(infd, (char *)&name, mapset, data_type,
+ windowSize, threshold,
+ nrows_less_one, ncols_less_one,
+ rowIdx, colIdx);
+
+#ifdef DEBUG
+ fprintf(stderr, "Stream pixels: ");
+ print_list(stderr, streamPixels, " ");
+ fprintf(stderr, "\n");
+ fflush(stderr);
+#endif
+ PointList_t *nearestStreamPixel =
+ find_nearest_point(streamPixels, colIdx, rowIdx);
+
+ if (NULL != nearestStreamPixel) {
+
+#ifdef DEBUG
+ double nearestValue;
+
+ void *tmpRow = Rast_allocate_buf(data_type);
+
+ int currCol = nearestStreamPixel->col;
+
+ int currRow = nearestStreamPixel->row;
+
+ fprintf(stderr, "Nearest pixel col: %d, row: %d\n", currCol, currRow);
+ fflush(stderr);
+
+ /* Get value of central cell */
+ Rast_get_row(infd, tmpRow, currRow, data_type);
+
+ switch (data_type) {
+ case FCELL_TYPE:
+ nearestValue = (double)((FCELL *) tmpRow)[currCol];
+ break;
+ case DCELL_TYPE:
+ nearestValue = (double)((DCELL *) tmpRow)[currCol];
+ break;
+ default:
+ nearestValue = (double)((CELL *) tmpRow)[currCol];
+ break;
+ }
+
+ fprintf(stderr, "Nearest stream pixel UAA value: %f\n", nearestValue);
+ fflush(stderr);
+ G_free(tmpRow);
+#endif
+
+ /* Get center of each column */
+ double nearestEasting =
+ Rast_col_to_easting(nearestStreamPixel->col + 0.5, &window);
+ double nearestNorthing =
+ Rast_row_to_northing(nearestStreamPixel->row + 0.5, &window);
+
+ /* Print snapped coordinates */
+ fprintf(stdout, _("%f%s%f\n"), nearestEasting, sep, nearestNorthing);
+ fflush(stdout);
+
+ }
+
+ /* Clean up */
+ destroy_list(streamPixels);
+
+ /* closing raster maps */
+ Rast_close(infd);
+
+ exit(EXIT_SUCCESS);
+}
+
+/*
+ * function definitions
+ */
+PointList_t *find_stream_pixels_in_window(int fd, char *name,
+ const char *mapset,
+ RASTER_MAP_TYPE dataType,
+ int windowSize, double threshold,
+ int nrows_less_one,
+ int ncols_less_one, int currRow,
+ int currCol)
+{
+ assert(DCELL_TYPE == dataType);
+ assert(windowSize & 1);
+
+ PointList_t *streamPixels = NULL;
+
+ DCELL centralValue, tmpValue;
+
+ double logCentralValue, logTmpValue;
+
+ void *tmpRow = Rast_allocate_buf(dataType);
+
+ /* Get value of central cell */
+ Rast_get_row(fd, tmpRow, currRow, dataType);
+
+ switch (dataType) {
+ case CELL_TYPE:
+ centralValue = (double)((CELL *) tmpRow)[currCol];
+ break;
+ case FCELL_TYPE:
+ centralValue = (double)((FCELL *) tmpRow)[currCol];
+ break;
+ case DCELL_TYPE:
+ centralValue = (double)((DCELL *) tmpRow)[currCol];
+ break;
+ }
+ if (centralValue <= 0)
+ centralValue = 1;
+ logCentralValue = log10(centralValue);
+#ifdef DEBUG
+ fprintf(stderr, "logCentralValue: %f\n", logCentralValue);
+ fflush(stderr);
+#endif
+
+ /* Determine threshold if need be */
+ if (-1.0 == threshold) {
+ struct FPRange *range =
+ (struct FPRange *)G_malloc(sizeof(struct FPRange));
+ if (Rast_read_fp_range(name, mapset, range) < 0) {
+ G_fatal_error(_("Unable to determine range of raster map <%s>"),
+ name);
+ }
+ double max = range->max;
+
+ G_free(range); // eggs or is it chicken?
+ if (max == centralValue)
+ return streamPixels;
+
+ double logMax = log10(max);
+
+ threshold = floor(logMax - logCentralValue);
+ if (threshold <= 0.0) {
+ threshold = 1.0;
+ }
+ else if (threshold > 2.0) {
+ threshold = 2.0;
+ }
+
+#ifdef DEBUG
+ fprintf(stderr, "logMax: %f\n", logMax);
+ fflush(stderr);
+#endif
+
+ }
+
+ G_verbose_message(_("Threshold: %f\n"), threshold);
+
+ /* Define window bounds */
+ int windowOffset = (windowSize - 1) / 2;
+
+ int minCol = currCol - windowOffset;
+
+ if (minCol < 0)
+ minCol = 0;
+ int maxCol = currCol + windowOffset;
+
+ if (maxCol > ncols_less_one)
+ maxCol = ncols_less_one;
+ int minRow = currRow - windowOffset;
+
+ if (minRow < 0)
+ minRow = 0;
+ int maxRow = currRow + windowOffset;
+
+ if (maxRow > nrows_less_one)
+ maxRow = nrows_less_one;
+#ifdef DEBUG
+ fprintf(stderr, "currCol: %d, currRow: %d\n", currCol, currRow);
+ fprintf(stderr, "min. col: %d, max. col: %d\n", minCol, maxCol);
+ fprintf(stderr, "min. row: %d, max. row: %d\n", minRow, maxRow);
+ fflush(stderr);
+#endif
+
+ /* Search for stream pixels within the window */
+ int row, col;
+
+ for (row = minRow; row <= maxRow; row++) {
+#ifdef DEBUG
+ fprintf(stderr, "row: %d\n", row);
+ fflush(stderr);
+#endif
+ /* Get the current row */
+ Rast_get_row(fd, tmpRow, row, dataType);
+
+ for (col = minCol; col <= maxCol; col++) {
+#ifdef DEBUG
+ fprintf(stderr, "\tcol: %d\n", col);
+ fflush(stderr);
+#endif
+ switch (dataType) {
+ case CELL_TYPE:
+ tmpValue = (double)((CELL *) tmpRow)[col];
+ break;
+ case FCELL_TYPE:
+ tmpValue = (double)((FCELL *) tmpRow)[col];
+ break;
+ case DCELL_TYPE:
+ tmpValue = (double)((DCELL *) tmpRow)[col];
+ break;
+ }
+ logTmpValue = log10(tmpValue);
+ /* Test for nearby pixels that are stream pixels when compared to the central pixel */
+#ifdef DEBUG
+ fprintf(stderr, "\t\tlogTmpValue: %f, logCentralValue: %f\n",
+ logTmpValue, logCentralValue);
+ fflush(stderr);
+#endif
+ if ((logTmpValue - logCentralValue) > threshold) {
+ /* Add to list of stream pixels */
+ if (NULL == streamPixels) {
+ streamPixels = create_list(col, row);
+ }
+ else {
+ append_point(streamPixels, col, row);
+ }
+ }
+ else if ((logCentralValue - logTmpValue) > threshold) {
+ /* Add to list of stream pixels */
+ if (NULL == streamPixels) {
+ streamPixels = create_list(currCol, currRow);
+ }
+ else {
+ append_point(streamPixels, currCol, currRow);
+ }
+ }
+ }
+ }
+ G_free(tmpRow);
+ return streamPixels;
+}
Added: grass-addons/grass7/raster/r.findtheriver/point_list.c
===================================================================
--- grass-addons/grass7/raster/r.findtheriver/point_list.c (rev 0)
+++ grass-addons/grass7/raster/r.findtheriver/point_list.c 2014-05-25 13:12:57 UTC (rev 60475)
@@ -0,0 +1,97 @@
+
+/****************************************************************************
+ *
+ * MODULE: r.findtheriver
+ * AUTHOR(S): Brian Miles - brian_miles at unc.edu
+ * with hints from: Glynn Clements - glynn gclements.plus.com
+ * PURPOSE: Finds the nearest stream pixel to a coordinate pair using
+ * an upstream accumulating area map.
+ *
+ * COPYRIGHT: (C) 2013 by the University of North Carolina at Chapel Hill
+ *
+ * This program is free software under the GNU General Public
+ * License (>=v2). Read the file COPYING that comes with GRASS
+ * for details.
+ *
+ *****************************************************************************/
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <assert.h>
+#include <math.h>
+
+#include <grass/config.h>
+#include <grass/gis.h>
+
+#include "point_list.h"
+
+PointList_t *create_list(int col, int row)
+{
+ PointList_t *head = (PointList_t *) G_malloc(sizeof(PointList_t));
+
+ head->col = col;
+ head->row = row;
+ head->next = NULL;
+ return head;
+}
+
+PointList_t *append_point(PointList_t * const head, int col, int row)
+{
+ assert(NULL != head);
+
+ PointList_t *tmp = head;
+
+ while (NULL != tmp->next)
+ tmp = tmp->next;
+ tmp->next = (PointList_t *) G_malloc(sizeof(PointList_t));
+ tmp->next->col = col;
+ tmp->next->row = row;
+ tmp->next->next = NULL;
+ return tmp->next;
+}
+
+void destroy_list(PointList_t * head)
+{
+ if (NULL != head) {
+ PointList_t *curr = head;
+
+ PointList_t *next = curr->next;
+
+ while (NULL != next) {
+ free(curr);
+ curr = next;
+ next = curr->next;
+ }
+ G_free(curr);
+ }
+}
+
+PointList_t *find_nearest_point(PointList_t * const head, int col, int row)
+{
+
+ PointList_t *nearest = NULL;
+
+ double tmpDistance, minDistance = HUGE_VAL;
+
+ PointList_t *curr = head;
+
+ while (NULL != curr) {
+ tmpDistance = sqrt(pow(col - curr->col, 2) + pow(row - curr->row, 2));
+ if (tmpDistance < minDistance) {
+ minDistance = tmpDistance;
+ nearest = curr;
+ }
+ curr = curr->next;
+ }
+ return nearest;
+}
+
+void print_list(FILE * fp, PointList_t * const head, const char *const sep)
+{
+ PointList_t *curr = head;
+
+ while (NULL != curr) {
+ fprintf(fp, "%d%s%d\n", curr->col, sep, curr->row);
+ curr = curr->next;
+ }
+}
Added: grass-addons/grass7/raster/r.findtheriver/point_list.h
===================================================================
--- grass-addons/grass7/raster/r.findtheriver/point_list.h (rev 0)
+++ grass-addons/grass7/raster/r.findtheriver/point_list.h 2014-05-25 13:12:57 UTC (rev 60475)
@@ -0,0 +1,39 @@
+
+/****************************************************************************
+ *
+ * MODULE: r.findtheriver
+ * AUTHOR(S): Brian Miles - brian_miles at unc.edu
+ * with hints from: Glynn Clements - glynn gclements.plus.com
+ * PURPOSE: Finds the nearest stream pixel to a coordinate pair using
+ * an upstream accumulating area map.
+ *
+ * COPYRIGHT: (C) 2013 by the University of North Carolina at Chapel Hill
+ *
+ * This program is free software under the GNU General Public
+ * License (>=v2). Read the file COPYING that comes with GRASS
+ * for details.
+ *
+ *****************************************************************************/
+
+#ifndef __POINT_LIST_H__
+#define __POINT_LIST_H__
+
+#include <stdio.h>
+
+typedef struct PointNode
+{
+ int col;
+ int row;
+ struct PointNode *next;
+} PointList_t;
+
+PointList_t *create_list(int col, int row);
+
+PointList_t *append_point(PointList_t * const head, int col, int row);
+
+void destroy_list(PointList_t * head);
+
+PointList_t *find_nearest_point(PointList_t * const head, int col, int row);
+
+void print_list(FILE * fp, PointList_t * const head, const char *const sep);
+#endif
Added: grass-addons/grass7/raster/r.findtheriver/r.findtheriver.html
===================================================================
--- grass-addons/grass7/raster/r.findtheriver/r.findtheriver.html (rev 0)
+++ grass-addons/grass7/raster/r.findtheriver/r.findtheriver.html 2014-05-25 13:12:57 UTC (rev 60475)
@@ -0,0 +1,46 @@
+<h2>DESCRIPTION</h2>
+
+<p><em>r.findtheriver</em> finds the nearest stream pixel to a coordinate pair
+using an upstream accumulating area (UAA) raster map. This is
+necessary because the coordinates for streamflow gages are often not
+perfectly registered to the topography represented by a digital
+elevation model (DEM) map. This presents a problem when trying to
+derive a watershed contributing area using r.water.outlet; if the
+streamflow gage does not fall on the stream as represented in the
+DEM, r.water.outlet can fail to derive the watershed area.</p>
+
+<p>The basic assumption is that the UAA for "stream" pixels will be much
+higher than for adjacent "non-stream" pixels.
+r.findtheriver attempts to "snap" the coordinates of the
+streamflow gage to the "true" stream location by first identifying
+stream pixels within a search window, and then selecting the stream
+pixel that is closest (cartesian distance) to the input gage
+coordinates. Stream pixels are identified by searching the UAA
+raster window for pixels that exceed a threshold. This is done by
+computing the log10 of the UAA value for the pixel corresponding to
+the gage coordinates and subtracting from it the log10 of each pixel
+in the window; for a given pixel if this difference is greater than
+the threshold, the pixel is deemed to be a stream pixel.</p>
+
+<p>r.findtheriver will automatically compute the window and threshold if
+they are not supplied by the user. The window is determined based on
+a THRESHOLD_DISTANCE / cell resolution of the UAA map. The threshold
+is determined by subtracting the log10 of the UAA value at the input
+gage coordinate from the log10 of the maximum UAA value of the map,
+and then rounding down to the nearest integer, in other words:
+threshold = floor( log(maxUAA) - log(gageUAA) ).</p>
+
+<p>The closest stream pixel is printed to standard output. If no stream
+pixels were found nothing is printed.</p>
+
+<h2>SEE ALSO</h2>
+
+<em><a href="r.water.outlet.html">r.water.outlet</a></em>
+
+
+<h2>AUTHORS</h2>
+
+<a href="mailto:brian_miles at unc edu">Brian Miles</a><br>
+Update to GRASS 7 by <a href="mailto:grass4u at gmail com">Huidae Cho</a>
+
+<p><i>Last changed: $Date: 2013-07-11 18:46:52 -0400 (Thu, 11 Jul 2013) $</i>
More information about the grass-commit
mailing list