[GRASS-SVN] r58637 - in grass-addons/grass7/vector: . v.stats

svn_grass at osgeo.org svn_grass at osgeo.org
Tue Jan 7 08:18:28 PST 2014


Author: zarch
Date: 2014-01-07 08:18:28 -0800 (Tue, 07 Jan 2014)
New Revision: 58637

Added:
   grass-addons/grass7/vector/v.stats/
   grass-addons/grass7/vector/v.stats/Makefile
   grass-addons/grass7/vector/v.stats/imp_csv.py
   grass-addons/grass7/vector/v.stats/v.stats.html
   grass-addons/grass7/vector/v.stats/v.stats.py
   grass-addons/grass7/vector/v.stats/vstats.py
Log:
Add a new module to combine shape and raster statistics

Added: grass-addons/grass7/vector/v.stats/Makefile
===================================================================
--- grass-addons/grass7/vector/v.stats/Makefile	                        (rev 0)
+++ grass-addons/grass7/vector/v.stats/Makefile	2014-01-07 16:18:28 UTC (rev 58637)
@@ -0,0 +1,10 @@
+MODULE_TOPDIR = ../..
+
+PGM = v.stats
+
+ETCFILES = vstats imp_csv
+
+include $(MODULE_TOPDIR)/include/Make/Script.make
+include $(MODULE_TOPDIR)/include/Make/Python.make
+
+default: script

