[GRASS-SVN] r52995 - in grass-addons/grass6/raster: . r.connectivity.distance
svn_grass at osgeo.org
svn_grass at osgeo.org
Thu Aug 30 11:27:59 PDT 2012
Author: sbl
Date: 2012-08-30 11:27:59 -0700 (Thu, 30 Aug 2012)
New Revision: 52995
Added:
grass-addons/grass6/raster/r.connectivity.distance/
grass-addons/grass6/raster/r.connectivity.distance/Makefile
grass-addons/grass6/raster/r.connectivity.distance/description.html
grass-addons/grass6/raster/r.connectivity.distance/r.connectivity.distance
Log:
A new addon for connectivity analysis based on graph theory
Property changes on: grass-addons/grass6/raster/r.connectivity.distance
___________________________________________________________________
Added: svn:ignore
+ :
Added: grass-addons/grass6/raster/r.connectivity.distance/Makefile
===================================================================
--- grass-addons/grass6/raster/r.connectivity.distance/Makefile (rev 0)
+++ grass-addons/grass6/raster/r.connectivity.distance/Makefile 2012-08-30 18:27:59 UTC (rev 52995)
@@ -0,0 +1,7 @@
+MODULE_TOPDIR = ../..
+
+PGM = r.connectivity.distance
+
+include $(MODULE_TOPDIR)/include/Make/Script.make
+
+default: script
Property changes on: grass-addons/grass6/raster/r.connectivity.distance/Makefile
___________________________________________________________________
Added: svn:executable
+ *
Added: grass-addons/grass6/raster/r.connectivity.distance/description.html
===================================================================
--- grass-addons/grass6/raster/r.connectivity.distance/description.html (rev 0)
+++ grass-addons/grass6/raster/r.connectivity.distance/description.html 2012-08-30 18:27:59 UTC (rev 52995)
@@ -0,0 +1,118 @@
+<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN">
+<html>
+<head>
+<title>GRASS GIS manual: r.connectivity.distance</title>
+<meta http-equiv="Content-Type" content="text/html; charset=iso-8859-1">
+<link rel="stylesheet" href="grassdocs.css" type="text/css">
+</head>
+<body bgcolor="white">
+
+<img src="grass_logo.png" alt="GRASS logo"><hr align=center size=6 noshade>
+
+<h2>NAME</h2>
+<em><b>r.connectivity.distance</b></em> - Compute cost-distance between all polygons (patches) of an input vector map within a user defined euclidean distance threshold
+<h2>KEYWORDS</h2>
+<dl>
+raster, vector, connectivity, distance, network
+</dl>
+
+<h2>SYNOPSIS</h2>
+<b>r.connectivity.distance</b><br>
+<b>r.connectivity.distance help</b><br>
+<b>r.connectivity.distance</b> [-<b>ste</b>] <b>patches</b>=<em>patches (input)</em> [<b>costs</b>=<em>costs (input)</em>] <b>pop_proxy</b>=<em>pop_proxy</em> <b>prefix</b>=<em>string</em> <b>folder</b>=<em>string</em> [<b>cutoff</b>=<em>float</em>] [<b>border_dist</b>=<em>integer</em>] [--<b>verbose</b>] [--<b>quiet</b>]
+
+<h3>Flags:</h3>
+<dl>
+<dt><b>-s</b></dt>
+<dd>Extract and save short paths and closest points into a vector map</dd>
+
+<dt><b>-t</b></dt>
+<dd>Rasterize patches with "all-touched" option using GDAL (can useful for narrow patches). This can be time consuming.</dd>
+
+<dt><b>-e</b></dt>
+<dd>Use euclidean distance (not cost distance)</dd>
+
+<dt><b>--verbose</b></dt>
+<dd>Verbose module output</dd>
+<dt><b>--quiet</b></dt>
+<dd>Quiet module output</dd>
+</dl>
+
+<h3>Parameters:</h3>
+<dl>
+<dt><b>folder</b>=<em>string</em></dt>
+<dd>Folder where all (non map) output from r.connectivity.* is stored (output)</dd>
+
+<dt><b>patches</b>=<em>string</em></dt>
+<dd>Name of input vector map containing patches (input)</dd>
+
+<dt><b>costs</b>=<em>string</em></dt>
+<dd>Name of input costs raster map (input)</dd>
+
+<dt><b>pop_proxy</b>=<em>string</em></dt>
+<dd>Column containig proxy for population size (attribute values must be numeric and neither values <= 0 nor NULL allowed!)</dd>
+
+<dt><b>prefix</b>=<em>string</em></dt>
+<dd>Prefix for output maps of the r.connectivity.*-tools. Also temporary data uses this prefix.</dd>
+
+<dt><b>cutoff</b>=<em>float</em></dt>
+<dd>Maximum euclidean search distance around patches in map units</dd>
+<dd>Default: <em>10000</em></dd>
+
+<dt><b>border_dist</b>=<em>integer</em></dt>
+<dd>Number of border cells used for distance measuring</dd>
+<dd>Default: <em>50</em></dd>
+</dl>
+
+<h2>DESCRIPTION</h2>
+
+Recently, graph-theory has been characterised as an efficient and useful tool for conservation planning (e.g. Bunn et al. 2000, Calabrese & Fagan 2004, Minor et al. 2008, Zetterberg et. al. 2010). As a part of the r.connectivity.* tool-chain, r.connectivity.distance is intended to make graph-theory more easily available to conservation planning.
+<br><br>
+r.connectivity.distance is the first tool of the r.connectivity.*-toolchain (followed by r.connectivity.network and r.connectivity.corridor).
+<br><br>
+r.connectivity.distance loops through all polygons in the input vector map and calculates the cost-distance to all the other polygons within a user-defined defined Euclidean distance threshold. It produces three csv-files containing an edge-list (connections between patches) and a vertex-list (patches) in two parts.
+<br><br>
+The edge-list csv-file structure (edges.csv) looks like this:
+<dl>
+<dd><em><b>from_patch,to_patch,cost-distance</b></em></dd>
+<dd><em>1,2,150</em></dd>
+<dd><em>1,3,75</em></dd>
+<dd><em>...</em></dd>
+</dl>
+<dl>
+The vertex-list consists of two parts:
+<dd>- one part containing the coordinates of the centroids of the input polygons (only for visualisation purposes (not for analysis) ("vertices_part_1.csv") and</dd>
+<dd>- one part containing a population proxy (values of an attribute column from the input patch vector map (both specified by the user)), representing a proxy for the amount of organisms potentially dispersing from a patch (e.g. habitat area) ("vertices_part_2.csv").</dd>
+</dl><br>
+The structure of the vertices_part_1.csv looks like this:
+<dl>
+<dd><em><b>patch_id;cetroid_x_coordinate;cetroid_y_coordinate</b></em></dd>
+<dd><em>1;796973.455617378;7874304.58988322</em></dd>
+<dd><em>2;1028575.90731092;7909142.41545215</em></dd>
+<dd><em>...</em></dd>
+</dl>
+<br>
+The structure of the vertices_part_2.csv looks like this:
+<dl>
+<dd><em><b>patch_id;population_proxy</b></em></dd>
+<dd><em>1;20.7652</em></dd>
+<dd><em>2;74.0114</em></dd>
+<dd><em>...</em></dd>
+</dl>
+<br>
+The map output of r.connectivity.distance is a (cost) distance raster map for every input polygon which later on are used in r.connectivity.corridors (together with output from r.connectivity.network) for corridor identification. All non map output is stored in a folder, which will be used throughout the entire r.connectivity toolchain. Output-files are: vertices_part_1.csv, vertices_part_2.csv, edges.csv and a log-file (r_connectivity.log). The log-file is used to store and document settings of the tools and to pass arguments to the following tools of the toolchain.
+<br><br>
+Distance is measured as border to border distance. The user can define the number of cells (n) along the border to be used for distance measuring. The distance from a (start) polygon to another (end) is measured as the n-th closest cell on the border of the other (end) polygon. An increased number of border cells used for distance measuring affects (increases) also the width of possible corridors computed with r.connectivity.corridor.
+
+
+<h2>SEE ALSO</h2>
+<a href="r.connectivity.network">r.connectivity.network</a>, <a href="r.connectivity.corridors">r.connectivity.corridors</a><br>
+
+
+<h2>AUTHOR</h2>
+Stefan Blumentrath, Norwegian Institute for Nature Research (NINA)
+
+<i>Last changed: $Date$</i>
+
+</body>
+</html>
Property changes on: grass-addons/grass6/raster/r.connectivity.distance/description.html
___________________________________________________________________
Added: svn:executable
+ *
Added: grass-addons/grass6/raster/r.connectivity.distance/r.connectivity.distance
===================================================================
--- grass-addons/grass6/raster/r.connectivity.distance/r.connectivity.distance (rev 0)
+++ grass-addons/grass6/raster/r.connectivity.distance/r.connectivity.distance 2012-08-30 18:27:59 UTC (rev 52995)
@@ -0,0 +1,555 @@
+#!/bin/sh
+#
+############################################################################
+#
+# MODULE: r.connectivity.distance
+# AUTHOR(S): Stefan Blumentrath <stefan dot blumentrath at nina dot no >
+# PURPOSE: Compute cost-distance between all polygons (patches) of an
+# input vector map within a user defined euclidean distance threshold
+#
+# Recently, graph-theory has been characterised as an efficient and
+# useful tool for conservation planning (e.g. Bunn et al. 2000,
+# Calabrese & Fagan 2004, Minor & Urban 2008, Zetterberg et. al. 2010).
+# As a part of the r.connectivity.* tool-chain, r.connectivity.distance
+# is intended to make graph-theory more easily available to conservation
+# planning.
+#
+# r.connectivity.distance is the first tool of the r.connectivity.*-toolchain
+# (followed by r.connectivity.network and r.connectivity.corridor).
+#
+# r.connectivity.distance loops through all polygons in the input vector
+# map and calculates the cost-distance to all the other polygons within
+# a user-defined defined Euclidean distance threshold. It produces three
+# csv-files containing an edge-list (connections between patches) and a
+# vertex-list (patches) in two parts.
+#
+# The edge-list csv-file structure (edges.csv) looks like this:
+#
+# from_patch,to_patch,cost-distance
+# 1,2,150
+# 1,3,75
+# ...
+#
+# The vertex-list consists of two parts:
+# - one part containing the coordinates of the centroids of
+# the input polygons (only for visualisation purposes (not for
+# analysis) ("vertices_part_1.csv") and
+# - one part containing a population proxy (values of an
+# attribute column from the input patch vector map (both
+# specified by the user)), representing a proxy for the
+# amount of organisms potentially dispersing from a patch
+# (e.g. habitat area) ("vertices_part_2.csv").
+#
+# The structure of the vertices_part_1.csv looks like this:
+#
+# patch_id;cetroid_x_coordinate;cetroid_y_coordinate
+# 1;796973.455617378;7874304.58988322
+# 2;1028575.90731092;7909142.41545215
+# ...
+#
+# The structure of the vertices_part_2.csv looks like this:
+#
+# patch_id;population_proxy
+# 1;20.7652
+# 2;74.0114
+# ...
+#
+# The map output of r.connectivity.distance is a (cost)
+# distance raster map for every input polygon which later on
+# are used in r.connectivity.corridors (together with output
+# from r.connectivity.network) for corridor identification.
+# All non map output is stored in a folder, which will be used
+# throughout the entire r.connectivity toolchain. Output-files
+# are: vertices_part_1.csv, vertices_part_2.csv, edges.csv and
+# a log-file (r_connectivity.log). The log-file is used to
+# store and document settings of the tools and to pass
+# arguments to the following tools of the toolchain.
+#
+# Distance is measured as border to border distance. The user
+# can define the number of cells (n) along the border to be
+# used for distance measuring. The distance from a (start)
+# polygon to another (end) is measured as the n-th closest
+# cell on the border of the other (end) polygon. An increased
+# number of border cells used for distance measuring affects
+# (increases) also the width of possible corridors computed
+# with r.connectivity.corridor.
+#
+# COPYRIGHT: (C) 2012 by the Norwegian Institute for Nature Research (NINA)
+#
+# This program is free software under the GNU General Public
+# License (>=v2). Read the file COPYING that comes with GRASS
+# for details.
+#
+#############################################################################
+#
+# REQUIREMENTS:
+# awk
+#
+#
+#%Module
+#% description: Compute cost-distance between all patches of an input vector map
+#%End
+#
+#%Option
+#% key: patches
+#% type: string
+#% required: yes
+#% multiple: no
+#% key_desc: patches (input)
+#% description: Name of input vector map containing patches
+#% gisprompt: old,vector,vector
+#%End
+#
+#%Option
+#% key: costs
+#% type: string
+#% required: no
+#% multiple: no
+#% key_desc: costs (input)
+#% description: Name of input costs raster map
+#% gisprompt: old,cell,raster
+#%End
+#
+#%Option
+#% key: pop_proxy
+#% type: string
+#% required: yes
+#% multiple: no
+#% key_desc: pop_proxy
+#% description: Column containig proxy for population size (not NULL and > 0)
+#%End
+#
+#%option
+#% key: prefix
+#% type: string
+#% description: Prefix for output maps of the r.connectivity.*-tools
+#% required : yes
+#%end
+#
+#%option
+#% key: folder
+#% type: string
+#% description: Folder where all (non map) output from r.connectivity.* is stored
+#% required : yes
+#%end
+#
+#%option
+#% key: cutoff
+#% type: double
+#% description: Maximum euclidean search distance around patches in map units
+#% required : no
+#% guisection: Settings
+#% answer : 10000
+#%end
+#
+#%option
+#% key: border_dist
+#% type: integer
+#% description: Number of border cells used for distance measuring
+#% required : no
+#% guisection: Settings
+#% answer : 50
+#%end
+#
+#%flag
+#% key: s
+#% description: Extract and save short paths and closest points into a vector map
+#% guisection: Settings
+#%end
+#
+#%flag
+#% key: t
+#% description: Rasterize patches with "all-touched" option using GDAL
+#% guisection: Settings
+#%end
+#
+#%flag
+#% key: e
+#% description: Use euclidean distance (not cost distance)
+#% guisection: Settings
+#%end
+#
+
+#Input:
+PATCHES="${GIS_OPT_PATCHES}"
+POP_PROXY="${GIS_OPT_POP_PROXY}"
+COSTS="${GIS_OPT_COSTS}"
+CUTOFF="${GIS_OPT_CUTOFF}"
+BORDER_DIST="${GIS_OPT_BORDER_DIST}"
+
+#Output
+PREFIX="${GIS_OPT_PREFIX}"
+FOLDER="${GIS_OPT_FOLDER}"
+
+#Settings
+S_FLAG=${GIS_FLAG_S}
+T_FLAG=${GIS_FLAG_T}
+E_FLAG=${GIS_FLAG_E}
+#
+#Check if script is started from within GRASS
+if [ -z "$GISBASE" ] ; then
+ echo "You must be in GRASS GIS to run this program." 1>&2
+ exit 1
+fi
+#
+# Pass command line arguments to gui and start it
+if [ "$1" != "@ARGS_PARSED@" ] ; then
+ exec g.parser "$0" "$@"
+fi
+#
+
+#Check if awk is installed
+if [ ! -x "`which awk`" ] ; then
+ g.message -e "awk is required, please install awk first"
+ exit 1
+fi
+
+#set environment so that awk works properly in all languages
+unset LC_ALL
+LC_NUMERIC=C
+export LC_NUMERIC
+
+#Check if patch vector map exists
+eval `g.findfile element=vector file=${PATCHES}`
+patchfile=$file
+if [ ! -r "$file" ] ; then
+ g.message -e "Cannot find vector map ${PATCHES}."
+ exit 1
+fi
+
+if [ $E_FLAG -ne 1 ] ; then
+ #Check if cost raster map exists
+ eval `g.findfile element=cell file=${COSTS}`
+ if [ ! -r "$file" ] ; then
+ g.message -e "Cannot find raster map ${COSTS}."
+ exit 1
+ fi
+fi
+
+#Check if cat column exists
+cat_column=$(v.db.connect -c map="${PATCHES}" | grep "cat" | cut -f2 -d'|')
+if [ ! "$cat_column" ] ; then
+ g.message -e "Cannot find the reqired column cat in vector map ${PATCHES}."
+ exit 1
+fi
+
+#Check if pop_proxy column exists
+pop_proxy_column=$(v.db.connect -c map="${PATCHES}" | grep "$POP_PROXY" | cut -f2 -d'|')
+if [ ! "$pop_proxy_column" ] ; then
+ g.message -e "Cannot find column ${POP_PROXY} in vector map ${PATCHES}."
+ exit 1
+fi
+
+#Check if pop_proxy column is numeric type
+pop_proxy_type_check=$(v.db.connect -c map="${PATCHES}" | grep "$POP_PROXY" | cut -f1 -d'|')
+if [ "$pop_proxy_type_check" != "INTEGER" -a "$pop_proxy_type_check" != "DOUBLE PRECISION" ] ; then
+ g.message -e "Column ${POP_PROXY} is of type $pop_proxy_type_check. Only numeric types (integer or double precision) allowed!"
+ exit 1
+fi
+
+#Check if pop_proxy column does not contain values <= 0
+pop_proxy_check=$(v.db.select -c map="${PATCHES}" columns="cat,$POP_PROXY" fs=' ' nv=-9999 | awk '{if($2 <= 0) print "Patch with category " $1 " has a population proxy value of " $2 ". Neither values <= 0 nor NULL allowed!"}')
+if [ -n "$pop_proxy_check" ] ; then
+ g.message -e "$pop_proxy_check"
+ exit 1
+fi
+
+#Check if FOLDER exists
+if [ ! -d "$FOLDER" ] ; then
+ g.message -w "Cannot find folder ${FOLDER}. This folder will be created now."
+ mkdir "$FOLDER"
+ if [ ! -d "$FOLDER" ] ; then
+ g.message -e "Could not create folder ${FOLDER}..."
+ exit 1
+ fi
+fi
+
+#pop_proxy
+# pop_proxy column has to be numeric, and neither values <= 0 nor NULL are allowed!)
+#
+# Pass command line arguments to log-file
+cat << EOF > ${FOLDER}/r_connectivity.log
+PREFIX=$PREFIX
+PATCHES=$PATCHES
+POP_PROXY="${GIS_OPT_POP_PROXY}"
+COSTS=$COSTS
+CUTOFF=$CUTOFF
+BORDER_DIST="${GIS_OPT_BORDER_DIST}"
+PREFIX=$PREFIX
+FOLDER=$FOLDER
+EUCL_DIST=$E_FLAG
+EOF
+
+
+########################################################################
+
+eval `g.region -ugp rast=$COSTS align=$COSTS`
+max_n=$n
+min_s=$s
+max_e=$e
+min_w=$w
+cost_nsres=$nsres
+cost_ewres=$ewres
+
+
+
+#Init network data
+#Remove existing network files (Only if networkdata should be produced completly new)
+echo "from_patch;to_patch;cost-distance" > ${FOLDER}/edges.csv
+echo "patch_id;cetroid_x_coordinate;cetroid_y_coordinate" > ${FOLDER}/vertices_part_1.csv
+echo "patch_id;population_proxy" > ${FOLDER}/vertices_part_2.csv
+
+#Get patch IDs (cat column)
+patch_list=$(v.db.select -c map=$PATCHES columns=cat)
+
+###Prepare boundaries for measuring boundary to boundary distance between patches
+g.message -v "Preparing boundaries..."
+
+###Extract centroids and export data as network vertices (output will not contain centroids of multipolygons)
+#Write centroids to nodes input text file to R (igraph)
+v.to.db -p --quiet map=$PATCHES option=coor | cut -f1,2,3 -d'|' --output-delimiter=';' >> ${FOLDER}/vertices_part_1.csv
+#Write centroids to nodes input text file to R (igraph)
+v.db.select -c map=$PATCHES columns="cat,${POP_PROXY}" fs=";" >> ${FOLDER}/vertices_part_2.csv
+
+#Rasterize patches
+if [ $T_FLAG -eq 1 ] ; then
+ #Rasterize patches with "all-touched" mode using GDAL
+ #Read region-settings
+ eval `g.region -ugp`
+
+ #Check if GDAL-GRASS plugin is installed
+ gdal_plugin=$(gdalinfo /home/stefan/grassdata/Vernevaluering/PERMANENT/vector/vernet_fjell_norge_sverige 2>&1 > /dev/null | head -n 1 | cut -f1 -d' ')
+
+ if [ "$gdal_plugin" != "ERROR" ] ; then
+ #With GDAL-GRASS plugin
+ #Locate file for patch vector map
+ eval `g.findfile element=vector file=${PATCHES}`
+ #Rasterize vector map with all-touched option
+ gdal_rasterize -l 1 -at -tr $ewres $nsres -te $w $s $e $n -ot Uint32 -a cat $file ${FOLDER}/patches_rast.tif -q
+ #
+ else
+ #Without GDAL-GRASS-plugin
+ g.message -w "Cannot find GDAL-GRASS plugin. Consider installing it in order to save time for all-touched rasterisation"
+ #Export patch vector map to temp-file in a GDAL-readable format (shp)
+ v.out.ogr -e -c --quiet input=$PATCHES type=area dsn=${FOLDER}/patches_vect.shp
+ #Rasterize vector map with all-touched option
+ gdal_rasterize -l patches_vect -at -tr $ewres $nsres -te $w $s $e $n -ot Uint32 -a cat ${FOLDER}/patches_vect.shp ${FOLDER}/patches_rast.tif -q
+ #Remove vector temp-file
+ rm -f ${FOLDER}/patches_vect.shp
+ rm -f ${FOLDER}/patches_vect.shx
+ rm -f ${FOLDER}/patches_vect.dbf
+ rm -f ${FOLDER}/patches_vect.prj
+ #Import rasterized containing patches
+ r.in.gdal -o --overwrite --quiet input=${FOLDER}/patches_rast.tif output=${PREFIX}_patches_pol
+ #Remove temporary geotiff
+ rm -f ${FOLDER}/patches_rast.tif
+ fi
+else
+ #Simple rasterisation (only area)
+ v.to.rast --quiet --overwrite input=$PATCHES output=${PREFIX}_patches_pol use=cat type=area
+fi
+#Extract boundaries from patch raster map
+r.mapcalc "${PREFIX}_patches_boundary=if(${PREFIX}_patches_pol,\
+ if((\
+ (${PREFIX}_patches_pol[-1,0]!=${PREFIX}_patches_pol)|||\
+ (${PREFIX}_patches_pol[0,1]!=${PREFIX}_patches_pol)|||\
+ (${PREFIX}_patches_pol[1,0]!=${PREFIX}_patches_pol)|||\
+ (${PREFIX}_patches_pol[0,-1]!=${PREFIX}_patches_pol)),${PREFIX}_patches_pol,null()),\
+ null())"
+
+#Init output vector maps if they are requested by user
+if [ $S_FLAG -eq 1 ] ; then
+ v.edit --quiet map=${PREFIX}_closest_stop_points tool=create
+ v.db.addtable --quiet map=${PREFIX}_closest_stop_points columns="cat integer, X integer,Y integer,FROM_P integer,TO_P integer,DIST double precision"
+
+ v.edit map=${PREFIX}_shortest_paths tool=create
+ v.db.addtable --quiet map=${PREFIX}_shortest_paths columns="cat integer, value integer, label varchar (10)"
+fi
+
+###Loop through patches
+for p in $patch_list
+do
+ g.message -v "Calculating connectivity-distances for patch number $p"
+ #Extract patch polygons
+ centroid_number=$(v.extract --overwrite --verbose input=$PATCHES output=${PREFIX}_patch_${p} where="cat=${p}" 2>&1 | grep "Number of centroids:" | cut -f4 -d' ')
+ #Extract centroids of convex hull for multipolygons and export data as network vertices
+ if [ $centroid_number -gt 1 ] ; then
+ v.hull -a -f --overwrite --quiet input=${PREFIX}_patch_${p} output=${PREFIX}_patch_${p}_ch 2>&1 > /dev/null
+ xy=$(v.to.db -p --quiet map=${PREFIX}_patch_${p}_ch option=coor | cut -f2,3 -d'|' --output-delimiter=';')
+ echo "${p};$xy" >> ${FOLDER}/vertices_part_1.csv
+ g.remove -f vect=${PREFIX}_patch_${p}_ch --quiet 2>&1 > /dev/null
+ g.message -v "Patch ${p} has $centroid_number centroids. Using centroid of the convex hull of all polygon parts."
+ #Remove node data for patches without centroid
+ elif [ $centroid_number -eq 0 ] ; then
+ g.remove -f vect=${PREFIX}_patch_${p} --quiet 2>&1 > /dev/null
+ echo "${p};no centroid" >> ${FOLDER}/unconsidered_patches.csv
+ sed -i "/^${p};/d" ${FOLDER}/vertices_part_1.csv
+ sed -i "/^${p};/d" ${FOLDER}/vertices_part_2.csv
+ g.message -w "Patch ${p} has no centroid and is therefor not considered in anaysis. Consider cleaning your input data..."
+ continue
+ fi
+
+ #Set region to start-patch
+ g.region -u --quiet vect=${PREFIX}_patch_${p} n=n+$cost_nsres s=s-$cost_nsres e=e+$cost_ewres w=w-$cost_ewres align=$COSTS save="${PREFIX}_${p}" --overwrite
+ WIND_OVERRIDE="${PREFIX}_${p}"
+ export WIND_OVERRIDE
+
+ #Prepare start patch
+ r.mapcalc "${PREFIX}_patch_${p}=if(${PREFIX}_patches_boundary==${p},1,null())"
+ #unset WIND_OVERRIDE
+
+ #Check if patch was rasterised (patches smaller raster resolution and close to larger patches may not be rasterised)
+ eval `r.info -r map=${PREFIX}_patch_${p}`
+ if [ "$min" = "NULL" ] ; then
+ g.remove -f vect=${PREFIX}_patch_${p} rast=${PREFIX}_patch_${p} --quiet 2>&1 > /dev/null
+ echo "${p};not rasterised" >> ${FOLDER}/unconsidered_patches.csv
+ sed -i "/^${p};/d" ${FOLDER}/vertices_part_1.csv
+ sed -i "/^${p};/d" ${FOLDER}/vertices_part_2.csv
+ g.message -w "Patch ${p} was not rasterised and is therefor not considered in anaysis. Consider adjusting region resolution..."
+ continue
+ fi
+
+ #Prepare stop patches
+ ############################################
+ #awk
+ ############################################
+ eval `g.region -ugp --quiet n=n+$CUTOFF s=s-$CUTOFF e=e+$CUTOFF w=w-$CUTOFF align=$COSTS`
+ n_test=$(echo $n $max_n | awk '{if($1 < $2) print $2}')
+ if [ $n_test ] ; then
+ n=$n
+ else
+ n=$max_n
+ fi
+ s_test=$(echo $s $min_s | awk '{if($1 > $2) print $2}')
+ if [ $s_test ] ; then
+ s=$s
+ else
+ s=$min_s
+ fi
+ e_test=$(echo $e $max_e | awk '{if($1 < $2) print $2}')
+ if [ $e_test ] ; then
+ e=$e
+ else
+ e=$max_e
+ fi
+ w_test=$(echo $w $min_w | awk '{if($1 > $2) print $2}')
+ if [ $w_test ] ; then
+ w=$w
+ else
+ w=$min_w
+ fi
+
+ g.region -u --quiet n=$n s=$s e=$e w=$w align=$COSTS save="${PREFIX}_${p}_buffer" --overwrite
+ WIND_OVERRIDE="${PREFIX}_${p}_buffer"
+ export WIND_OVERRIDE
+
+ if [ $E_FLAG -eq 1 ] ; then
+ r.buffer --overwrite --quiet input=${PREFIX}_patch_${p} output=MASK distances=$CUTOFF
+ r.mapcalc "${PREFIX}_patch_${p}_neighbours_contur=if(${PREFIX}_patches_boundary==${p},null(),${PREFIX}_patches_boundary)"
+ g.remove -f rast=MASK --quiet
+
+ #Calculate euclidean distance
+ r.grow.distance --overwrite --verbose input=${PREFIX}_patch_${p} distance=${PREFIX}_patch_${p}_eucl_dist
+
+ #Export distance at boundaries
+ r.stats -1 -n -g --quiet input=${PREFIX}_patch_${p}_neighbours_contur,${PREFIX}_patch_${p}_eucl_dist fs=";" | sort -n -t ';' -k 3 -k 4 > ${FOLDER}/tmp
+
+ #Remove temporary map data for patch
+ g.remove -f vect=${PREFIX}_patch_${p} rast=${PREFIX}_patch_${p},${PREFIX}_patch_${p}_neighbours_contur --quiet 2>&1 > /dev/null
+
+ else
+ #Create buffer around start-patch as a mask for cost distance analysis
+ r.buffer --overwrite --quiet input=${PREFIX}_patch_${p} output=MASK distances=$CUTOFF
+ r.mapcalc "${PREFIX}_patch_${p}_neighbours_contur=if(${PREFIX}_patches_boundary==${p},null(),${PREFIX}_patches_boundary)"
+ g.remove -f rast=MASK --quiet
+
+ #Calculate cost distance
+ r.cost -k -n --overwrite --quiet input=$COSTS output="${PREFIX}_patch_${p}_cost_dist" start_rast="${PREFIX}_patch_${p}"
+
+ #Export distance at boundaries
+ r.stats -1 -n -g --quiet input=${PREFIX}_patch_${p}_neighbours_contur,${PREFIX}_patch_${p}_cost_dist fs=";" | sort -n -t ';' -k 3 -k 4 > ${FOLDER}/tmp
+
+ #Remove temporary map data for patch
+ g.remove -f vect=${PREFIX}_patch_${p} rast=${PREFIX}_patch_${p},${PREFIX}_patch_${p}_neighbours_contur --quiet 2>&1 > /dev/null
+
+ fi
+ #Loop through connections from the start-patch
+ connections=$(cat ${FOLDER}/tmp | cut -f3 -d';' | sort -n -u)
+
+ for c in $connections
+ do
+ #Find closest points on neigbours patches to temporary file
+ cat ${FOLDER}/tmp | awk -v patch=$c -v start=$p -F ";" '{ if($3 == patch) print $1 ";" $2 ";" start ";" $3 ";" $4}' > ${FOLDER}/tmp_part
+
+ #Save closest points on neighbour patches to temporary file
+ cat ${FOLDER}/tmp_part | head -n 1 | cut -f1-5 -d';' >> ${FOLDER}/closest.points
+
+ #Save edges to network dataset
+ #order_length=$(cat ${FOLDER}/tmp_part | wc -l)
+ #relevant_border_length=`expr $border_length \* $BORDER_DIST / 100`
+ if [ $BORDER_DIST -le 1 ] ; then
+ cat ${FOLDER}/tmp_part | awk -F ";" '{ print $3 ";" $4 ";" $5 }' | head -n 1 | tail -n 1 >> ${FOLDER}/edges.csv
+ else
+ cat ${FOLDER}/tmp_part | awk -F ";" '{ print $3 ";" $4 ";" $5 }' | head -n $BORDER_DIST | tail -n 1 >> ${FOLDER}/edges.csv
+ fi
+ done
+
+
+ #Save closest points and shortest paths through cost raster as vector map (r.drain limited to 1000 points) if requested by user
+ if [ $S_FLAG -eq 1 ] ; then
+ points_n=$(cat ${FOLDER}/closest.points | wc -l)
+ if [ $points_n -gt 1000 ] ; then
+ tiles=`expr $points_n / 1000`
+ rest=`expr $points_n % 1000`
+ if [ $rest -ne 0 ] ; then
+ tiles=`expr $tiles + 1`
+ fi
+ n=0
+ #Init closest points file for start-patch
+ v.edit --quiet map=patch_${p}_closest_points tool=create
+ v.db.addtable --quiet map=patch_${p}_closest_points columns="cat integer, X integer,Y integer,FROM_P integer,TO_P integer,DIST double precision"
+ #Init cost paths file for start-patch
+ v.edit --quiet map="patch_${p}_cost_pathes" tool=create
+ v.db.addtable --quiet map="patch_${p}_cost_pathes" columns="cat integer, value integer, label varchar (10)"
+ while [ $n -lt $tiles ]
+ do
+ start=`expr $n \* 1000`
+ tail=`expr $points_n - $start`
+ n=`expr $n + 1`
+ cat ${FOLDER}/closest.points | tail -n $tail | head -n 1000 > closest.points_part
+ #Import closest points for start-patch in 1000er blocks
+ v.in.ascii -n -r --overwrite --quiet input=${FOLDER}/closest.points_part output=patch_${p}_closest_points_part_${n} fs=";" columns="X integer,Y integer,FROM_P integer,TO_P integer,DIST double precision"
+ v.patch -a -e --overwrite --quiet input=patch_${p}_closest_points_part_${n} output=patch_${p}_closest_points
+ #Extract shortest paths for start-patch in 1000er blocks
+ r.drain --overwrite --quiet input="patch_${p}_cost_dist" output="patch_${p}_cost_pathes_part_${n}" vector_points="patch_${p}_closest_points_part_${n}"
+ r.to.vect --overwrite --quiet input="patch_${p}_cost_pathes_part_${n}" output="patch_${p}_cost_pathes_part_${n}" feature=line
+ v.patch -a -e --overwrite --quiet input="patch_${p}_cost_pathes_part_${n}" output="patch_${p}_cost_pathes"
+ done
+ #Patch closest points for start-patch to complete closest points map
+ v.patch -a -e --overwrite --quiet input=patch_${p}_closest_points output=${PREFIX}_closest_stop_points
+ #Patch shortest paths for start-patch to complete shortest paths map
+ v.patch -a -e --overwrite --quiet input="patch_${p}_cost_pathes" output=${PREFIX}_shortest_paths
+ #Remove temporary map data
+ g.mremove --quiet -f rast=*_part_${n} vect=*_part_${n}
+ rm -f ${FOLDER}/closest.points_part
+ else
+ #Import closest points for start-patch
+ v.in.ascii -n -r --overwrite --quiet input=${FOLDER}/closest.points output=patch_${p}_closest_points fs=";" columns="X integer,Y integer,FROM_P integer,TO_P integer,DIST double precision"
+ #Patch closest points for start-patch to complete closest points map
+ v.patch -a -e --overwrite --quiet input=patch_${p}_closest_points output=${PREFIX}_closest_stop_points
+ #Extract shortest paths for start-patch
+ r.drain --overwrite --quiet input="patch_${p}_cost_dist" output="patch_${p}_cost_pathes" vector_points="patch_${p}_closest_points"
+ r.to.vect --overwrite --quiet input="patch_${p}_cost_pathes" output="patch_${p}_cost_pathes" feature=line
+ #Patch shortest paths for start-patch to complete shortest paths map
+ v.patch -a -e --overwrite --quiet input="patch_${p}_cost_pathes" output=${PREFIX}_shortest_paths
+ fi
+ fi
+ #Reset region to original settings and remove saved temporary regions
+ unset WIND_OVERRIDE
+ g.remove region="${PREFIX}_${p},${PREFIX}_${p}_buffer" --q
+
+done
+
+#Remove external temporary data
+rm -f ${FOLDER}/tmp
+rm -f ${FOLDER}/tmp_part
Property changes on: grass-addons/grass6/raster/r.connectivity.distance/r.connectivity.distance
___________________________________________________________________
Added: svn:executable
+ *
More information about the grass-commit
mailing list