[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