[GRASS-user] Python script: focal filter based on r.cost
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 =
Re: [GRASS-user] Report from ongoing GRASS GIS Community Sprint in Prague
Il 28/05/2012 21:17, Markus Metz ha scritto: not to java SEXTANTE. While on it (the blaming), I think it is time for Paolo to apologize for his rude [1] comments [2][3][4][5][6][7] earlier in the thread. I'm more than happy to apologize if I have hurt the sensibility of anybody. However, I think there is some misunderstanding, let's try to clarify: -SEXTANTE has two versions now (which are like 2 different projects) -The QGIS version is still unstable, but is a much better software (cleaner design, different things, etc) -It would be great to coordinate both versions, but it has to be done manually. currently, Victor Olaya is putting time on the QGIS version and not in the gvSIG one, so they might get uncoordinated, which Victor do not see as a terrible problem, anyway. -SEXTANTE for gvSIG has more algorithms and is more prepared for that, while SEXTANTE for QGIS has more backends, and is better suited for it. [2] http://osgeo-org.1560.n6.nabble.com/Report-from-ongoing-GRASS-GIS-Community-Sprint-in-Prague-tp4977136p4977385.html the link is wrong, the original source code is at http://code.google.com/p/sextante/source/browse/trunk/soft/bindings/qgis-plugin/src From the link I provided, you can easily find the source code at: http://plugins.qgis.org/plugins/sextante/ so it should be, according to Victor choice: http://sextante.googlecode.com/svn/trunk/ [3] http://osgeo-org.1560.n6.nabble.com/Report-from-ongoing-GRASS-GIS-Community-Sprint-in-Prague-tp4977136p4977411.html Paolo is appropriating python SEXTANTE as the QGIS Python plugin called sextante. See [2] I really do not understand this (see above). What is Python SEXTANTE? Anyway, sending a link and appropriating code seems different acts. [4] http://osgeo-org.1560.n6.nabble.com/Report-from-ongoing-GRASS-GIS-Community-Sprint-in-Prague-tp4977136p4977430.html Paolo is again appropriating python SEXTANTE (Please do nont send bug reports concerning it to other bugtrackers.) as the QGIS Python plugin called sextante. See [2] Bugs concerning the SEXTANTE QGIS Python plugin should be sent to: http://hub.qgis.org/projects/sextante/issues I do not see how this is appropriating anything. [5] http://osgeo-org.1560.n6.nabble.com/Report-from-ongoing-GRASS-GIS-Community-Sprint-in-Prague-tp4977136p4977473.html Paolo: I can correct you :) 1. Paolo was not addressed. 2. The issue has been resolved, no need for any further correction 3. I, not Paolo, gave an explanation and examples for the issue discussed (there is a complete re-write of SEXTANTE in Python). As I explained in a previous mail, Victor Olaya is the person who can best correct anybody on this. I took the liberty of talking on his behalf just because we are in close contact, and he has not subscribed to this list. [6] http://osgeo-org.1560.n6.nabble.com/Report-from-ongoing-GRASS-GIS-Community-Sprint-in-Prague-tp4977136p4977489.html Paolo: ... and join us. Paolo is not a developer of SEXTANTE, but very active for QGIS. Therefore us does obviously not mean SEXTANTE but QGIS. Paolo is again appropriating python SEXTANTE (Please do nont send bug reports concerning it to other bugtrackers.) as the QGIS Python plugin called sextante. See [2] Yes, us means QGIS. We (QGIS) already decided to integrate SEXTANTE into QGIS master as soon as code freeze is over (we are now releasing QGIS 1.8), so the distinction is blurred. [7] http://osgeo-org.1560.n6.nabble.com/Report-from-ongoing-GRASS-GIS-Community-Sprint-in-Prague-tp4977136p4977500.html At this stage of the python version of SEXTANTE, the main objective should be to get a stable version of python SEXTANTE. That is, a generic python SEXTANTE, not just something that works for QGIS and nothing else (with the SEXTANTE/GRASS interface not even working properly, see main comment). Paolo is again appropriating python SEXTANTE (Please do nont send bug reports concerning it to other bugtrackers.) as the QGIS Python plugin called sextante. See [2] Where do you find a generic python SEXTANTE? Having said that, if this still hurts, I apologize once more: http://www.youtube.com/watch?v=m7mIy97_rlo All the best. -- Paolo Cavallini - Faunalia www.faunalia.eu Full contact details at www.faunalia.eu/pc ___ grass-user mailing list grass-user@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/grass-user
Re: [GRASS-user] Identifying Why Vector Map Lables Do Not Display
On 27/05/12 23:43, Rich Shepard wrote: I've a vector map of areas created from the raster original using r.to.vect with the '-v' option. The vector map's attribute table has cat values which is sufficient for me to identify specific polygons. However, when displaying the map the lables do not appear despite the attr and cat options being selected. You have to display the centroids to be able to see labels, i.e. d.vect census_wake2000 display=shape,cat type=point,line,boundary,area,face will not display the cat values, while d.vect census_wake2000 display=shape,cat type=point,line,boundary,area,face,centroid will. No need to use attr if all you want to see are the categories of the features (unless your attribute table's cat column does not correspond to the category values of the features, but represents a separate attribute, then attr is needed, but cat not.) Moritz ___ grass-user mailing list grass-user@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/grass-user
Re: [GRASS-user] v.dissolve renames the attribute column
On 28/05/12 12:30, Markus Metz wrote: On Mon, May 28, 2012 at 11:52 AM, spiderplant0spiderpla...@gmail.com wrote: Hi. It doesn't add any new information but I wanted the .dbf file to be there so I can open it up with other tools and check the contents for debug purposes etc. And I wanted to change the column to the correct name for consistency with downstream processing scripts. Try v.dissolve input=pass2@mapset1 layer=1 column=fttcp output=pass3 v.db.addtable map=pass3 v.db.renamecol map=pass3 column=cat,fftcp IIUC, if you do that you will lose the connection between the map and the attribute table if that connection is defined with the cat column as the key column. If you want to keep the connection, but use the name fftcp instead of cat for the key column, you will have to do something like this: v.dissolve input=pass2@mapset1 layer=1 column=fttcp output=pass3 v.db.addtable map=pass3 v.db.addcolumn map=pass3 col=fftcp int v.db.update map=pass3 col=fftcp qcolumn=cat v.db.connect -o map=pass3 table=pass3 key=fftcp v.db.dropcolumn map=pass3 col=cat If you want that to be done by v.dissolve, that is, v.reclass, automatically, you could file an enhancement ticket for v.reclass to add a new flag, say -t, to always create a table for the output where the original column name is used when the column type is integer. Maybe adding the option to v.db.addtable to use another name for the key column would also be good idea. Moritz ___ grass-user mailing list grass-user@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/grass-user
Re: [GRASS-user] v.dissolve renames the attribute column
Moritz Lennert wrote: On 28/05/12 12:30, Markus Metz wrote: On Mon, May 28, 2012 at 11:52 AM, spiderplant0spiderpla...@gmail.com wrote: Hi. It doesn't add any new information but I wanted the .dbf file to be there so I can open it up with other tools and check the contents for debug purposes etc. And I wanted to change the column to the correct name for consistency with downstream processing scripts. Try v.dissolve input=pass2@mapset1 layer=1 column=fttcp output=pass3 v.db.addtable map=pass3 v.db.renamecol map=pass3 column=cat,fftcp IIUC, if you do that you will lose the connection between the map and the attribute table if that connection is defined with the cat column as the key column. Oops, right, thanks for pointing this out. Markus If you want to keep the connection, but use the name fftcp instead of cat for the key column, you will have to do something like this: v.dissolve input=pass2@mapset1 layer=1 column=fttcp output=pass3 v.db.addtable map=pass3 v.db.addcolumn map=pass3 col=fftcp int v.db.update map=pass3 col=fftcp qcolumn=cat v.db.connect -o map=pass3 table=pass3 key=fftcp v.db.dropcolumn map=pass3 col=cat If you want that to be done by v.dissolve, that is, v.reclass, automatically, you could file an enhancement ticket for v.reclass to add a new flag, say -t, to always create a table for the output where the original column name is used when the column type is integer. Maybe adding the option to v.db.addtable to use another name for the key column would also be good idea. Moritz ___ grass-user mailing list grass-user@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/grass-user
Re: [GRASS-user] Identifying Why Vector Map Lables Do Not Display [FIXED]
On Tue, 29 May 2012, Moritz Lennert wrote: You have to display the centroids to be able to see labels, i.e. Moritz, You lead me to the solution. There is one default choice that must be maintained (the centroid checkbox on the Selections tab) and one non-default that needs to be selected: the the 'Display category number of features' on the Required tab. I had consistently left the centroids selected for display, but on the Required tab had checked the 'Display selected attribute from attrcol' box rather than the 'Display category number of features' box. Because the category numbers column appears on the Labels tab as a choice for display I interpreted it as an attribute column and missed the second-from-the-top choice on the Required tab. Thank you very much, Rich ___ grass-user mailing list grass-user@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/grass-user
Re: [GRASS-user] wxgui, wxpython Pronunciation
On Mon, 28 May 2012, Hamish wrote: the wx of wxWidgets (nee wxWindows) apparently comes from W for MS Windows, and X for UNIX's X-Windows: http://www.wxwidgets.org/about/history.htm w-x-widgets ? Hamish, Everyone I know pronouces them as you show aboveabove: w-x-python, w-x-gui, w-x-windows. Rich ___ grass-user mailing list grass-user@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/grass-user
Re: [GRASS-user] Identifying Why Vector Map Lables Do Not Display [FIXED]
On 29/05/12 14:58, Rich Shepard wrote: On Tue, 29 May 2012, Moritz Lennert wrote: You have to display the centroids to be able to see labels, i.e. Moritz, You lead me to the solution. There is one default choice that must be maintained (the centroid checkbox on the Selections tab) and one non-default that needs to be selected: the the 'Display category number of features' on the Required tab. I had consistently left the centroids selected for display, but on the Required tab had checked the 'Display selected attribute from attrcol' box rather than the 'Display category number of features' box. Because the category numbers column appears on the Labels tab as a choice for display I interpreted it as an attribute column and missed the second-from-the-top choice on the Required tab. Well, actually, it should work both ways, i.e. with the d.vect census_wake2000 display=shape,attr type=point,line,boundary,area,face,centroid attrcol=cat d.vect census_wake2000 display=shape,cat type=point,line,boundary,area,face,centroid should give the same result (assuming that all features that have a category value have a corresponding entry in the cat column of the attribute table). Note that in grass7.0 the use of centroids is not necessary anymore. Moritz ___ grass-user mailing list grass-user@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/grass-user
