#!/usr/bin/env python #**************************************************************** #* #* MODULE: tin.sterro_riporto, v 1.0.0 #* #* AUTHOR(S): Antonio Alliegro #* #* PURPOSE: Calcolo sterro e riporto di un'area mediante #* campionamento #* #* COPYRIGHT: (C) 2011 Antonio Alliegro Civil Engineer #* Salerno, Italy #* antonioall(at)libero.it #* #* First Version: 2011/09/28 #* Last Version: 2011/09/28 #* #* This program is free software under the #* GNU General Public License (>=v2). #* Read the file COPYING that comes with GRASS #* for details. #* #**************************************************************** #%module #% description: Calcolo sterro e riporto di un'area mediante campionamento #% keywords: TIN, vector #%end #%option #% key: map #% type: string #% gisprompt: old,vector,vector #% key_desc: vector #% description: Copertura area di calcolo #% required: yes #%end #%option #% key: sf #% type: string #% gisprompt: old,vector,vector #% key_desc: tin #% description: Copertura TIN stato di fatto #% required: yes #%end #%option #% key: sp #% type: string #% gisprompt: old,vector,vector #% key_desc: tin #% description: Copertura TIN di progetto #% required: yes #%end #%option #% key: delta #% type: double #% answer: 1.0 #% key_desc: float #% description: Delta campionamento #% required: yes #%end import os, sys import grass.script as grass from grass.lib.gis import * from grass.lib.vector import * def main(): area = options['map'] statof = options['sf'] statop = options['sp'] delta = float(options['delta']) sterro = 0.0 riporto = 0.0 h = 0.0 ns = 0 nr = 0 # Inializza la libreria GRASS G_gisinit('') # Controlla se le coperture sono nel percorso di ricerca mapset = G_find_vector2(area, "") if not mapset: grass.fatal("Il file vettoriale <%s> non e' stato trovato" % area) mapset_f = G_find_vector2(statof, "") if not mapset_f: grass.fatal("Il file vettoriale <%s> non e' stato trovato" % statof) mapset_p = G_find_vector2(statop, "") if not mapset_p: grass.fatal("Il file vettoriale <%s> non e' stato trovato" % statop) # Definisce la struttura delle coperture map_info = pointer(Map_info()) map_info_f = pointer(Map_info()) map_info_p = pointer(Map_info()) # Setta le coperture al livello 2 (level 2: topology) Vect_set_open_level(2) # Apre i files vettoriali Vect_open_old(map_info, area, mapset) Vect_open_old(map_info_f, statof, mapset_f) Vect_open_old(map_info_p, statop, mapset_p) # Verifica se il TIN e' un vettoriale 3D if Vect_is_3d(map_info_f): grass.message("Il file vettoriale <%s> e' 3D" % statof) else: Vect_close(map_info) Vect_close(map_info_f) Vect_close(map_info_p) grass.fatal("Il file vettoriale <%s> non e' 3D" % statof) if Vect_is_3d(map_info_p): grass.message("Il file vettoriale <%s> e' 3D" % statop) else: Vect_close(map_info) Vect_close(map_info_f) Vect_close(map_info_p) grass.fatal("Il file vettoriale <%s> non e' 3D" % statop) # Verifica se le coperture hanno tipologia areale nareas = Vect_get_num_areas(map_info) if nareas == 0: Vect_close(map_info) Vect_close(map_info_f) Vect_close(map_info_p) grass.fatal("Il file vettoriale <%s> non ha tipologia areale" % area) nareas = Vect_get_num_areas(map_info_f) if nareas == 0: Vect_close(map_info) Vect_close(map_info_f) Vect_close(map_info_p) grass.fatal("Il file vettoriale <%s> non e' un TIN" % statof) nareas = Vect_get_num_areas(map_info_p) if nareas == 0: Vect_close(map_info) Vect_close(map_info_f) Vect_close(map_info_p) grass.fatal("Il file vettoriale <%s> non e' un TIN" % statop) # Estrazione valori della regione corrente window = pointer(Cell_head()) G_get_window(window) nrows = window.contents.rows ncols = window.contents.cols xref = window.contents.west yref = window.contents.south xres = window.contents.ew_res yres = window.contents.ns_res # Calcolo numero righe e colonne nuova griglia nx = int(nrows * xres / delta) + 1 ny = int(ncols * yres / delta) + 1 # Blocco elaborazione sterro/riporto grass.message("Passo 1/3: Pre-elaborazione TIN stato di fatto...") zf = c_double() zp = c_double() G_percent(0, nx, 2) Vect_tin_get_z(map_info_f, xref, yref, byref(zf), None, None) grass.message("...Fatto") grass.message("Passo 2/3: Pre-elaborazione TIN di progetto...") G_percent(0, nx, 2) Vect_tin_get_z(map_info_p, xref, yref, byref(zp), None, None) grass.message("...Fatto") grass.message("Passo 3/3: Calcolo volumi di sterro e riporto...") for i in range(nx): for j in range(ny): x = xref + j * delta y = yref + i * delta if Vect_point_in_area(map_info, 1, x, y) == 1: Vect_tin_get_z(map_info_f, x, y, byref(zf), None, None) Vect_tin_get_z(map_info_p, x, y, byref(zp), None, None) h = zp.value - zf.value if h < 0.0: sterro += abs(h) ns += 1 elif h > 0.0: riporto += h nr += 1 G_percent(i, nx, 2) grass.message("...Fatto") # Risultati volumes = sterro * delta * delta volumer = riporto * delta * delta volumec = volumes - volumer grass.message("-------------------------------------------") grass.message("Delta campionamento =%11.3f" % delta) grass.message("Numero celle totali =%11d" % (nx * ny)) grass.message("Numero celle sterro =%11d" % ns) grass.message("Numero celle riporto =%11d" % nr) grass.message("Volume sterro =%15.3f" % volumes) grass.message("Volume riporto =%15.3f" % volumer) if volumec < 0.0: grass.message("Volume eccedente =%15.3f riporto" % abs(volumec)) elif volumec == 0.0: grass.message("Volume eccedente =%15.3f compenso" % volumec) else: grass.message("Volume eccedente =%15.3f sterro" % volumec) # Chiusura coperture Vect_close(map_info) Vect_close(map_info_f) Vect_close(map_info_p) return 0 #End main if __name__ == "__main__": options, flags = grass.parser() sys.exit(main())