Added: grass-addons/grass7/vector/v.stats/imp_csv.py
===================================================================
--- grass-addons/grass7/vector/v.stats/imp_csv.py	                        (rev 0)
+++ grass-addons/grass7/vector/v.stats/imp_csv.py	2014-01-07 16:18:28 UTC (rev 58637)
@@ -0,0 +1,164 @@
+# -*- coding: utf-8 -*-
+from __future__ import print_function, division
+import numpy as np
+
+import grass.lib.gis as glg
+
+NPY2COLTYPE = {'<i8': 'INTEGER',
+               '<f8': 'DOUBLE'}
+
+# shp column names
+ALLSHPN = ['area_id',    'cat',   'nisles',     'x_extent', 'y_extent',
+           'iperimeter', 'iarea', 'icompact',   'ifd',      'perimeter',
+           'area',       'boundarea', 'aratio', 'compact',  'fd']
+
+# shp column types
+ALLSHPT = ['<i8', '<i8', '<i8', '<f8', '<f8',
+           '<f8', '<f8', '<f8', '<f8', '<f8',
+           '<f8', '<f8', '<f8', '<f8', '<f8']
+
+# rst column names
+ALLRSTN = ['zone', 'label', 'all_cells', 'non_null_cells', 'null_cells',
+           'sum', 'sum_abs', 'min', 'max', 'range',
+           'mean', 'mean_of_abs', 'variance', 'stddev', 'coeff_var',
+           'skewness', 'kurtosis', 'variance2', 'stddev2', 'coeff_var2',
+           'skewness2', 'kurtosis2', 'first_quart', 'median', 'third_quart',
+           'perc_90', 'mode', 'occurrences']
+
+# rst column type
+ALLRSTT = ['<i8', '<i8', '<i8', '<i8', '<i8',
+           '<f8', '<f8', '<f8', '<f8', '<f8',
+           '<f8', '<f8', '<f8', '<f8', '<f8',
+           '<f8', '<f8', '<f8', '<f8', '<f8',
+           '<f8', '<f8', '<f8', '<f8', '<f8',
+           '<f8', '<f8', '<f8']
+
+SKIPCSV = ['label', 'non_null_cells', 'null_cells',
+           'mean_of_abs', 'sum', 'sum_abs']
+
+SKIPSHP = ['area_id', ]
+
+#SKIPSHP = ['areas', 'red__n', 'red__min', 'red__max', 'red__range', 'red__sum',
+#           'red__mean', 'red__stddev', 'red__variance', 'red__cf_var',
+#           'red__first_quartile', 'red__median', 'red__third_quartile',
+#           'red__percentile_90', 'green__n', 'green__min', 'green__max',
+#           'green__range', 'green__mean', 'green__stddev', 'green__variance',
+#           'green__cf_var', 'green__sum', 'green__first_quartile',
+#           'green__median', 'green__third_quartile', 'green__percentile_90',
+#           'blue__n', 'blue__min', 'blue__max', 'blue__range', 'blue__mean',
+#           'blue__stddev', 'blue__variance', 'blue__cf_var', 'blue__sum',
+#           'blue__first_quartile', 'blue__median', 'blue__third_quartile',
+#           'blue__percentile_90']
+
+#----------------------------------------
+
+
+def getrow(shpdata, rstdata, dim):
+    drow = np.zeros((dim, ))
+    lenght = len(shpdata)
+    for i, rows in enumerate(zip(shpdata, *[rst[:, 1:] for rst in rstdata])):
+        glg.G_percent(i, lenght, 2)
+        start = 0
+        for row in rows:
+            end = start + row.shape[0]
+            drow[start:end] = row
+            start = end
+        yield drow
+
+
+def getcols(names, types, skipnames=None, prefix='',
+            gettype=lambda x: x):
+    cols = []
+    skipnames = skipnames if skipnames else []
+    for cname, ctype in zip(names, types):
+        if cname not in skipnames:
+            cols.append((prefix + cname, gettype(ctype)))
+    return cols
+
+
+def gettablecols(prefixes, allshpn, allshpt, skipshp,
+                 allrstn, allrstt, skiprst):
+    gettype = lambda x: NPY2COLTYPE[x]
+    # define the new columns
+    cols = getcols(allshpn, allshpt, skipshp, gettype=gettype)
+    for prfx in prefixes:
+        cols += getcols(allrstn, allrstt, skiprst,
+                        prefix=prfx, gettype=gettype)
+    return cols
+
+
+def update_cols(tab, shpcsv, rstcsv, prefixes=None,
+                skipshp=None, skiprst=None,
+                shpcat=0, rstcat=0,
+                allshpn=ALLSHPN, allrstn=ALLRSTN,
+                allshpt=ALLSHPT, allrstt=ALLRSTT, overwrite=False):
+    prefixes = prefixes if prefixes else [csv[:-4] for csv in rstcsv]
+    skipshp = skipshp if skipshp else []
+    skiprst = skiprst if skiprst else []
+    useshpcols = [allshpn.index(c) for c in allshpn if c not in skipshp]
+    userstcols = [allrstn.index(c) for c in allrstn if c not in skiprst]
+
+    print("Start loading data from:")
+    print("    - %s." % shpcsv, end='')
+    shpdata = np.genfromtxt(shpcsv, delimiter=';', names=True,
+                            usecols=useshpcols,
+                            dtype=getcols(allshpn, allshpt, skipshp))
+    shpdata.sort(order='cat')
+    print(' Done.')
+
+    rstdata = []
+    for rst in rstcsv:
+        print("    - %s." % rst, end='')
+        rstdata.append(np.genfromtxt(rst, delimiter='|', names=True,
+                                     usecols=userstcols,
+                                     missing_values=('nan', '-nan'),
+                                     dtype=getcols(allrstn, allrstt,
+                                                   userstcols)))
+        rstdata[-1].sort(order='zone')
+        print(' Done.')
+
+    print("Cheking categories and zones correspondance:")
+    for i, rst in enumerate(rstdata):
+        print("    - <%s>." % rstcsv[i], end='')
+#        if rstcsv[i] == 'median5_contr.csv':
+#            import ipdb; ipdb.set_trace()
+        if not (shpdata['cat'] == rst['zone']).all():
+            msg = "The categories and the zones are not equal, in <%s>"
+            raise ValueError(msg % rstcsv[i])
+        print(' Ok.')
+
+    print("Conversion from record array to array")
+    shpdata = np.array(shpdata.tolist())
+    print('.', end='')
+    for i, rst in enumerate(rstdata):
+        rstdata[i] = np.array(rst.tolist())
+        print('.', end='')
+    print('Done.')
+    # create the new table
+    if tab.exist():
+        tab.drop(force=True)
+    cols = gettablecols(prefixes, allshpn, allshpt, skipshp,
+                        allrstn, allrstt, skiprst)
+    tab.create(cols)
+    cur = tab.conn.cursor()
+    print("Merge shape table with raster csv.")
+    tab.insert(getrow(shpdata, rstdata, len(cols)), cursor=cur, many=True)
+    tab.conn.commit()
+    print("%d rows inserted." % tab.n_rows())
+    return tab
+
+
+#link = None
+#with VectorTopo(VSEG, mode='r') as vect:
+#    #link = update_cols(vect.table, CSV, PREFIX, allcsvcols=ALL,
+#    #                   skipcsv=SKIPCSV, skipshp=SKIPSHP)
+#    link = update_cols(vect.table, ECSV, PREFIX, allcsvcols=ALL,
+#                       skipcsv=SKIPCSV, skipshp=ESKIPSHP)
+#    link.layer = vect.layer + 1
+#
+##----------------------------
+#with Vector(VSEG, mode='rw') as vect:
+#    link = Link(layer=NEW_LAYER, name=NEW_NAME_LAYER, table=NEW_TABLE_NAME)
+#    vect.dblinks.add(link)
+
+

