:: Python :: tin.hmed.basin :: |
|---|
|
Lo script tin.hmed.basin.py
(*) è una applicazione estensiva dell'algoritmo
implementato in v.quota.py. Lo script, permette
di calcolare la quota media, minima e massima di un bacino idrografico e per estensione di una
qualsiasi entità topologica areale sovrapponibile ad una copertura TIN.
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.hmed.basin.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 bacino idrografico e del 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. 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.hmed.basin, v 1.0.0
#*
#* AUTHOR(S): Antonio Alliegro
#*
#* PURPOSE: Calcolo hmed, hmin, hmax di un bacino mediante
#* campionamento
#*
#* COPYRIGHT: (C) 2011 Antonio Alliegro Civil Engineer
#* Salerno, Italy
#* antonioall(at)libero.it
#*
#* First Version: 2011/09/14
#* Last Version: 2011/09/14
#*
#* 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 hmed, hmin, hmax di un bacino 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: basin
#% type: string
#% gisprompt: old,vector,vector
#% key_desc: vector
#% description: Nome copertura bacino
#% 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
hmax = -hmin
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)
hmed += z.value
if z.value < hmin:
hmin = z.value
if z.value > hmax:
hmax = 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 = %.3f" % delta)
grass.message("Numero celle totali = %d" % (nx * ny))
grass.message("Numero celle bacino = %d" % n)
grass.message("hmed =%15.3f" % hmed)
grass.message("hmin =%15.3f" % hmin)
grass.message("hmax =%15.3f" % hmax)
# 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.hmed.basin.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:
|
|
