[GRASS-SVN] r51665 - in grass-addons/grass7/vector: . v.surf.icw
svn_grass at osgeo.org
svn_grass at osgeo.org
Tue May 22 06:40:42 EDT 2012
Author: hamish
Date: 2012-05-22 03:40:41 -0700 (Tue, 22 May 2012)
New Revision: 51665
Added:
grass-addons/grass7/vector/v.surf.icw/
grass-addons/grass7/vector/v.surf.icw/v.surf.icw.html
grass-addons/grass7/vector/v.surf.icw/v.surf.icw.py
Removed:
grass-addons/grass7/vector/v.surf.icw/description.html
grass-addons/grass7/vector/v.surf.icw/qgis-toolbox/
grass-addons/grass7/vector/v.surf.icw/v.surf.icw
Log:
part 1 of python port (still a bit more to do before it's working)
Deleted: grass-addons/grass7/vector/v.surf.icw/description.html
===================================================================
--- grass-addons/grass6/vector/v.surf.icw/description.html 2011-12-07 09:28:27 UTC (rev 49599)
+++ grass-addons/grass7/vector/v.surf.icw/description.html 2012-05-22 10:40:41 UTC (rev 51665)
@@ -1,115 +0,0 @@
-<h2>DESCRIPTION</h2>
-
-Inverse cost weighting is like inverse distance weighted (IDW) interpolation,
-but uses <i>cost</i> instead of shortest Euclidean distance. In this way
-solid barriers and molasses zones may be correctly taken into account.
-<p>
-Input data points do not need to have direct line of sight to each other.
-This is interpolation "as the fish swims", not "as the crow flies", and can
-see around headlands or across land-bridges without polluting over barriers
-which the natural value (or flightless bird) can not cross.
-<p>
-It was initially written to interpolate water chemistry in two parallel
-arms of a fjord, but may just as well be used for population abundance
-constrained by topography, or in studies of archeologic technology transfer.
-
-
-<h2>NOTES</h2>
-
-In the simplest case, the cost map will just be a mask raster with values
-of 1 in areas to interpolate, and NULL in impenetrable areas.
-Fancier cost maps can be used, for example, to make it more expensive for
-a measured pollutant to diffuse upstream in an estuary, or to make it more
-expensive for a stone tool technology to cross waterways.
-<p>
-Since generating cost maps can take a long time, it is recommended to keep
-the raster region relatively small and limit the number of starting points
-to less than 100.
-<p>
-Higher values of <b>friction</b> will help limit unconstrained boundary
-effects at the edges of your coverage, but will incur more of a stepped
-transition between sites.
-<p>
-The <b>post_mask</b>, if given, is applied after the interpolation is
-complete. A common use for that might be to only present data within
-a certain distance (thus confidence) of an actual sampling station.
-In that case the <em>r.cost</em> module can be used to create the mask.
-<p>
-This module writes lots of temporary files and so can be rather disk I/O
-intensive. If you are running it many times in a big loop you may want to
-set up a RAM-disk to put the mapset in (on UNIX symlinking back into the
-location is ok), or adjust the disk-cache flushing timer to be slightly
-longer than one iteration of the script.
-<br>To do this on Linux you can tune the
- <tt>/proc/sys/vm/dirty_expire_centisecs</tt> kernel control. The default
-is to keep files in memory a maximum of 30 seconds before writing them to
-disk, although if you are short on free RAM the kernel may write to disk
-long before that timeout is reached.
-<!-- Linux I/O tuning hints:
- http://www.westnet.com/~gsmith/content/linux-pdflush.htm
- http://tldp.org/LDP/sag/html/buffer-cache.html
- -->
-
-
-<h2>EXAMPLE</h2>
-
-In this example we'll generate some fake island barriers from the
-standard North Carolina GRASS dataset, then interpolate a continuous
-variable between given point stations. We'll use rainfall, but for the
-purposes of the exercise pretend it is a nutrient concentration in our
-fake wayerway. Point data stations outside of the current region or
-in NULL areas of the <i>cost_map</i> will be ignored.
-
-<div class="code"><pre>
-# set up a fake sea with islands:
-g.region n=276000 s=144500 w=122000 e=338500 res=500
-r.mapcalc "pseudo_elev = elev_state_500m - 1100"
-r.colors pseudo_elev color=etopo2
-r.mapcalc "navigable_mask = if(pseudo_elev < 0, 1, null())"
-
-# pick a data column from the points vector:
-v.info -c precip_30ynormals
-
-# run the interpolation:
-v.surf.icw input=precip_30ynormals column=annual output=annual_interp.3 \
- cost_map=navigable_mask friction=3 --quiet
-
-# equalize colors to show maximum detail:
-r.colors -e annual_interp.3 color=bcyr
-
-# display results in an Xmonitor:
-d.erase black
-d.rast -o annual_interp.3
-d.vect precip_30ynormals fcolor=red icon=basic/circle
-d.legend annual_interp.3 color=white at=48.4,94.8,3.4,6.0
-</pre></div>
-
-
-<!-- todo
-<h2>REFERENCES</h2>
-
- * Hamish's MfE reports / published confernce proceedings?
- * Ben D's journal article (Deutsch)
--->
-
-
-<h2>SEE ALSO</h2>
-
-<em>
-<a href="v.surf.idw.html">v.surf.idw</a><br>
-<a href="v.surf.rst.html">v.surf.rst</a><br>
-<a href="v.surf.bspline.html">v.surf.bspline</a><br>
-<a href="r.cost.html">r.cost</a><br>
-<a href="r.surf.idw.html">r.surf.idw</a><br>
-<a href="r.surf.idw2.html">r.surf.idw2</a><br>
-</em>
-
-
-<h2>AUTHOR</h2>
-
-Hamish Bowman<br>
-<i>Department of Marine Science,<br>
-Dunedin, New Zealand</i>
-
-<p>
-<i>Last changed: $Date$</i>
Deleted: grass-addons/grass7/vector/v.surf.icw/v.surf.icw
===================================================================
--- grass-addons/grass6/vector/v.surf.icw/v.surf.icw 2011-12-07 09:28:27 UTC (rev 49599)
+++ grass-addons/grass7/vector/v.surf.icw/v.surf.icw 2012-05-22 10:40:41 UTC (rev 51665)
@@ -1,407 +0,0 @@
-#!/bin/bash
-#############################################################################
-#
-# MODULE: v.surf.icw
-# version $Id$
-#
-# AUTHOR: M. Hamish Bowman, Dunedin, New Zealand
-# Originally written aboard the NZ DoC ship M/V Renown,
-# Bligh Sound, Fiordland National Park, November 2003
-# With thanks to Franz Smith and Steve Wing for support
-#
-# PURPOSE: Like IDW interpolation, but distance is cost to get to any
-# other site.
-#
-# COPYRIGHT: (c) 2003-2011 Hamish Bowman
-# This program is free software under the GNU General Public
-# License (>=v2). Read the file COPYING that comes with GRASS
-# for details.
-#
-#############################################################################
-#
-# Description:
-# Non-euclidean, non-polluting IDW from areas separated by null cells
-# e.g.
-# - two parallel lakes, connected at one end (water chemistry in arms of a fjord)
-# - Population abundance constrained by topography (Kiwis crossing a land-bridge)
-#
-# w(d)= 1/d^p where p is user definable, usually 2.
-#
-# Notes:
-# A cost surface containing molasses barrier data may be used as well.
-# Input data points need not have direct line of sight to each other.
-# Try and keep the number of input sites to under a few dozen, as the
-# process is very computationally expensive. You might consider creating
-# a few hundred MB RAM-disk for a temporary mapset to avoid thrashing your
-# hard drive (r.cost is heavy on disk IO).
-
-
-#%Module
-#% description: IDW interpolation, but distance is cost to get to any other site
-#%End
-#%option
-#% key: input
-#% type: string
-#% gisprompt: old,vector,vector
-#% description: Name of existing vector points map containing seed data
-#% required : yes
-#%end
-#%option
-#% key: column
-#% type: string
-#% description: Column name in points map that contains data values
-#% required : yes
-#%end
-#%option
-#% key: output
-#% type: string
-#% gisprompt: new,cell,raster
-#% description: Name for output raster map
-#% required : yes
-#%end
-#%option
-#% key: cost_map
-#% type: string
-#% gisprompt: old,cell,raster
-#% description: Name of existing raster map containing cost information
-#% required : yes
-#%end
-#%option
-#% key: friction
-#% type: double
-#% description: Friction of distance, (the 'n' in 1/d^n)
-#% answer: 2
-#% options: 1-6
-#% required : no
-#%end
-#%option
-#% key: layer
-#% type: integer
-#% answer: 1
-#% description: Layer number of data in points map
-#% required: no
-#%end
-#%option
-#% key: where
-#% type: string
-#% label: WHERE conditions of SQL query statement without 'where' keyword
-#% description: Example: income < 1000 and inhab >= 10000
-#% required : no
-#%end
-
-##%option
-##% key: max_cost
-##% type: integer
-##% description: Optional maximum cumulative cost before setting weight to zero
-##% required : no
-##%end
-
-#%option
-#% key: post_mask
-#% type: string
-#% gisprompt: old,cell,raster
-#% description: Name of existing raster map to be used as post-processing MASK
-#% required : no
-#%end
-#%flag
-#% key: r
-#% description: Use (d^n)*log(d) instead of 1/(d^n) for radial basis function
-#%end
-#%flag
-#% key: v
-#% description: Verbose mode (to be removed soon; use --verbose instead)
-#%end
-
-
-
-if [ -z "$GISBASE" ] ; then
- echo "You have to be in GRASS to use this." 1>&2
- exit 1
-fi
-
-if [ "$1" != "@ARGS_PARSED@" ] ; then
- exec "$GISBASE/bin/g.parser" "$0" "$@"
-fi
-
-if [ -z "$GIS_OPT_POST_MASK" ] ; then
- eval `g.findfile element=cell file=MASK`
- if [ -n "$name" ] ; then
- echo "A MASK already exists; remove it before using the post_mask option." 1>&2
- exit 1
- fi
-fi
-
-echo
-echo "v.surf.icw -- Inverse Cost Weighted Interpolation"
-echo "Processing $GIS_OPT_INPUT -> $GIS_OPT_OUTPUT, column=$GIS_OPT_COLUMN, Cf=$GIS_OPT_FRICTION"
-if [ "$GIS_FLAG_R" -eq 1 ] ; then
- echo "Using (d^n)*log(d) radial basis function."
-fi
-echo "------------------------------------------------------------------------"
-
-FRICTION="$GIS_OPT_FRICTION"
-#if [ "$FRICTION" -lt 1 ] || [ "$FRICTION" -gt 6 ]; then
-# echo "Friction of Distance out of range [1-6]" 1>&2
-# exit 1
-#fi
-
-# adjust so that tiny numbers don't hog all the FP precision space
-#DIVISOR=""
-#if [ "$FRICTION" -eq 4 ] ; then DIVISOR=" / 10.0" ; fi
-#if [ "$FRICTION" -eq 5 ] ; then DIVISOR=" / 100.0" ; fi
-#if [ "$FRICTION" -eq 6 ] ; then DIVISOR=" / 500.0" ; fi
-# fp version:
-DIVISOR=`echo "$FRICTION" | awk '{if ( $1 > 4 ) {print " / " 0.01 * $1^6 }}'`
-
-
-POINTS_FILE="$GIS_OPT_INPUT"
-
-
-if [ $GIS_FLAG_V -eq 1 ] || \
- ( [ -n "$GRASS_VERBOSE" ] && [ "$GRASS_VERBOSE" -gt 1 ] ) ; then
- VERBOSE="-v"
-else
- VERBOSE=""
-fi
-
-### Check that we have the column and it is the correct type
-COL_NAME=`v.info -c "$POINTS_FILE" layer="$GIS_OPT_LAYER" 2> /dev/null | grep "|${GIS_OPT_COLUMN}$" | cut -f2 -d'|'`
-if [ -z "$COL_NAME" ] ; then
- echo "Data column [$GIS_OPT_COLUMN] not found in [$POINTS_FILE]." 1>&2
- exit
-fi
-COL_TYPE=`v.info -c "$POINTS_FILE" layer="$GIS_OPT_LAYER" 2> /dev/null | grep "|${GIS_OPT_COLUMN}$" | cut -f1 -d'|'`
-if [ "$COL_TYPE" != "DOUBLE PRECISION" ] && [ "$COL_TYPE" != "INTEGER" ] ; then
- echo "Data column [$COL_NAME] must be numeric." 1>&2
- exit
-fi
-
-### cleanse cost area mask to a flat =1 for my porpoises
-GRASS_OVERWRITE=1 \
- r.mapcalc "tmp_icw_area_$$ = if(\"$GIS_OPT_COST_MAP\", 1, null())"
-###
-
-########################################################################
-## Commence crunching ..
-TMP_POINTS=`g.tempfile pid=$$`
-if [ $? -ne 0 ] || [ -z "$TMP_POINTS" ] ; then
- echo "ERROR: unable to create temporary files" 1>&2
- exit 1
-fi
-
-#crop out only points in region
-# GRASS 6.3+: use v.out.ascii -r ?
-v.in.region output=tmp_icw_region_$$ > /dev/null
-v.select ainput="$POINTS_FILE" alayer="$GIS_OPT_LAYER" atype=point \
- binput=tmp_icw_region_$$ btype=area output=tmp_icw_points_sel_$$
-
-if [ -n "$GIS_OPT_WHERE" ] ; then
- v.extract in=tmp_icw_points_sel_$$ layer="$GIS_OPT_LAYER" \
- where="$GIS_OPT_WHERE" out=tmp_icw_points_$$
-else
- g.rename vect=tmp_icw_points_sel_$$,tmp_icw_points_$$
-fi
-
-v.out.ascii tmp_icw_points_$$ > "$TMP_POINTS"
-#db.select tmp_icw_points_$$ > "$TMP_TABLE"
-
-# find n
-#N="`v.info tmp_icw_points_$$ | grep "Number of points" | cut -f2 -d: | awk '{print $1}'`"
-N="`cat "$TMP_POINTS" | wc -l`"
-
-#figure out which column cats are in (3 or 4 depending on if there is a z coordinate)
-CATCOL="`cat "$TMP_POINTS" | head -1 | tr '|' '\n' | wc -l`"
-
-
-if [ "$N" -gt 50 ] ; then
- echo "WARNING: Computation is expensive! Please consider fewer points." 1>&2
- echo "Press enter to continue.." 1>&2
- read -t 30
-fi
-
-# gen cost maps for each site in range
-echo "Generating cost maps.." 1>&2
-NUM=1
-for POS in `cat "$TMP_POINTS"` ; do
-
-# echo "POS=[$POS]"
- EASTING=`echo "$POS" | cut -f1 -d"|"`
- NORTHING=`echo "$POS" | cut -f2 -d"|"`
- CAT=`echo "$POS" | cut -f"$CATCOL" -d"|"`
- DATA_VALUE="`v.db.select -c tmp_icw_points_$$ column="$COL_NAME" where="cat=${CAT}"`"
-
- echo "Site $NUM of $N e=$EASTING n=$NORTHING data=$DATA_VALUE"
-
- if [ -z "$DATA_VALUE" ] ; then
- echo " Skipping, no data here." 1>&2
- N=`expr $N - 1`
- continue
- fi
-
- if [ -z "`r.what input=tmp_icw_area_$$ east_north=$EASTING,$NORTHING | grep -v "*"`" ] ; then
- echo " Skipping, site lays outside of cost_map." 1>&2
- N=`expr $N - 1`
- continue
- fi
-
-# echo "$EASTING $NORTHING $DATA_VALUE" | v.in.ascii output=tmp_idw_cost_site_$$ fs=space --o 2> /dev/null
- r.cost $VERBOSE -k in=tmp_icw_area_$$ output=cost_site.$NUM \
- coordinate=$EASTING,$NORTHING
- #max_cost="$GIS_OPT_MAX_COST" : commented out until r.null cleansing/continue code is sorted out
- #start_points=tmp_idw_cost_site_$$
-
- # so the divisor exists and the weighting is huge at the exact sample spots
- # more efficient to reclass to 1?
- GRASS_OVERWRITE=1 \
- r.mapcalc "cost_site.$NUM = if(cost_site.$NUM == 0, 0.1, cost_site.$NUM)"
- # r.to.vect then r.patch output
-# v.to.rast in=tmp_idw_cost_site_29978 out=tmp_idw_cost_val_$$ use=val val=10
-
- if [ "$GIS_FLAG_R" -eq 0 ] ; then
- # exp(3,2) is 3^2 etc. as is pow(3,2)
- # r.mapcalc "1by_cost_site_sqrd.$NUM = 1.0 / exp(cost_site.$NUM , $FRICTION)"
- EXPRESSION="1.0 / pow(cost_site.$NUM $DIVISOR, $FRICTION )"
- else
- # use log10() or ln() ?
- EXPRESSION="1.0 / ( pow(cost_site.$NUM, $FRICTION) * log (cost_site.$NUM) )"
- fi
-# echo "r.mapcalc expression is: [$EXPRESSION]"
- GRASS_OVERWRITE=1 \
- r.mapcalc "1by_cost_site_sqrd.$NUM = $EXPRESSION"
-
-# r.patch in=1by_cost_site_sqrd.${NUM},tmp_idw_cost_val_$$ out=1by_cost_site_sqrd.${NUM} --o
-# g.remove rast=cost_site.$NUM
-
- NUM="`expr $NUM + 1`"
-done
-
-# 3) find sum(cost^2)
-echo 1>&2
-echo "Finding sum of squares.." 1>&2
-INPUT_MAPS=1by_cost_site_sqrd.1
-NUM=2
-while [ $NUM -le $N ] ; do
- INPUT_MAPS="$INPUT_MAPS,1by_cost_site_sqrd.$NUM"
- NUM="`expr $NUM + 1`"
-done
-
-if [ -n "$GIS_OPT_POST_MASK" ] ; then
- echo "Setting post_mask [$GIS_OPT_POST_MASK]" 1>&2
- GRASS_OVERWRITE=1 \
- r.mapcalc "MASK = \"$GIS_OPT_POST_MASK\""
-fi
-
-echo "Summation of cost weights.." 1>&2
-r.series method=sum in="$INPUT_MAPS" out=sum_of_1by_cost_sqs
-
-if [ -n "$GIS_OPT_POST_MASK" ] ; then
- echo "Removing post_mask [$GIS_OPT_POST_MASK]" 1>&2
- g.remove MASK | grep REMOVE
-fi
-
-# 4) ( 1/di^2 / sum(1/d^2) ) * ai
-echo 1>&2
-echo "Creating partial weights.." 1>&2
-NUM=1
-for POS in `cat "$TMP_POINTS"` ; do
-
- EASTING=`echo "$POS" | cut -f1 -d"|"`
- NORTHING=`echo "$POS" | cut -f2 -d"|"`
- CAT=`echo "$POS" | cut -f"$CATCOL" -d"|"`
- DATA_VALUE="`v.db.select -c tmp_icw_points_$$ column="$COL_NAME" where="cat=${CAT}"`"
-
- echo "Site $NUM of $N data value = $DATA_VALUE" 1>&2
-
- if [ -z "$DATA_VALUE" ] ; then
- echo " Skipping, no data here." 1>&2
- continue
- fi
-
- if [ -z "`r.what input=tmp_icw_area_$$ east_north=$EASTING,$NORTHING | grep -v "*"`" ] ; then
- echo " Skipping, site lays outside of cost_map." 1>&2
- continue
- fi
-
- GRASS_OVERWRITE=1 \
- r.mapcalc "partial_icw.$NUM = \
- ( $DATA_VALUE * 1by_cost_site_sqrd.$NUM ) / sum_of_1by_cost_sqs"
-#"( $DATA_VALUE / $N ) * (1.0 - ( cost_sq_site.$NUM / sum_of_cost_sqs ))"
-#"( cost_sq_site.$NUM / sum_of_cost_sqs ) * ( $DATA_VALUE / $N )"
-
-# g.remove rast=1by_cost_site_sqrd.$NUM
- NUM="`expr $NUM + 1`"
-
- if [ $NUM -gt $N ] ; then
- break
- fi
-
-done
-
-INPUT_MAPS=partial_icw.1
-NUM=2
-while [ $NUM -le $N ] ; do
- INPUT_MAPS=$INPUT_MAPS,partial_icw.$NUM
- NUM="`expr $NUM + 1`"
-done
-
-echo 1>&2
-echo "Calculating final values.." 1>&2
-r.series method=sum in=$INPUT_MAPS out="$GIS_OPT_OUTPUT"
-
-#TODO: r.patch in v.to.rast of values at exact seed site locations. currently set to null
-
-r.colors "$GIS_OPT_OUTPUT" rule=bcyr
-
-r.support "$GIS_OPT_OUTPUT" title="Inverse cost-weighted interpolation" history=""
-r.support "$GIS_OPT_OUTPUT" history="v.surf.icw interpolation:"
-r.support "$GIS_OPT_OUTPUT" history=" input map=$GIS_OPT_INPUT attribute column=$GIS_OPT_COLUMN"
-r.support "$GIS_OPT_OUTPUT" history=" cost map=$GIS_OPT_COST_MAP coefficient of friction=$FRICTION"
-if [ $GIS_FLAG_R -eq 1 ] ; then
- r.support "$GIS_OPT_OUTPUT" history=" (d^n)*log(d) as radial basis function"
-fi
-if [ -n "$GIS_OPT_POST_MASK" ] ; then
- r.support "$GIS_OPT_OUTPUT" history=" post-processing mask=$GIS_OPT_POST_MASK"
-fi
-if [ -n "$GIS_OPT_WHERE" ] ; then
- r.support "$GIS_OPT_OUTPUT" history=" SQL query= WHERE \"$GIS_OPT_WHERE\""
-fi
-# save layer #? to metadata? command line hist?
-
-
-# 5) rm cost and cost_sq maps, tmp_icw_points, etc
-cleanup()
-{
- echo "Cleanup.." 1>&2
- g.mremove -f rast=cost_site.*
- g.mremove -f rast=1by_cost_site_sqrd.*
- g.mremove -f rast=partial_icw.*
- #g.mremove -f rast=tmp_icw_*_$$
- #g.mremove -f vect=tmp_icw_*_$$
- g.remove rast=sum_of_1by_cost_sqs
- #g.remove vect=tmp_idw_cost_site_$$
- g.remove rast=tmp_icw_area_$$
- g.remove vect=tmp_icw_region_$$,tmp_icw_points_$$
- # check if it exists
- eval `g.findfile element=vector file="tmp_icw_points_sel_$$"`
- if [ -e "$file" ] ; then
- g.remove vect=tmp_icw_points_sel_$$
- fi
- rm -f "$TMP_POINTS"
-}
-
-# TODO: trap ^C
-# what to do in case of user break:
-#exitprocedure()
-#{
-# g.message -e 'User break!'
-# cleanup
-# exit 1
-#}
-# shell check for user break (signal list: trap -l)
-#trap "exitprocedure" 2 3 15
-
-cleanup
-
-# 6) done!
-echo "Done! Results written to <${GIS_OPT_OUTPUT}>." 1>&2
-exit 0
Copied: grass-addons/grass7/vector/v.surf.icw/v.surf.icw.html (from rev 49599, grass-addons/grass6/vector/v.surf.icw/description.html)
===================================================================
--- grass-addons/grass7/vector/v.surf.icw/v.surf.icw.html (rev 0)
+++ grass-addons/grass7/vector/v.surf.icw/v.surf.icw.html 2012-05-22 10:40:41 UTC (rev 51665)
@@ -0,0 +1,113 @@
+<h2>DESCRIPTION</h2>
+
+Inverse cost weighting is like inverse distance weighted (IDW) interpolation,
+but uses <i>cost</i> instead of shortest Euclidean distance. In this way
+solid barriers and molasses zones may be correctly taken into account.
+<p>
+Input data points do not need to have direct line of sight to each other.
+This is interpolation "as the fish swims", not "as the crow flies", and can
+see around headlands or across land-bridges without polluting over barriers
+which the natural value (or flightless bird) can not cross.
+<p>
+It was initially written to interpolate water chemistry in two parallel
+arms of a fjord, but may just as well be used for population abundance
+constrained by topography, or in studies of archeologic technology transfer.
+
+
+<h2>NOTES</h2>
+
+In the simplest case, the cost map will just be a mask raster with values
+of 1 in areas to interpolate, and NULL in impenetrable areas.
+Fancier cost maps can be used, for example, to make it more expensive for
+a measured pollutant to diffuse upstream in an estuary, or to make it more
+expensive for a stone tool technology to cross waterways.
+<p>
+Since generating cost maps can take a long time, it is recommended to keep
+the raster region relatively small and limit the number of starting points
+to less than 100.
+<p>
+Higher values of <b>friction</b> will help limit unconstrained boundary
+effects at the edges of your coverage, but will incur more of a stepped
+transition between sites.
+<p>
+The <b>post_mask</b>, if given, is applied after the interpolation is
+complete. A common use for that might be to only present data within
+a certain distance (thus confidence) of an actual sampling station.
+In that case the <em>r.cost</em> module can be used to create the mask.
+<p>
+This module writes lots of temporary files and so can be rather disk I/O
+intensive. If you are running it many times in a big loop you may want to
+set up a RAM-disk to put the mapset in (on UNIX symlinking back into the
+location is ok), or adjust the disk-cache flushing timer to be slightly
+longer than one iteration of the script.
+<br>To do this on Linux you can tune the
+ <tt>/proc/sys/vm/dirty_expire_centisecs</tt> kernel control. The default
+is to keep files in memory a maximum of 30 seconds before writing them to
+disk, although if you are short on free RAM the kernel may write to disk
+long before that timeout is reached.
+<!-- Linux I/O tuning hints:
+ http://www.westnet.com/~gsmith/content/linux-pdflush.htm
+ http://tldp.org/LDP/sag/html/buffer-cache.html
+ -->
+
+
+<h2>EXAMPLE</h2>
+
+In this example we'll generate some fake island barriers from the
+standard North Carolina GRASS dataset, then interpolate a continuous
+variable between given point stations. We'll use rainfall, but for the
+purposes of the exercise pretend it is a nutrient concentration in our
+fake wayerway. Point data stations outside of the current region or
+in NULL areas of the <i>cost_map</i> will be ignored.
+
+<div class="code"><pre>
+# set up a fake sea with islands:
+g.region n=276000 s=144500 w=122000 e=338500 res=500
+r.mapcalc "pseudo_elev = elev_state_500m - 1100"
+r.colors pseudo_elev color=etopo2
+r.mapcalc "navigable_mask = if(pseudo_elev < 0, 1, null())"
+
+# pick a data column from the points vector:
+v.info -c precip_30ynormals
+
+# run the interpolation:
+v.surf.icw input=precip_30ynormals column=annual output=annual_interp.3 \
+ cost_map=navigable_mask friction=3 --quiet
+
+# equalize colors to show maximum detail:
+r.colors -e annual_interp.3 color=bcyr
+
+# display results in an Xmonitor:
+d.erase black
+d.rast -o annual_interp.3
+d.vect precip_30ynormals fcolor=red icon=basic/circle
+d.legend annual_interp.3 color=white at=48.4,94.8,3.4,6.0
+</pre></div>
+
+
+<h2>REFERENCES</h2>
+
+ * Hamish's MfE reports / published confernce proceedings?
+ * Ben D's journal article (Deutsch)
+
+
+<h2>SEE ALSO</h2>
+
+<em>
+<a href="v.surf.idw.html">v.surf.idw</a><br>
+<a href="v.surf.rst.html">v.surf.rst</a><br>
+<a href="v.surf.bspline.html">v.surf.bspline</a><br>
+<a href="r.cost.html">r.cost</a><br>
+<a href="r.surf.idw.html">r.surf.idw</a><br>
+<a href="r.surf.idw2.html">r.surf.idw2</a><br>
+</em>
+
+
+<h2>AUTHOR</h2>
+
+Hamish Bowman<br>
+<i>Department of Marine Science,<br>
+Dunedin, New Zealand</i>
+
+<p>
+<i>Last changed: $Date$</i>
Copied: grass-addons/grass7/vector/v.surf.icw/v.surf.icw.py (from rev 49599, grass-addons/grass6/vector/v.surf.icw/v.surf.icw)
===================================================================
--- grass-addons/grass7/vector/v.surf.icw/v.surf.icw.py (rev 0)
+++ grass-addons/grass7/vector/v.surf.icw/v.surf.icw.py 2012-05-22 10:40:41 UTC (rev 51665)
@@ -0,0 +1,475 @@
+#!/usr/bin/env python
+#############################################################################
+#
+# MODULE: v.surf.icw
+# version $Id$
+#
+# AUTHOR: M. Hamish Bowman, Dunedin, New Zealand
+# Originally written aboard the NZ DoC ship M/V Renown,
+# Bligh Sound, Fiordland National Park, November 2003
+# With thanks to Franz Smith and Steve Wing for support
+# Ported to Python for GRASS 7 December 2011
+#
+# PURPOSE: Like IDW interpolation, but distance is cost to get to any
+# other site.
+#
+# COPYRIGHT: (c) 2003-2011 Hamish Bowman
+# This program is free software under the GNU General Public
+# License (>=v2). Read the file COPYING that comes with GRASS
+# for details.
+#
+#############################################################################
+#
+# Description:
+# Non-euclidean, non-polluting IDW from areas separated by null cells
+# e.g.
+# - two parallel lakes, connected at one end (water chemistry in arms of a fjord)
+# - Population abundance constrained by topography (Kiwis crossing a land-bridge)
+#
+# w(d)= 1/d^p where p is user definable, usually 2.
+#
+# Notes:
+# A cost surface containing molasses barrier data may be used as well.
+# Input data points need not have direct line of sight to each other.
+# Try and keep the number of input sites to under a few dozen, as the
+# process is very computationally expensive. You might consider creating
+# a few hundred MB RAM-disk for a temporary mapset to avoid thrashing your
+# hard drive (r.cost is heavy on disk IO).
+
+
+#%Module
+#% description: IDW interpolation, but distance is cost to get to any other site.
+#%End
+#%option
+#% key: input
+#% type: string
+#% gisprompt: old,vector,vector
+#% description: Name of existing vector points map containing seed data
+#% required : yes
+#%end
+#%option
+#% key: column
+#% type: string
+#% description: Column name in points map that contains data values
+#% required : yes
+#%end
+#%option
+#% key: output
+#% type: string
+#% gisprompt: new,cell,raster
+#% description: Name for output raster map
+#% required : yes
+#%end
+#%option
+#% key: cost_map
+#% type: string
+#% gisprompt: old,cell,raster
+#% description: Name of existing raster map containing cost information
+#% required : yes
+#%end
+#%option
+#% key: friction
+#% type: double
+#% description: Friction of distance, (the 'n' in 1/d^n)
+#% answer: 2
+#% options: 1-6
+#% required : no
+#%end
+#%option
+#% key: layer
+#% type: integer
+#% answer: 1
+#% description: Layer number of data in points map
+#% required: no
+#%end
+#%option
+#% key: where
+#% type: string
+#% label: WHERE conditions of SQL query statement without 'where' keyword
+#% description: Example: income < 1000 and inhab >= 10000
+#% required : no
+#%end
+
+##%option
+##% key: max_cost
+##% type: integer
+##% description: Optional maximum cumulative cost before setting weight to zero
+##% required : no
+##%end
+
+#%option
+#% key: post_mask
+#% type: string
+#% gisprompt: old,cell,raster
+#% description: Name of existing raster map to be used as post-processing MASK
+#% required : no
+#%end
+#%flag
+#% key: r
+#% description: Use (d^n)*log(d) instead of 1/(d^n) for radial basis function
+#%end
+#%flag
+#% key: v
+#% description: Verbose mode (to be removed soon; use --verbose instead)
+#%end
+
+
+import sys
+import os
+from grass.script import core as grass
+
+def main():
+
+ pts_input = options['input']
+ output = options['output']
+ cost_map = options['cost_map']
+ post_mask = options['post_mask']
+ column = options['column']
+ friction = options['friction']
+ layer = options['layer']
+ where = options['where']
+
+ pid = str(os.getpid())
+ tmp_base = 'tmp_icw_' + pid + '_'
+
+ # does map exist?
+ if not grass.find_file(pts_input, element = 'vector')['file']:
+ grass.fatal(_("Vector map <%s> not found") % pts_input)
+
+ if grass.findfile(post_mask):
+ grass.fatal(_("A MASK already exists; remove it before using the post_mask option."))
+
+ grass.verbose(_("v.surf.icw -- Inverse Cost Weighted Interpolation"))
+ grass.verbose(_("Processing %s -> %s, column=%s, Cf=%g") % (pts_input, output, column, friction))
+
+ if flags['r']:
+ grass.verbose(_("Using (d^n)*log(d) radial basis function."))
+
+ grass.verbose("------------------------------------------------------------------------")
+
+ # adjust so that tiny numbers don't hog all the FP precision space
+ #if friction = 4:
+ # divisor = 10.0
+ #if friction = 5:
+ # divisor = 100.0
+ #if friction = 6:
+ # divisor = 500.0
+ # fp version:
+ if friction > 4:
+ divisor = 0.01 * pow(friction, 6)
+ else:
+ divisor = 1
+
+ # Check that we have the column and it is the correct type
+ try:
+ coltype = grass.vector_columns(pts_input, layer)[column]
+ except KeyError:
+ grass.fatal(_("Data column <%s> not found in vector points map <%s>") % (column, pts_input))
+
+ if coltype['type'] not in ('INTEGER', 'DOUBLE PRECISION'):
+ grass.fatal(_("Data column must be numberic"))
+
+ # cleanse cost area mask to a flat =1 for my porpoises
+ grass.mapcalc("$result = if($cost_map, 1, null())",
+ result = tmp_base + 'area', cost_map = cost_map,
+ overwrite = True)
+
+
+ ########################################################################
+ ## Commence crunching ..
+ tmp_points = grass.tempfile()
+
+ #crop out only points in region
+ grass.run_command('v.in.region', output = tmp_base + 'region', quiet = True)
+ grass.run_command('v.select',
+ ainput = pts_file, alayer = layer, atype = 'point',
+ binput = tmp_base + 'region', btype = 'area',
+ output = tmp_base + 'points_sel')
+
+ if where:
+ grass.run_command('v.select', input = tmp_base + 'points_sel',
+ layer = layer, where = where,
+ output = tmp_base + 'points')
+ else:
+ grass.run_command('g.rename', rast = (tmp_base + 'points_sel',
+ tmp_base + 'points'), quiet = True)
+
+ # use v.out.ascii -r instead?
+ if where:
+ grass.run_command('v.out.ascii', input = pts_file, output = tmp_points,
+ where = where, flags = 'r')
+ else:
+ grass.run_command('v.out.ascii', input = pts_file, output = tmp_points,
+ flags = 'r')
+
+ #db.select tmp_icw_points_$$ > "$TMP_TABLE"
+
+ # count number of starting points (n)
+ n = grass.vector_info_topo(tmp_base + 'points')['points']
+ #n = sum(1 for line in open(tmp_points))
+
+ #figure out which column cats are in (3 or 4 depending on if there is a z coordinate)
+ cat_col = grass.vector_info_topo(tmp_base + 'points')['map3d'] + 3
+
+ if n > 200:
+ grass.warning(_("Computation is expensive! Please consider fewer points or get ready to wait a while ..."))
+ sleep 5
+
+ # generate cost maps for each site in range
+ grass.verbose(_("Generating cost maps.."))
+
+ num = 1
+
+####
+#see osm python script
+
+infile = tmp_points
+if not os.path.exists(infile):
+ grass.fatal("Unable to read input data")
+inf_id = file(infile)
+print("input file=[%s]" % infile)
+lines = inf_id.readlines()
+for i in range(len(lines)):
+ lines[i] = lines[i].rstrip('\n')
+for line in lines:
+ bits = line.split('|')
+bits[0], bits[1], bits[2], float(bits[3]), ....
+inf_id.close()
+
+or
+
+ while True:
+ line = inf_old.readline()
+ #.strip()
+ if not line:
+ break
+
+ if 'node id=' not in line:
+ continue
+
+newid_oldid_lonlat = [[] for i in range(4)]
+
+for i in range(len(old_id_lonlat[0])):
+ if abs(new_lon - old_id_lonlat[1][i]) < thresh_lon and \
+ abs(new_lat - old_id_lonlat[2][i]) < thresh_lat:
+ newid_oldid_lonlat[0].append(new_id)
+ newid_oldid_lonlat[1].append(old_id_lonlat[0][i])
+ newid_oldid_lonlat[2].append(old_id_lonlat[1][i])
+ newid_oldid_lonlat[3].append(old_id_lonlat[2][i])
+
+bits = string.split(line,'"')
+
+####
+
+for POS in `cat "$TMP_POINTS"` ; do
+
+# echo "POS=[$POS]"
+ EASTING=`echo "$POS" | cut -f1 -d"|"`
+ NORTHING=`echo "$POS" | cut -f2 -d"|"`
+ CAT=`echo "$POS" | cut -f"$CATCOL" -d"|"`
+ DATA_VALUE="`v.db.select -c tmp_icw_points_$$ column="$COL_NAME" where="cat=${CAT}"`"
+
+ data_value = grass.vector_select(map name)['values'][cat]
+
+ grass.message(_("Site %d of %d e=%.4f n=%.4f data=%.3g"),
+ num, n, easting, northing, data_value)
+
+ if not data_value:
+ grass.messsage(_(" Skipping, no data here."))
+ n = n - 1
+ continue
+
+ if [ -z "`r.what input=tmp_icw_area_$$ east_north=$EASTING,$NORTHING | grep -v "*"`" ] ; then
+ grass.messsage(_(" Skipping, point lays outside of cost_map."))
+ n = n - 1
+ continue
+ fi
+
+ # echo "$EASTING $NORTHING $DATA_VALUE" | v.in.ascii output=tmp_idw_cost_site_$$ fs=space --o 2> /dev/null
+
+ cost_site_name = tmp_base + 'cost_site.' + str(num)
+ grass.run_command('r.cost', flags = 'k', input = tmp_base + 'area',
+ output = cost_site_name, coordinate = easting + ',' + northing)
+ #max_cost="$GIS_OPT_MAX_COST" : commented out until r.null cleansing/continue code is sorted out
+ #start_points=tmp_idw_cost_site_$$
+
+ # we do this so the divisor exists and the weighting is huge at the exact sample spots
+ # more efficient to reclass to 1?
+ grass.mapcalc("$cost_n_cleansed = if($cost_n == 0, 0.1, $cost_n",
+ cost_n_cleansed = cost_site_name + '.cleansed',
+ cost_n = cost_site_name, overwrite = True)
+ grass.run_command('g.remove', rast = cost_site_name, quiet = True)
+ grass.run_command('g.rename', rast = cost_site_name + '.cleansed'
+ + ',' + cost_site_name, quiet = True)
+
+
+ # r.to.vect then r.patch output
+ # v.to.rast in=tmp_idw_cost_site_29978 out=tmp_idw_cost_val_$$ use=val val=10
+
+ if not flags['r']:
+ # exp(3,2) is 3^2 etc. as is pow(3,2)
+ # r.mapcalc "1by_cost_site_sqrd.$NUM = 1.0 / exp(cost_site.$NUM , $FRICTION)"
+# EXPRESSION="1.0 / pow(cost_site.$NUM $DIVISOR, $FRICTION )"
+ expr = '1.0 / pow($cost_n / ' + str(divisor) + ', $friction)'
+ else
+ # use log10() or ln() ?
+# EXPRESSION="1.0 / ( pow(cost_site.$NUM, $FRICTION) * log (cost_site.$NUM) )"
+ expr = '1.0 / ( pow($cost_n, $friction) * log($cost_n) )"'
+
+ grass.debug("r.mapcalc expression is: [%s]" % expr)
+ grass.mapcalc("$result.$num = " + expr,
+ result = tmp_base + '1by_cost_site_sqrd',
+ num = str(num), cost_n = cost_site_name,
+ friction = friction, overwrite = True)
+
+ # r.patch in=1by_cost_site_sqrd.${NUM},tmp_idw_cost_val_$$ out=1by_cost_site_sqrd.${NUM} --o
+ # g.remove rast=cost_site.$NUM
+
+ num = num + 1
+ done
+
+
+ # Step 3) find sum(cost^2)
+ grass.message("")
+ grass.message(_("Finding sum of squares.."))
+
+INPUT_MAPS=1by_cost_site_sqrd.1
+NUM=2
+while [ $NUM -le $N ] ; do
+ INPUT_MAPS="$INPUT_MAPS,1by_cost_site_sqrd.$NUM"
+ NUM="`expr $NUM + 1`"
+done
+
+ if post_mask:
+ grass.message(_("Setting post_mask <%s>"), post_mask)
+ grass.mapcalc("MASK = $maskmap", maskmap = post_mask, overwrite = True)
+
+ grass.message(_("Summation of cost weights.."))
+ grass.run_command('r.series', method = 'sum', input = input_maps,
+ output = tmp_base + 'sum_of_1by_cost_sqs')
+
+ if post_mask:
+ grass.message(_("Removing post_mask <%s>"), post_mask)
+ grass.run_command('g.remove', rast = 'MASK', quiet = True)
+
+ # Step 4) ( 1/di^2 / sum(1/d^2) ) * ai
+ grass.message("")
+ grass.message(_("Creating partial weights.."))
+
+ num = 1
+for POS in `cat "$TMP_POINTS"` ; do
+
+ EASTING=`echo "$POS" | cut -f1 -d"|"`
+ NORTHING=`echo "$POS" | cut -f2 -d"|"`
+ CAT=`echo "$POS" | cut -f"$CATCOL" -d"|"`
+ DATA_VALUE="`v.db.select -c tmp_icw_points_$$ column="$COL_NAME" where="cat=${CAT}"`"
+
+ grass.message(_("Site %d of %d, data value = %.3f") % (num, n, data_value))
+
+ if not data_value:
+ grass.message(_(" Skipping, no data here."))
+ continue
+
+ if [ -z "`r.what input=tmp_icw_area_$$ east_north=$EASTING,$NORTHING | grep -v "*"`" ] ; then
+ grass.message(_(" Skipping, site lays outside of cost_map."))
+ continue
+ fi
+
+ partial_n = tmp_base + 'partial.' + str(num)
+ grass.mapcalc("$partial_n = ($data * $1by_cost_sq) / $sum_of_1by_cost_sqs",
+ partial_n = partial_n, data = data_value, overwrite = True)
+ r.mapcalc "$result = ( $data_value * $one_by_cost_site_sqrd ) / \
+ $sum_of_1by_cost_sqs",
+ data_value = data_value,
+ result = tmp_base + 'partial_icw.' + str(num),
+ one_by_cost_site_sqrd = tmp_base + '1by_cost_site_sqrd' + str(num),
+ sum_of_1by_cost_sqs = tmp_base + 'sum_of_1by_cost_sqs', overwrite = True)
+ #"( $DATA_VALUE / $N ) * (1.0 - ( cost_sq_site.$NUM / sum_of_cost_sqs ))"
+ #"( cost_sq_site.$NUM / sum_of_cost_sqs ) * ( $DATA_VALUE / $N )"
+
+ # g.remove rast=1by_cost_site_sqrd.$NUM
+
+ num = num + 1
+ if num > n:
+ break
+ # done
+
+INPUT_MAPS=partial_icw.1
+NUM=2
+while [ $NUM -le $N ] ; do
+ INPUT_MAPS=$INPUT_MAPS,partial_icw.$NUM
+ NUM="`expr $NUM + 1`"
+done
+
+ grass.message("")
+ grass.message(_("Calculating final values.."))
+
+ grass.run_command('r.series', method = 'sum', input = input_maps,
+ output = output)
+
+ #TODO: r.patch in v.to.rast of values at exact seed site locations. currently set to null
+
+ grass.run_command('r.colors', map = output, color = 'bcyr')
+ grass.run_command('r.support', map = output, history = '',
+ title = 'Inverse cost-weighted interpolation')
+ grass.run_command('r.support', map = output,
+ history = 'v.surf.icw interpolation:')
+ grass.run_command('r.support', map = output,
+ history = ' input map=' + input + ' attribute column=' + column)
+ grass.run_command('r.support', map = output,
+ history = ' cost map=' + cost_map + ' coefficient of friction=' + friction)
+ if flags['r']:
+ grass.run_command('r.support', map = output,
+ history = ' (d^n)*log(d) as radial basis function')
+ if post_mask:
+ grass.run_command('r.support', map = output,
+ history = ' post-processing mask=' + post_mask)
+ if where:
+ grass.run_command('r.support', map = output,
+ history = ' SQL query= WHERE ' + where)
+
+ # save layer #? to metadata? command line hist?
+
+
+ # Step 5) rm cost and cost_sq maps, tmp_icw_points, etc
+cleanup:
+ grass.verbose(_("Cleanup.."))
+ grass.run_command('g.mremove', flags = 'f', rast = tmp_base + '*', quiet = True)
+
+ g.mremove -f rast=cost_site.*
+ g.mremove -f rast=1by_cost_site_sqrd.*
+ g.mremove -f rast=partial_icw.*
+ #g.mremove -f rast=tmp_icw_*_$$
+ #g.mremove -f vect=tmp_icw_*_$$
+ g.remove rast=sum_of_1by_cost_sqs
+ #g.remove vect=tmp_idw_cost_site_$$
+ g.remove rast=tmp_icw_area_$$
+ g.remove vect=tmp_icw_region_$$,tmp_icw_points_$$
+ # check if it exists
+ eval `g.findfile element=vector file="tmp_icw_points_sel_$$"`
+ if [ -e "$file" ] ; then
+ g.remove vect=tmp_icw_points_sel_$$
+ fi
+ rm -f "$TMP_POINTS"
+}
+
+# TODO: trap ^C
+# what to do in case of user break:
+#exitprocedure()
+#{
+# g.message -e 'User break!'
+# cleanup
+# exit 1
+#}
+# shell check for user break (signal list: trap -l)
+#trap "exitprocedure" 2 3 15
+
+cleanup
+
+ # Step 6) done!
+ grass.message(_("Done! Results written to <%s>." % output)
+
+
+if __name__ == "__main__":
+ options, flags = grass.parser()
+ main()
Property changes on: grass-addons/grass7/vector/v.surf.icw/v.surf.icw.py
___________________________________________________________________
Added: svn:executable
+ *
Added: svn:mime-type
+ text/x-python
Added: svn:keywords
+ Id
Added: svn:eol-style
+ native
More information about the grass-commit
mailing list