Added: grass-addons/grass7/vector/v.stats/v.stats.html
===================================================================
Added: grass-addons/grass7/vector/v.stats/v.stats.py
===================================================================
--- grass-addons/grass7/vector/v.stats/v.stats.py	                        (rev 0)
+++ grass-addons/grass7/vector/v.stats/v.stats.py	2014-01-07 16:18:28 UTC (rev 58637)
@@ -0,0 +1,206 @@
+#!/usr/bin/env python
+# -- coding: utf-8 --
+#
+############################################################################
+#
+# MODULE:	    i.segment.hierarchical
+#
+# AUTHOR(S):   Pietro Zambelli (University of Trento)
+#
+# COPYRIGHT:	(C) 2013 by 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: Vector stats
+#%  keywords: vector
+#%  keywords: statistics
+#%  keywords: shape
+#%  overwrite: yes
+#%End
+#%option G_OPT_V_MAP
+#%  key: vector
+#%  description: Name of input vector map
+#%  required: yes
+#%end
+#%option
+#%  key: layer
+#%  type: integer
+#%  description: Vector layer
+#%  multiple: no
+#%  required: no
+#%  answer: 1
+#%end
+#%option G_OPT_R_INPUTS
+#%  key: rasters
+#%  description: Name of input raster maps
+#%  multiple: yes
+#%  required: no
+#%end
+#%option G_OPT_R_INPUT
+#%  key: zones
+#%  description: Name of raster zones map
+#%  multiple: no
+#%  required: no
+#%end
+#%option
+#%  key: rprefix
+#%  description: Raster prefixes
+#%  multiple: yes
+#%  required: no
+#%end
+#%option
+#%  key: skipshape
+#%  type: string
+#%  multiple: yes
+#%  description: Skip shape columns
+#%  required: no
+#%  answer: area_id
+#%end
+#%option
+#%  key: skipunivar
+#%  type: string
+#%  multiple: yes
+#%  description: Skip shape columns
+#%  required: no
+#%  answer: label,non_null_cells,null_cells,mean_of_abs,sum,sum_abs
+#%end
+#%option
+#%  key: shpcsv
+#%  type: string
+#%  multiple: no
+#%  description: CSV with the vector statistics of the shape
+#%  required: no
+#%end
+#%option
+#%  key: rstcsv
+#%  type: string
+#%  multiple: yes
+#%  description: CSV with the statistics of the raster maps
+#%  required: no
+#%end
+#%option
+#%  key: rstpercentile
+#%  type: integer
+#%  multiple: no
+#%  description: Raster percentile to use
+#%  required: no
+#%  answer: 90
+#%end
+#%option
+#%  key: newlayer
+#%  type: integer
+#%  description: New vector layer that will be add to the vector map
+#%  multiple: no
+#%  required: no
+#%  answer: 2
+#%end
+#%option
+#%  key: newlayername
+#%  type: string
+#%  description: New vector layer that will be add to the vector map
+#%  multiple: no
+#%  required: no
+#%end
+#%option
+#%  key: newtabname
+#%  type: string
+#%  description: New vector layer that will be add to the vector map
+#%  multiple: no
+#%  required: no
+#%end
+#%flag
+#%  key: r
+#%  description: Read from existing CSV files
+#%end
+
+#-----------------------------------------------------
+"""
+# convert segments to vector
+r.to.vect input=seg_0.05 at pietro output=seg005 type=area
+v.category input=seg005 layer=1,2,3,4,5 output=seg_005 type=area option=transfer
+v.to.rast input=seg_005 output=vseg_005
+
+"""
+import sys
+import os
+
+from grass.script.core import parser
+
+
+from grass.pygrass.functions import get_lib_path, get_mapset_raster
+from grass.pygrass.vector import VectorTopo, Vector
+from grass.pygrass.vector.table import Link
+
+path = get_lib_path("v.stats", "")
+if path is None:
+    raise ImportError("Not able to find the path %s directory." % path)
+
+sys.path.append(path)
+
+from vstats import get_shp_csv, get_zones, get_rst_csv
+from imp_csv import update_cols
+
+
+def main(opt, flg):
+    #
+    # Set check variables
+    #
+    overwrite = True
+    rasters = opt['rasters'].split(',') if opt['rasters'] else []
+    rprefix = opt['rprefix'].split(',') if opt['rprefix'] else []
+
+    split = lambda x: x.split('@') if '@' in x else (x, '')
+
+    vname, vmset = split(opt['vector'])
+    shpcsv = opt['shpcsv'] if opt['shpcsv'] else vname + '.csv'
+    rstcsv = (opt['rstcsv'].split(',') if opt['rstcsv']
+              else [split(rst)[0] + '.csv' for rst in rasters])
+    zones = opt['zones'] if opt['zones'] else vname + '_zones'
+    if rasters:
+        if rprefix and len(rasters) != len(rprefix):
+            raise
+        if len(rasters) != len(rstcsv):
+            raise
+        prefixes = rprefix if rprefix else rasters
+
+    skipshp = opt['skipshape'].split(',') if opt['skipshape'] else []
+    skiprst = opt['skipunivar'].split(',') if opt['skipunivar'] else []
+    layer = int(opt['layer'])
+    newlayer = int(opt['newlayer'])
+    newlayername = (opt['newlayername'] if opt['newlayername']
+                    else vname + '_stats')
+    newtabname = opt['newtabname'] if opt['newtabname'] else vname + '_stats'
+    rstpercentile = float(opt['rstpercentile'])
+
+    #
+    # compute
+    #
+    if not os.path.exists(shpcsv):
+        get_shp_csv(opt['vector'], shpcsv, overwrite)
+    if not get_mapset_raster(zones):
+        get_zones(opt['vector'], zones, layer)
+    if not rstcsv or not os.path.exists(rstcsv[0]):
+        get_rst_csv(rasters, zones, rstcsv, rstpercentile, overwrite)
+
+    newlink = Link(newlayer, newlayername, newtabname)
+    newtab = newlink.table()
+
+    with Vector(vname, vmset, mode='r', layer=layer) as vct:
+        mode = 'r' if newlink in vct.dblinks else 'rw'
+
+    with VectorTopo(vname, vmset, mode=mode, layer=layer) as vct:
+        update_cols(newtab, shpcsv, rstcsv, prefixes, skipshp, skiprst)
+
+        if mode == 'rw':
+            # add the new link
+            vct.dblinks.add(newlink)
+            vct.build()
+
+
+if __name__ == "__main__":
+    main(*parser())


