:: Python :: tin.sterro_riporto :: |
|---|
|
Lo script tin.sterro_riporto.py
(*) è una estensione dell'algoritmo sviluppato in
v.quota.py. L'applicativo, è in grado di calcolare
i volumi di sterro e riporto di una qualsiasi entità topologica areale sovrapponibile a due
copertura TIN rappresentanti una lo stato di fatto e l'altra quello di progetto.
Per apprendere come realizzare un TIN con Grass, si faccia riferimento alle seguenti pagine: Costruzione TIN da Punti Quotati, Costruzione TIN da Curve di Livello e Costruzione TIN da Entità Quotate. Per utilizzare lo script tin.sterro_riporto.py (*), occorre inserire nella finestra di dialogo un file vettoriale a topologia areale e due coperture TIN, nonché il valore della spaziatura di campinamento delle altezze puntuali (Delta campionamento). L'algoritmo prevede l'apertura dei file delle tre coperture se esistenti e, conseguentemente mediante un doppio ciclo for, procede al calcolo delle differenze di altezza nei punti campionati mediante l'applicazione della funzione Vect_tin_get_z(map_info, x, y, byref(z), None, None) ai due TIN controllando, preventivamente, che gli stessi punti ricadano all'interno dell'area mediante la funzione Vect_point_in_area(map_info_b, 1, x, y) inclusa in una istruzione if. Mediante un'altra istruzione if posta all'interno del precedente if, si procede a sommare le differenza di altezza nelle variabili sterro e riporto a seconda che il valore calcolato della differenza di altezza sia rispettivamente negativo o positivo. A chiusura dei due cicli vengono infine ottenuti i volumi di sterro e riporto moltiplicando queste stesse variabili per la superficie occupata da una cella (delta * delta). Ovviamente la precisione dei risultati ed i tempi di elaborazione sono proporzionali al valore della spaziatura di campionamento impostata (Delta campionamento). Una griglia con 10.000.000 di celle richiede un tempo di elaborazione di circa 6 minuti.
#!/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())
Visualizza il file sorgente di tin.sterro_riporto.py NOTA (*): La libreria grass.lib è disponibile a partire dalla versione 6.4.2 di Grass. Pertanto il presente applicativo non è utilizzabile nelle versioni 6.4.1, 6.4.0 e precedenti.
Argomenti correlati:
|
|
