[GRASS-SVN] r30348 - in grass-addons/vector: . v.surf.icw
svn_grass at osgeo.org
svn_grass at osgeo.org
Mon Feb 25 18:11:53 EST 2008
Author: hamish
Date: 2008-02-25 18:11:53 -0500 (Mon, 25 Feb 2008)
New Revision: 30348
Added:
grass-addons/vector/v.surf.icw/
grass-addons/vector/v.surf.icw/v.surf.icw
Log:
new: IDW interp correctly respecting hard full barriers and soft cost barriers
Added: grass-addons/vector/v.surf.icw/v.surf.icw
===================================================================
--- grass-addons/vector/v.surf.icw/v.surf.icw (rev 0)
+++ grass-addons/vector/v.surf.icw/v.surf.icw 2008-02-25 23:11:53 UTC (rev 30348)
@@ -0,0 +1,380 @@
+#!/bin/bash
+#
+###########################################################################
+#
+# v.surf.icw
+# version 30 May 2006
+#
+# Like IDW interpolation, but distance is cost to get to any other site
+# Based on s.surf.icw.sh for GRASS 5
+#
+# (c) 2003-2006 Hamish Bowman
+#
+#
+# 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
+#
+###########################################################################
+# Non-euclidian, non-polluting IDW from areas separated by null cells
+# e.g.
+# - two parallel lakes, connected at one end (water chemisty 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.
+#
+#
+
+#%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: output
+#% type: string
+#% gisprompt: new,cell,raster
+#% description: Name of new raster file to be used for output
+#% required : yes
+#%end
+
+#%option
+#% key: column
+#% type: string
+#% description: Column name in vector points map that contains data values
+#% required : yes
+#%end
+
+#%option
+#% key: cost_map
+#% type: string
+#% gisprompt: old,cell,raster
+#% description: Name of existing raster file containing cost information
+#% required : yes
+#%end
+
+#%option
+#% key: post_mask
+#% type: string
+#% gisprompt: old,cell,raster
+#% description: Name of existing raster file to be used as post-proccessing MASK
+#% required : no
+#%end
+
+#%option
+#% key: friction
+#% type: integer
+#% description: Friction of distance, (the 'n' in 1/d^n)
+#% answer: 2
+#% options: 1-6
+#% required : no
+#%end
+
+### needed?
+#%option
+#% key: null_string
+#% type: string
+#% description: Text string that represents "no data" in input data
+#% answer: nan
+#% required : no
+#%end
+# in the above, s.out.ascii changes a "%NaN" to "nan" so we use lower case
+
+#%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
+#%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"
+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
+
+unset 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
+
+#COLMN=$(( 4 + $GIS_OPT_index )) # skip over e,n,|,cat ### wrong for 3D+
+# maybe it's better to do for i=1:index sed delete line upto & incl. "%" then |cut -f1 ?
+
+NULL="$GIS_OPT_null_string"
+POINTS_FILE="$GIS_OPT_input"
+
+unset VERBOSE
+if [ $GIS_FLAG_v -eq 1 ] ; then VERBOSE="-v" ; fi
+
+#figure out which column number we are after
+#COLMN="`db.columns "$POINTS_FILE" | grep -n "^${GIS_OPT_column}$" | cut -f1 -d:`"
+#if [ -z "$COLMN" ] || [ "$COLMN" -le 0 ] ; then
+# echo "Data column [$GIS_OPT_column] not found in [$POINTS_FILE]." 1>&2
+# exit
+#fi
+
+### Check that we have the column and it is the correct type
+COL_NAME=`v.info -c "$POINTS_FILE" 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" 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
+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
+#s.out.ascii -d "$GIS_OPT_sites" > /tmp/idw_tmp.$$
+#s.out.ascii -d "$GIS_OPT_sites" > "$TMP_SITES"
+
+#crop out only points in region
+v.in.region output=tmp_icw_region_$$ > /dev/null
+v.select ainput="$POINTS_FILE" binput=tmp_icw_region_$$ \
+ atype=point btype=area output=tmp_icw_points_$$
+
+v.out.ascii tmp_icw_points_$$ > "$TMP_POINTS"
+#db.select tmp_icw_points_$$ > "$TMP_TABLE"
+
+# find n
+#N=`wc -l /tmp/idw_tmp.$$ | awk '{print $1}'`
+#N=`wc -l "$TMP_SITES" | awk '{print $1}'`
+#N=`s.out.ascii "$SITES_FILE" | wc -l`
+#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/idw_tmp.$$ | sed -e 's/ /|/g'` ; do
+#for POS in `cat "$TMP_SITES" | sed -e 's/ /|/g'` ; do
+#for POS in `s.out.ascii -d "$SITES_FILE" | sed -e 's/ /|/g'` ; do
+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=`echo "$POS" | cut -f"$COLMN" -d"|"`
+# DATA_VALUE="`db.select -c tmp_icw_points_$$ sql="select * from tmp_icw_points_$$ where cat = $CAT" | cut -f"$COLMN" -d'|'`"
+ DATA_VALUE="`db.select -c tmp_icw_points_$$ sql="select $COL_NAME from tmp_icw_points_$$ where cat=${CAT}"`"
+
+ echo site $NUM of $N e=$EASTING n=$NORTHING data=$DATA_VALUE
+
+ if [ -z "$DATA_VALUE" ] || [ "$DATA_VALUE" = "$NULL" ] ; 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_$$ out=cost_site.$NUM coordinate=$EASTING,$NORTHING
+ #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?
+ 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
+
+ # 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 )"
+ if [ "$GIS_FLAG_r" -eq 1 ] ; then
+ echo "Using (d^n)*log(d) radial basis function." 1>&2
+ EXPRESSION="pow(cost_site.$NUM, $FRICTION) * log (cost_site.$NUM)"
+ fi
+# echo "r.mapcalc expression is: [$EXPRESSION]"
+
+ 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
+
+#r.mapcalc resolution_mask="if(fjords_mask at all == 9 )"
+#r.mapcalc MASK="$FJORD"_mask # @all
+if [ -n "$GIS_OPT_post_mask" ] ; then
+ echo "Setting post_mask [$GIS_OPT_post_mask]" 1>&2
+ 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/idw_tmp.$$ | sed -e 's/ /|/g'` ; do
+#for POS in `cat "$TMP_SITES" | sed -e 's/ /|/g'` ; do
+#for POS in `v.out.ascii -d "$SITES_FILE" | sed -e 's/ /|/g'` ; do
+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="`db.select -c tmp_icw_points_$$ sql="select * from tmp_icw_points_$$ where cat = $CAT" | cut -f"$COLMN" -d'|'`"
+ DATA_VALUE="`db.select -c tmp_icw_points_$$ sql="select $COL_NAME from tmp_icw_points_$$ where cat=${CAT}"`"
+
+ echo "site $NUM of $N data value = $DATA_VALUE" 1>&2
+
+ if [ -z "$DATA_VALUE" ] || [ "$DATA_VALUE" = "$NULL" ] ; 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
+
+ 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"
+
+r.colors "$GIS_OPT_output" rule=bcyr
+if [ 0 -eq 1 ] ; then
+r.colors "$GIS_OPT_output" col=rules << EOF
+ 0 magenta
+ 10 red
+ 20 brown
+ 25 yellow
+ 30 green
+ 35 blue
+ 40 white
+EOF
+fi
+
+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
+
+
+# 5) rm cost and cost_sq maps, site_tmp, /tmp/idw_tmp.$$, etc
+#exit ######
+echo "Cleanup .." 1>&2
+g.mremove -f rast=cost_site.* | grep REMOVE
+g.mremove -f rast=1by_cost_site_sqrd.* | grep REMOVE
+g.mremove -f rast=partial_icw.* | grep REMOVE
+#g.mremove -f rast=tmp_icw_*_$$
+#g.mremove -f vect=tmp_icw_*_$$
+g.remove rast=sum_of_1by_cost_sqs | grep REMOVE
+#g.remove vect=tmp_idw_cost_site_$$
+g.remove rast=tmp_icw_area_$$ | grep REMOVE
+g.remove vect=tmp_icw_region_$$,tmp_icw_points_$$
+#rm /tmp/idw_tmp.$$
+#rm "$TMP_SITES"
+rm -f "$TMP_POINTS"
+
+# 6) done!
+echo "Done! Wrote results to <${GIS_OPT_output}>." 1>&2
+exit 0
Property changes on: grass-addons/vector/v.surf.icw/v.surf.icw
___________________________________________________________________
Name: svn:executable
+ *
More information about the grass-commit
mailing list