Property changes on: grass-addons/grass7/vector/v.stats/v.stats.py
___________________________________________________________________
Added: svn:executable
   + *

Added: grass-addons/grass7/vector/v.stats/vstats.py
===================================================================
--- grass-addons/grass7/vector/v.stats/vstats.py	                        (rev 0)
+++ grass-addons/grass7/vector/v.stats/vstats.py	2014-01-07 16:18:28 UTC (rev 58637)
@@ -0,0 +1,36 @@
+# -*- coding: utf-8 -*-
+"""
+Created on Mon Oct 21 12:54:42 2013
+
+ at author: pietro
+"""
+from __future__ import print_function, division
+import os
+from grass.pygrass.modules import Module
+
+
+def get_shp_csv(vector, csv=None, overwrite=False):
+    vasts = Module('v.area.stats')
+    csv = vector + '.csv' if csv is None else csv
+    if os.path.exists(csv) and overwrite:
+        os.remove(csv)
+    vasts(map=vector, output=csv, overwrite=overwrite)
+    return csv
+
+
+def get_zones(vector, zones, layer=1, overwrite=False):
+    v2rast = Module('v.to.rast', input=vector, layer=str(layer), type='area',
+                    output=zones, overwrite=overwrite, rows=65536, use='cat')
+    rclr = Module("r.colors", map=zones, color="random")
+
+
+def get_rst_csv(rasters, zones, csvfiles, percentile=90., overwrite=False):
+    procs = []
+    for rast, csv in zip(rasters, csvfiles):
+        procs.append(Module('r.univar2', map=rast, zones=zones,
+                            percentile=percentile, output=csv,
+                            overwrite=overwrite, flags='e', finish_=False))
+    # wait the end of all process
+    for proc in procs:
+        proc.popen.wait()
+    return csvfiles



More information about the grass-commit mailing list