:: Python :: tin.lago :: |
|---|
|
Lo script tin.lago.py
(*) è una ulteriore applicazione estensiva dell'algoritmo
implementato in v.quota.py. Lo script, permette
di calcolare la profondità media e massima di un lago nonché il suo volume avendo in precedenza
delimitato il contorno dello stesso mediante una entità topologica areale.
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. Tornando ora, allo script tin.lago.py (*), il suo utilizzo è semplice, la finestra di dialogo richiede l'immissione di una copertura TIN, di un file vettoriale a topologia areale racchiudente il contorno del lago, la quota della superficie dell'acqua ed il valore della spaziatura di campinamento delle altezze puntuali (Delta campionamento). Lo schema seguito nel codice è il seguente, vengono aperti entrambi i file se esisteni e, mediante un doppio ciclo for, si procede al calcolo delle altezze nei punti campionati mediante la funzione Vect_tin_get_z(map_info, x, y, byref(z), None, None) controllando, preventivamente, che questi ricadano all'interno del bacino mediante la funzione Vect_point_in_area(map_info_b, 1, x, y) inclusa in una istruzione if e che inoltre abbiano quota minore o uguale a quella della superficie del lago. 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 3 minuti.
#!/usr/bin/env python
#****************************************************************
#*
#* MODULE: tin.lago, v 1.0.0
#*
#* AUTHOR(S): Antonio Alliegro
#*
#* PURPOSE: Calcolo profondita' media e massima di un lago
#* mediante campionamento
#*
#* COPYRIGHT: (C) 2012 Antonio Alliegro Civil Engineer
#* Salerno, Italy
#* antonioall(at)libero.it
#*
#* First Version: 2012/03/12
#* Last Version: 2012/03/12
#*
#* 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 profondita' media e massima di un lago mediante campionamento
#% keywords: TIN, vector
#%end
#%option
#% key: map
#% type: string
#% gisprompt: old,vector,vector
#% key_desc: tin
#% description: Nome copertura TIN
#% required: yes
#%end
#%option
#% key: lago
#% type: string
#% gisprompt: old,vector,vector
#% key_desc: vector
#% description: Nome copertura contorno lago
#% required: yes
#%end
#%option
#% key: quota
#% type: double
#% answer: 400.0
#% key_desc: float
#% description: Quota superfice lago (m s.m.)
#% 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():
input = options['map']
basin = options['basin']
delta = float(options['delta'])
hmin = 99999.9
hmed = 0.0
n = 0
# Inializza la libreria GRASS
G_gisinit('')
# Controlla se le coperture sono nel percorso di ricerca
mapset = G_find_vector2(input, "")
if not mapset:
grass.fatal("Il file vettoriale <%s> non e' stato trovato" % input)
mapset_b = G_find_vector2(basin, "")
if not mapset_b:
grass.fatal("Il file vettoriale <%s> non e' stato trovato" % basin)
# Definisce la struttura delle coperture
map_info = pointer(Map_info())
map_info_b = 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, input, mapset)
Vect_open_old(map_info_b, basin, mapset_b)
# Verifica se il TIN e' un vettoriale 3D
if Vect_is_3d(map_info):
grass.message("Il file vettoriale <%s> e' 3D" % input)
else:
grass.fatal("Il file vettoriale <%s> non e' 3D" % input)
# Verifica se il vettoriale del bacino ha tipologia areale
nareas = Vect_get_num_areas(map_info_b)
if nareas == 0:
grass.fatal("Il file vettoriale <%s> non ha tipologia areale" % basin)
# 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
grass.message("Passo 1/2: Pre-elaborazione TIN...")
z = c_double()
G_percent(0, nx, 2)
Vect_tin_get_z(map_info, xref, yref, byref(z), None, None)
grass.message("...Fatto")
grass.message("Passo 2/2: Calcolo quote...")
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_b, 1, x, y) == 1:
Vect_tin_get_z(map_info, x, y, byref(z), None, None)
if z.value <= quota:
hmed += z.value
if z.value < hmin:
hmin = z.value
n += 1
G_percent(i, nx, 2)
grass.message("...Fatto")
# Risultati
if n > 0:
hmed = hmed / n
grass.message("----------------------------------")
grass.message("Delta campionamento = %12.3f" % delta)
grass.message("Numero celle totali = %12d" % (nx * ny))
grass.message("Numero celle lago = %14d" % n)
grass.message("Profondita' media =%15.3f" % (quota - hmed))
grass.message("Profondita' massima =%13.3f" % (quota - hmin))
grass.message("Volume acqua lago =%15.3f" % ((quota - hmed) * n * delta * delta))
# Chiusura coperture
Vect_close(map_info)
Vect_close(map_info_b)
return 0
#End main
if __name__ == "__main__":
options, flags = grass.parser()
sys.exit(main())
Visualizza il file sorgente di tin.lago.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:
|
|
