[GRASS-dev] port of v.points.cog to python
Patrice Dumas
pertusus at free.fr
Sun Mar 1 08:01:15 PST 2015
Hello,
Here is a rewrite of the v.points.cog module in python. I tried
translating code only, keeping the code organization, variables names as it
was previously to help those who would want to review the differences
with the shell script. I also did the few changes required for changes
in grass7, but there weren't that much since python function already
abstract some commands.
I did a svn copy of the grass6 module before doing the modifications.
I cannot (and don't want to...) claim any copyright on a code translation.
I didn't test thoroughly, but since it is code translation, I don't expect
big surprises.
--
Pat
-------------- next part --------------
Index: grass7/vector/Makefile
===================================================================
--- grass7/vector/Makefile (r?vision 64773)
+++ grass7/vector/Makefile (copie de travail)
@@ -35,6 +35,7 @@
v.out.ply \
v.out.png \
v.ply.rectify \
+ v.points.cog \
v.stats \
v.surf.icw \
v.surf.mass \
Index: grass7/vector/v.points.cog/description.html
===================================================================
--- grass7/vector/v.points.cog/description.html (r?vision 64773)
+++ grass7/vector/v.points.cog/description.html (copie de travail)
@@ -1,53 +0,0 @@
-<h2>DESCRIPTION</h2>
-
-
-<em>v.points.cog</em> condenses points or centroids sharing
-a common attribute into a single point in a new vector map at their
-average position (center of gravity).
-<p>
-For this to work well your clusters of points must be gregarious
-(Gaussian distribution) - if two groups of points habitate in two
-corners of the map the output point will fall in the center and
-match niether well.
-<p>
-If needed you can use <em>v.digit</em> to adjust point positions
-created by this module, <!-- repel 2 points if v.distance < minimum? -->
-or use this module as a preprocessing step before running
-<em>v.label.sa</em> to deal with overlapping labels automatically.
-
-<!--<p> appropriate for lat/lon? (weighted avg for y by cos(y)?) -->
-
-
-<h2>EXAMPLE</h2>
-
-Create single points at the average of all points in the
-<tt>bugsites</tt> map, and place a single label there.
-
-<div class="code"><pre>
-v.points.cog in=bugsites out=bug_cog column=str1
-
-d.vect -c bugsites color=none icon=basic/circle
-
-d.vect bug_cog disp=attr attrcol=str1 lcolor=black \
- lsize=14 xref=center yref=center bgcolor=white
-</pre></div>
-
-
-<h2>SEE ALSO</h2>
-<em>
-<a href="v.label.sa.html">v.label.sa</a><br>
-<a href="v.label.html">v.label</a><br>
-<a href="v.digit.html">v.digit</a>
-</em>
-
-
-<h2>AUTHOR</h2>
-
-Hamish Bowman<br>
-<i>Dept. Marine Science<br>
-University of Otago<br>
-Dunedin, New Zealand</i>
-<br>
-
-<p>
-<i>Last changed: $Date$</i>
Index: grass7/vector/v.points.cog/v.points.cog
===================================================================
--- grass7/vector/v.points.cog/v.points.cog (r?vision 64773)
+++ grass7/vector/v.points.cog/v.points.cog (copie de travail)
@@ -1,142 +0,0 @@
-#!/bin/sh
-#
-############################################################################
-#
-# MODULE: v.points.cog
-#
-# AUTHOR(S): Hamish Bowman
-#
-# PURPOSE: Condense points or centroids sharing a common attribute into a single point
-#
-# COPYRIGHT: (c) 2010 Hamish Bowman, and the GRASS Development Team
-#
-# This program is free software under the GNU General Public
-# License (>=v2). Read the file COPYING that comes with GRASS
-# for details.
-#
-#############################################################################
-#%Module
-#% description: Condense points or centroids sharing a common attribute into a single point.
-#% keywords: vector, cluster
-#%End
-#%option
-#% key: input
-#% type: string
-#% gisprompt: old,vector,vector
-#% description: Name of input vector map
-#% required : yes
-#%end
-#%option
-#% key: output
-#% type: string
-#% gisprompt: new,vector,vector
-#% description: Name for output vector map
-#% required : yes
-#%end
-#%option
-#% key: column
-#% type: string
-#% description: Column containing common attribute
-#% required : yes
-#%end
-#%option
-#% key: layer
-#% type: integer
-#% answer: 1
-#% description: Layer number
-#% required: no
-#%end
-
-##%option
-##% key: type
-##% type: string
-##% description: Feature type(s)
-##% options: point,centroid
-##% answer: point
-##% required: no
-##% multiple: yes
-##%end
-
-
-if [ -z "$GISBASE" ] ; then
- echo "You must be in GRASS GIS to run this program." 1>&2
- exit 1
-fi
-
-if [ "$1" != "@ARGS_PARSED@" ] ; then
- exec g.parser "$0" "$@"
-fi
-
-MAP="$GIS_OPT_INPUT"
-COLUMN="$GIS_OPT_COLUMN"
-LAYER="$GIS_OPT_LAYER"
-
-# check for input map
-eval `g.findfile element=vector file="$MAP"`
-if [ ! "$file" ] ; then
- g.message -e "Vector map <$MAP> does not exist."
- exit 1
-fi
-
-# check for column
-if [ `v.info -c "$MAP" layer="$LAYER" --quiet | cut -f2 -d'|' | grep -c "^$COLUMN$"` -ne 1 ] ; then
- g.message -e "Column <$COLUMN> not found."
- exit 1
-fi
-
-# get column details so we can recreate it.
-# The db.* modules need special care when querying from @another mapset
-DB=`v.db.connect "$MAP" -g layer="$LAYER" fs='|' | cut -f4 -d'|'`
-BASENAME=`echo "$MAP" | sed -e 's/@.*//'`
-db.describe -c table="$BASENAME" database="$DB" > /dev/null
-if [ $? -ne 0 ] ; then
- g.message -e "Unable to describe table"
- exit 1
-fi
-COLUMN_DESC=`db.describe -c table="$BASENAME" database="$DB" | grep " $COLUMN:" | cut -f3- -d:`
-
-
-if [ `echo "$COLUMN_DESC" | grep -c CHARACTER` -eq 1 ] ; then
- COLUMN_TYPE="string"
- COLUMN_LEN=`echo "$COLUMN_DESC" | cut -f2 -d:`
- COLUMN_DEFN="varchar($COLUMN_LEN)"
-else
- COLUMN_TYPE="number"
- COLUMN_DEFN=`echo "$COLUMN_DESC" | cut -f1 -d:`
-fi
-
-# cheap hack to avoid conflict
-if [ "$COLUMN" = "cat" ] ; then
- OUT_COLUMN="cat_"
-else
- OUT_COLUMN="$COLUMN"
-fi
-
-
-(
-IFS='|'
-for ITEM in `v.db.select "$MAP" -c column="$COLUMN" layer="$LAYER" | sort | uniq | tr '\n' '|'` ; do
- #echo "[$ITEM]"
- if [ "$COLUMN_TYPE" = "string" ] ; then
- WHERE_STR="$COLUMN = '$ITEM'"
- else
- WHERE_STR="$COLUMN = $ITEM"
- fi
-
- v.out.ascii "$MAP" column="$COLUMN" layer="$LAYER" where="$WHERE_STR" | \
- awk -F'|' \
- 'BEGIN { sum_x=0; sum_y=0; i=0 }
- { sum_x += $1; sum_y += $2; i++ }
- END { if(i>0) { printf("%.15g|%.15g|%s\n", sum_x/i, sum_y/i, $4) } }'
-done
-) | v.in.ascii out="$GIS_OPT_OUTPUT" columns="x double, y double, $OUT_COLUMN $COLUMN_DEFN"
-#echo columns="x double, y double, $OUT_COLUMN $COLUMN_DEFN"
-
-retval=$?
-
-# cleanup cheap hack
-if [ $COLUMN = "cat" ] ; then
- v.db.dropcol map="$GIS_OPT_OUTPUT" column="cat_" --quiet
-fi
-
-exit $retval
Index: grass7/vector/v.points.cog/v.points.cog.py
===================================================================
--- grass7/vector/v.points.cog/v.points.cog.py (r?vision 64773)
+++ grass7/vector/v.points.cog/v.points.cog.py (copie de travail)
@@ -1,4 +1,4 @@
-#!/bin/sh
+#!/usr/bin/env python
#
############################################################################
#
@@ -57,86 +57,74 @@
##% multiple: yes
##%end
+import sys
+import grass.script as grass
-if [ -z "$GISBASE" ] ; then
- echo "You must be in GRASS GIS to run this program." 1>&2
- exit 1
-fi
+def main():
-if [ "$1" != "@ARGS_PARSED@" ] ; then
- exec g.parser "$0" "$@"
-fi
+ map = options['input']
+ column = options['column']
+ layer = options['layer']
+ output = options['output']
-MAP="$GIS_OPT_INPUT"
-COLUMN="$GIS_OPT_COLUMN"
-LAYER="$GIS_OPT_LAYER"
+ result = grass.find_file(map, element='vector')
+ if len(result['name']) == 0:
+ grass.fatal(_("Vector map <%s> does not exist") % map)
-# check for input map
-eval `g.findfile element=vector file="$MAP"`
-if [ ! "$file" ] ; then
- g.message -e "Vector map <$MAP> does not exist."
- exit 1
-fi
+ if column not in grass.vector_columns(map, layer).keys():
+ grass.fatal(_("Column <%s>, layer %s not found") % (map, layer))
-# check for column
-if [ `v.info -c "$MAP" layer="$LAYER" --quiet | cut -f2 -d'|' | grep -c "^$COLUMN$"` -ne 1 ] ; then
- g.message -e "Column <$COLUMN> not found."
- exit 1
-fi
+ # get column details so we can recreate it.
+ # The db.* modules need special care when querying from @another mapset
+ map_connection_info = grass.vector_db(map)[int(layer)]
+ db = map_connection_info['database']
+ basename = map.split("@")[0]
-# get column details so we can recreate it.
-# The db.* modules need special care when querying from @another mapset
-DB=`v.db.connect "$MAP" -g layer="$LAYER" fs='|' | cut -f4 -d'|'`
-BASENAME=`echo "$MAP" | sed -e 's/@.*//'`
-db.describe -c table="$BASENAME" database="$DB" > /dev/null
-if [ $? -ne 0 ] ; then
- g.message -e "Unable to describe table"
- exit 1
-fi
-COLUMN_DESC=`db.describe -c table="$BASENAME" database="$DB" | grep " $COLUMN:" | cut -f3- -d:`
+ map_description = grass.db_describe(basename, database=db)
+ if 'cols' not in map_description:
+ grass.fatal(_("Unable to describe table"))
+ for column_desc in map_description['cols']:
+ if column_desc[0] == column:
+ break
+ if 'CHARACTER' in column_desc[1]:
+ column_type = "string"
+ column_len = column_desc[2]
+ column_defn = "varchar(%d)" % (column_len)
+ else:
+ column_type = "number"
+ column_defn = column_desc[1]
-if [ `echo "$COLUMN_DESC" | grep -c CHARACTER` -eq 1 ] ; then
- COLUMN_TYPE="string"
- COLUMN_LEN=`echo "$COLUMN_DESC" | cut -f2 -d:`
- COLUMN_DEFN="varchar($COLUMN_LEN)"
-else
- COLUMN_TYPE="number"
- COLUMN_DEFN=`echo "$COLUMN_DESC" | cut -f1 -d:`
-fi
+ # cheap hack to avoid conflict
+ if column == 'cat':
+ out_column = "cat_"
+ else:
+ out_column = column
-# cheap hack to avoid conflict
-if [ "$COLUMN" = "cat" ] ; then
- OUT_COLUMN="cat_"
-else
- OUT_COLUMN="$COLUMN"
-fi
+ average_points_positions = ''
+ for item in sorted(set(grass.vector_db_select(map, columns=column, layer=int(layer))['values'])):
+ if column_type == "string":
+ where_str = "%s = '%s'" % (column, item)
+ else:
+ where_str = "%s = %s" % (column, str(item))
+
+ out_ascii_output = grass.read_command('v.out.ascii', input=map, output='-', column=column, layer=layer, where=where_str).splitlines()
+ if len(out_ascii_output) > 1:
+ sum_x = 0.
+ sum_y = 0.
+ for item_point in out_ascii_output:
+ position = item_point.split('|')
+ sum_x += float(position[0])
+ sum_y += float(position[1])
+ average_points_positions += "%.15g|%.15g|%s\n" % (sum_x/len(out_ascii_output), sum_y/len(out_ascii_output), position[-1])
+ retval = grass.write_command('v.in.ascii', stdin=average_points_positions, input='-', output=output, columns="x double, y double, %s %s" % (out_column, column_defn))
+ # cleanup cheap hack
+ if column == 'cat':
+ grass.run_command('v.db.dropcol', map=output, column='cat_', quiet=True)
-(
-IFS='|'
-for ITEM in `v.db.select "$MAP" -c column="$COLUMN" layer="$LAYER" | sort | uniq | tr '\n' '|'` ; do
- #echo "[$ITEM]"
- if [ "$COLUMN_TYPE" = "string" ] ; then
- WHERE_STR="$COLUMN = '$ITEM'"
- else
- WHERE_STR="$COLUMN = $ITEM"
- fi
+ sys.exit(retval)
- v.out.ascii "$MAP" column="$COLUMN" layer="$LAYER" where="$WHERE_STR" | \
- awk -F'|' \
- 'BEGIN { sum_x=0; sum_y=0; i=0 }
- { sum_x += $1; sum_y += $2; i++ }
- END { if(i>0) { printf("%.15g|%.15g|%s\n", sum_x/i, sum_y/i, $4) } }'
-done
-) | v.in.ascii out="$GIS_OPT_OUTPUT" columns="x double, y double, $OUT_COLUMN $COLUMN_DEFN"
-#echo columns="x double, y double, $OUT_COLUMN $COLUMN_DEFN"
-
-retval=$?
-
-# cleanup cheap hack
-if [ $COLUMN = "cat" ] ; then
- v.db.dropcol map="$GIS_OPT_OUTPUT" column="cat_" --quiet
-fi
-
-exit $retval
+if __name__ == "__main__":
+ options, flags = grass.parser()
+ main()
More information about the grass-dev
mailing list