Hi,
I try to create a tiny script that allows to create a focal filter based on
r.cost. This means that rather using a moving window I use r.cost to create a
"buffer" around a cell which is then used for getting certain statistics of the
neighborhood of that cell (e.g. mean, max, sum) and create a new raster based
on these statistical values. I'd like to use r.cost as it also allows to create
distance buffers along rasterized lines (the main purpose for me). One example
can be: "Get the maximum raster value which is up- and downstream from a cell
in a river distance of 500 m"...
So far I think I have a plan how it can work and I tried to write the script in
python. Anyway there are some problems where I am not sure how I can solve
them...
1) I don't really know how I can loop over the points (raster cells) and update
the points with the statistical value (e.g. mean). I try following:
grass.write_command("db.execute", stdin = 'UPDATE point_tmp%s SET univar=%d
WHERE cat = %d' % (str(os.getpid()),stat,cat))
but I get a "time out" (WARNING: Busy SQLITE db, already waiting for 10
seconds, 20 sec...). Probably it is related to the fact the two processes are
using the same database table (loop over the points and update)...
2) I don't know what is the best way to get a variable based on the input,
which is like defining a global variable in a if-loop, but I am not sure if
that is a good way?
Best is to look at the script which is attached. Theoretically it should be
possible to be used on every raster. But as it loops over each cell it takes
really long time to process large maps.
Maybe anyone wants to help me further developing this tiny python script. Any
help is very much appreciated. Thank you!
Cheers
/johannes
--
NEU: FreePhone 3-fach-Flat mit kostenlosem Smartphone!
Jetzt informieren: http://mobile.1und1.de/?ac=OM.PW.PW003K20328T7073a
#!/usr/bin/env python
#
############################################################################
#
# SCRIPT: Creating focal statistical raster based on r.cost distance
#
# AUTHOR(S): Johannes Radinger
# DATE: 2012-05-03
#
#############################################################################
#%Module
#% description: Calculating focal raster based on a r.cost distance of a raster cell
#% keywords: raster
#%End
#%option
#% key: input
#% type: string
#% gisprompt: old,cell,raster
#% description: input raster
#% required: yes
#%end
#%option
#% key: stat
#% type: string
#% multiple:no
#% options:mean,max,sum
#% description: Desired focal statistic for output raster
#% required: yes
#%end
#%option
#% key: distance
#% type: double
#% required: yes
#% multiple: no
#% description: Distance for filter/buffer
#%End
import sys
import os
import atexit
import sqlite3
import grass.script as grass
tmp_map_rast = None
tmp_map_vect = None
def cleanup():
if (tmp_map_rast or tmp_map_vect):
grass.run_command("g.remove",
rast = tmp_map_rast,
vect = tmp_map_vect,
quiet = True)
def main():
#defining temporary maps
tmp_map_rast = ["MASK","tmp"+str(os.getpid()),"distance_mask"+str(os.getpid())]
tmp_map_vect = ["point_tmp"+str(os.getpid())]
#database-connection
env = grass.gisenv()
gisdbase = env['GISDBASE']
location = env['LOCATION_NAME']
mapset = env['MAPSET']
grass.run_command("db.connect",
driver = "sqlite",
database = os.path.join(gisdbase, location, mapset,'sqlite.db'))
database = sqlite3.connect(os.path.join(gisdbase, location, mapset, 'sqlite.db'))
db = database.cursor()
# Getting resolution of raster map
res = int(grass.read_command("g.region",
flags = "m").strip().split('\n')[4].split('=')[1])
# Defining variables
input = [options['input']]
distance = options['distance']
# Remove any mask before processing
grass.run_command("g.remove",
rast = "MASK")
for i in input:
grass.mapcalc("$tmp = if($raster,$res,null())", # creating temporary map
tmp = "tmp"+str(os.getpid()),
res = res,
raster = i)
# make points from raster
grass.run_command("r.to.vect",
overwrite = True,
input = "tmp"+str(os.getpid()),
output = "point_tmp"+str(os.getpid()),
feature = "point")
# Get coordinates for each point
grass.run_command("v.db.addcol",
map = "point_tmp"+str(os.getpid()),
columns = "X DOUBLE, Y DOUBLE, univar DOUBLE")
grass.run_command("v.to.db",
map = "point_tmp"+str(os.getpid()),
type = "point",
option = "coor",
columns = "X,Y")
#update database connection
db = database.cursor()
# loop over points
for cat,X,Y in db.execute('SELECT cat, X, Y FROM point_tmp%s' % str(os.getpid())):
coors = str(X)+","+str(Y)
# creating focal distance buffer as MASK (using r.cost and a distance/buffer)
grass.run_command("r.cost",
flags = 'n',
overwrite = True,
max_cost = distance,
input = "tmp"+str(os.getpid()),
output = "distance_mask"+str(os.getpid()),
coordinate = coors)
# creating MASK based on r.cost output
grass.mapcalc("$MASK= if(!isnull($distance_mask),1,null())",
MASK = "MASK",
distance_mask = "distance_mask"+str(os.getpid()))
# Getting focal statistics (e.g. mean) based on distance_mask
univar = grass.read_command("r.univar", map = i)
stat = float(univar.split('\n')[10].split(':')[1])
# Here are still problems to set the statistics conditional based on the input
# Getting different statistics based on input
#if str(options['stat']) is "sum":
# stat = float(univar.split('\n')[15].split(':')[1])
#elif str(options['stat']) is "mean":
# stat = float(univar.split('\n')[10].split(':')[1])
# print stat
#elif str(options['stat']) is "max":
# stat = float(univar.split('\n')[8].split(':')[1])
# print stat
#update point map with new stat value
grass.write_command("db.execute", stdin = 'UPDATE point_tmp%s SET univar=%d WHERE cat = %d' % (str(os.getpid()),stat,cat))
grass.run_command("g.remove",rast = ["MASK","distance_mask%d" % os.getpid()])
grass.run_command("v.to.rast",
input = "point_"+str(os.getpid()),
output = "d"+str(distance)+"_"+i,
use = "univar")
grass.run_command("g.remove",
rast = tmp_map_rast,
vect = tmp_map_vect)
return 0
if __name__ == "__main__":
options, flags = grass.parser()
atexit.register(cleanup)
sys.exit(main())
_______________________________________________
grass-user mailing list
[email protected]
http://lists.osgeo.org/mailman/listinfo/grass-user
