[GRASS-SVN] r56364 - in grass-addons/grass6/raster: . r.findtheriver
svn_grass at osgeo.org
svn_grass at osgeo.org
Wed May 22 11:27:41 PDT 2013
Author: selimnairb
Date: 2013-05-22 11:27:41 -0700 (Wed, 22 May 2013)
New Revision: 56364
Added:
grass-addons/grass6/raster/r.findtheriver/
grass-addons/grass6/raster/r.findtheriver/Makefile
grass-addons/grass6/raster/r.findtheriver/README
grass-addons/grass6/raster/r.findtheriver/README.md
grass-addons/grass6/raster/r.findtheriver/ReadMe-OSX.rtf
grass-addons/grass6/raster/r.findtheriver/description.html
grass-addons/grass6/raster/r.findtheriver/main.c
grass-addons/grass6/raster/r.findtheriver/point_list.c
grass-addons/grass6/raster/r.findtheriver/point_list.h
Log:
Added module r.findtheriver, useful for snapping streamflow gage coordinates to a stream network
Added: grass-addons/grass6/raster/r.findtheriver/Makefile
===================================================================
--- grass-addons/grass6/raster/r.findtheriver/Makefile (rev 0)
+++ grass-addons/grass6/raster/r.findtheriver/Makefile 2013-05-22 18:27:41 UTC (rev 56364)
@@ -0,0 +1,18 @@
+# fix this relative to include/
+# or use absolute path to the GRASS source code
+MODULE_TOPDIR = ../..
+GRASS_VERSION_STR = 6.4
+
+PGM = r.findtheriver
+
+COMPILE_FLAGS = -g -Wall -Werror-implicit-function-declaration -fno-common
+LIBES = $(GISLIB)
+DEPENDENCIES = $(GISDEP)
+
+include $(MODULE_TOPDIR)/include/Make/Module.make
+
+default: cmd
+
+deploy:
+ cp dist.*/bin/$(PGM) /Library/GRASS/$(GRASS_VERSION_STR)/Modules/bin
+ cp dist.*/docs/html/* /Library/GRASS/$(GRASS_VERSION_STR)/Modules/docs/html
\ No newline at end of file
Added: grass-addons/grass6/raster/r.findtheriver/README
===================================================================
--- grass-addons/grass6/raster/r.findtheriver/README (rev 0)
+++ grass-addons/grass6/raster/r.findtheriver/README 2013-05-22 18:27:41 UTC (rev 56364)
@@ -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/grass6/raster/r.findtheriver/README.md
===================================================================
--- grass-addons/grass6/raster/r.findtheriver/README.md (rev 0)
+++ grass-addons/grass6/raster/r.findtheriver/README.md 2013-05-22 18:27:41 UTC (rev 56364)
@@ -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/grass6/raster/r.findtheriver/ReadMe-OSX.rtf
===================================================================
--- grass-addons/grass6/raster/r.findtheriver/ReadMe-OSX.rtf (rev 0)
+++ grass-addons/grass6/raster/r.findtheriver/ReadMe-OSX.rtf 2013-05-22 18:27:41 UTC (rev 56364)
@@ -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/grass6/raster/r.findtheriver/description.html
===================================================================
--- grass-addons/grass6/raster/r.findtheriver/description.html (rev 0)
+++ grass-addons/grass6/raster/r.findtheriver/description.html 2013-05-22 18:27:41 UTC (rev 56364)
@@ -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><br>
+<EM><A HREF="http://grass.itc.it/devel/index.php#prog">GRASS Programmer's Manual</A></EM>
+
+
+<H2>AUTHOR</H2>
+
+Brian Miles - brian_miles at unc.edu
+
+<p><i>Last changed: $Date: 2013-02-09 10:43:00 -0500 (Sat, 09 Feb 2013) $</i>
Added: grass-addons/grass6/raster/r.findtheriver/main.c
===================================================================
--- grass-addons/grass6/raster/r.findtheriver/main.c (rev 0)
+++ grass-addons/grass6/raster/r.findtheriver/main.c 2013-05-22 18:27:41 UTC (rev 56364)
@@ -0,0 +1,433 @@
+/**
+\file main.c
+
+\brief r.findtheriver - Finds the nearest stream pixel to a
+coordinate pair using an upstream accumulating area map.
+
+(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
+
+\date May, 20 2013
+*/
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <assert.h>
+#include <math.h>
+
+#include <grass/config.h>
+#include <grass/gis.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, char *mapset, RASTER_MAP_TYPE dataType,
+ int windowSize, double threshold,
+ int nrows_less_one, int ncols_less_one,
+ int currRow, int currCol);
+
+/*
+ * global variables
+ */
+int quiet;
+
+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 */
+ 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 Flag *flag1; /* flags */
+ struct Option *optInput, *optWindow, *optThreshold, *optE, *optN, *optSep;
+ double E, N;
+ struct Cell_head window;
+ char *buff;
+ 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();
+ module->keywords = _("raster, keyword2, keyword3");
+ 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\n\t pixels. Must be an odd integer. If not supplied, window will be\n\t 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\n\t non-stream pixels. If not supplied, threshold will be inferred from \n\t 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\n\t in length.");
+
+ /* Define the different flags */
+ flag1 = G_define_flag();
+ flag1->key = 'q';
+ flag1->description = _("Quiet");
+
+ /* options and flags parser */
+ if (G_parser(argc, argv))
+ exit(EXIT_FAILURE);
+
+ if (G_get_window(&window) < 0) {
+ G_asprintf(&buff, _("Unable to read current window parameters"));
+ G_fatal_error(buff);
+ }
+
+ /* stores options and flags to variables */
+ strncpy(name, optInput->answer, GNAME_MAX);
+ quiet = (flag1->answer);
+
+ /* returns NULL if the map was not found in any mapset,
+ * mapset name otherwise */
+ mapset = G_find_cell2(name, "");
+ if (mapset == NULL) {
+ G_fatal_error(_("Raster map <%s> not found"), name);
+ }
+
+ /* Get raster metadata */
+ if (G_get_cellhd(name, mapset, &cellhd) < 0) {
+ G_fatal_error(_("Unable to read file header of <%s>"), name);
+ }
+
+ 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);
+ }
+ if ( !quiet ) {
+ G_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);
+ }
+ if ( !quiet ) {
+ G_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 = G_raster_map_type(name, mapset);
+
+ /* Open the raster - returns file descriptor (>0) */
+ if ((infd = G_open_cell_old(name, mapset)) < 0)
+ G_fatal_error(_("Unable to open raster map <%s>"), name);
+
+ G_get_set_window(&window);
+ nrows = G_window_rows();
+ ncols = G_window_cols();
+ nrows_less_one = nrows - 1;
+ ncols_less_one = ncols - 1;
+
+ rowIdx = (int)G_northing_to_row(N, &window);
+ colIdx = (int)G_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 = G_allocate_raster_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 */
+ if (G_get_raster_row(infd, tmpRow, currRow, data_type) < 0) {
+ G_fatal_error(_("Unable to read raster row %d"), currRow);
+ }
+
+ 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 = G_col_to_easting(nearestStreamPixel->col+0.5, &window);
+ double nearestNorthing = G_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 */
+ G_close_cell(infd);
+
+ exit(EXIT_SUCCESS);
+}
+
+/*
+ * function definitions
+ */
+PointList_t *find_stream_pixels_in_window(int fd, char *name, 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 = G_allocate_raster_buf(dataType);
+
+ /* Get value of central cell */
+ if (G_get_raster_row(fd, tmpRow, currRow, dataType) < 0) {
+ G_fatal_error(_("Unable to read raster row %d"), currRow);
+ }
+ 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 ( G_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
+
+ }
+
+ if ( !quiet ) {
+ G_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 */
+ if (G_get_raster_row(fd, tmpRow, row, dataType) < 0) {
+ G_fatal_error(_("Unable to read raster row %d"), row);
+ }
+ 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);
+ }
+ }
+ }
+ }
+ G_free(tmpRow);
+ return streamPixels;
+}
Added: grass-addons/grass6/raster/r.findtheriver/point_list.c
===================================================================
--- grass-addons/grass6/raster/r.findtheriver/point_list.c (rev 0)
+++ grass-addons/grass6/raster/r.findtheriver/point_list.c 2013-05-22 18:27:41 UTC (rev 56364)
@@ -0,0 +1,83 @@
+/****************************************************************************
+ *
+ * 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/grass6/raster/r.findtheriver/point_list.h
===================================================================
--- grass-addons/grass6/raster/r.findtheriver/point_list.h (rev 0)
+++ grass-addons/grass6/raster/r.findtheriver/point_list.h 2013-05-22 18:27:41 UTC (rev 56364)
@@ -0,0 +1,33 @@
+/****************************************************************************
+ *
+ * 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
More information about the grass-commit
mailing list