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
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()
_______________________________________________ grass-dev mailing list grass-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/grass-dev