[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