Programmazione & GIS

:: Python :: tin.to.raster ::
Nello script v.quota.py si calcola la quota di un singolo punto appartenente ad un vettoriale 3D. Questa tipologia di copertura, se è realizzata mediante facce triangolari, è detta TIN, per apprendere come realizzarla utilizzando Grass, si faccia riferimento alle seguenti pagine: Costruzione TIN da Punti Quotati, Costruzione TIN da Curve di Livello e Costruzione TIN da Entità Quotate.

Si analizza ora, lo script tin.to.raster.py (*) sviluppato per convertire i TIN in coperture vettoriali; all'utente è sufficiente indicare il nome di un TIN e il nome di un nuovo raster nelle cui celle, durante l'elaborazione, vengono inseriti i valori delle quote calcolate. L'applicativo preleva in automatico i valori delle impostazioni della regione corrente e quindi il numero di righe e colonne, la risoluzione lungo x e y, ecc. Noti questi valori, in un doppio ciclo for dei quali, il più l'esterno incrementa le ordinate e il più interno le ascisse, si calcola la quota da assegnare alla singola cella mediante la funzione Vect_tin_get_z(map_info, x, y, byref(z), None, None).

Si tenga presente che la conversione richiede diversi minuti, ad esempio, per un raster di circa 10.000.000 di celle il tempo di elaborazione si aggira intorno ai 12 minuti.

tin.to.raster.py 

#!/usr/bin/env python
#****************************************************************
#*
#* MODULE:     tin.to.raster, v 1.0.3
#*
#* AUTHOR(S):  Antonio Alliegro
#*
#* PURPOSE:    Converte coperture TIN in raster
#*
#* COPYRIGHT:  (C) 2011-2012 Antonio Alliegro Civil Engineer
#*             Salerno, Italy
#*             antonioall(at)libero.it
#*
#*             First Version: 2011/09/08
#*             Last  Version: 2012/09/29
#*
#*             This program is free software under the
#*             GNU General Public License (>=v2).
#*             Read the file COPYING that comes with GRASS
#*             for details.
#*
#****************************************************************

#%module
#% description: Converte coperture TIN in raster
#% keywords: TIN, vector, raster
#%end

#%option
#% key: map
#% type: string
#% gisprompt: old,vector,vector
#% key_desc: tin
#% description: Nome copertura TIN
#% required: yes
#%end

#%option
#% key: output
#% type: string
#% gisprompt: new,cell,raster
#% key_desc: raster
#% description: Nome copertura raster
#% 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']
    output = options['output']

    
    # Inializza la libreria GRASS
    G_gisinit('')

    
    # Controlla se il vettoriale e' nel percorso di ricerca
    mapset = G_find_vector2(input, "")
    if not mapset:
        grass.fatal("Il file vettoriale <%s> non e' stato trovato" % input)

    # Definisce la struttura della copertura
    map_info = pointer(Map_info())

    # Setta la copertura al livello 2 (level 2: topology)
    Vect_set_open_level(2)

    # Apre il file vettoriale
    Vect_open_old(map_info, input, mapset)

    # Verifica se il vettoriale e' 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)

    # Allocazione del buffer di output utilizzando i 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

    outrast = []
    for i in range(nrows):
        outrast[i:] = [G_allocate_d_raster_buf()]

    # Creazione raster
    outfd = G_open_raster_new(output, DCELL_TYPE)
    if outfd < 0:
        grass.fatal("Impossibile creare il file raster <%s>" % output)

    # Inserisce valori nulli nelle celle
    grass.message("Passo 1/4: Inserimento valori nulli nelle celle...")
    for i in range(nrows):
        G_set_d_null_value(outrast[i], ncols)
        G_percent(i, nrows, 2)
    grass.message("...Fatto")

    # Blocco principale
    grass.message("Passo 2/4: Pre-elaborazione TIN...")
    z = c_double()
    G_percent(0, nrows, 2)
    Vect_tin_get_z(map_info, xref, yref, byref(z), None, None)
    grass.message("...Fatto")

    grass.message("Passo 3/4: Conversione TIN in raster...")
    for i in range(nrows):
        for j in range(ncols):
            x = xref + j * xres
            y = yref + i * yres
            if Vect_tin_get_z(map_info, x, y, byref(z), None, None) > 0:
                outrast[i][j] = z
        G_percent(i, nrows, 2)
    grass.message("...Fatto")

    grass.message("Passo 4/4: Scrittura file raster...")
    for i in range(nrows - 1, -1, -1):
        if G_put_d_raster_row(outfd, outrast[i]) < 0:
            grass.fatal("Errore nella scrittura della copertura raster <%s>" % output)
        G_percent(nrows - i, nrows, 2)
    grass.message("...Fatto")

    # Rilascio buffer
    for i in range(nrows):
        G_free(outrast[i])

    # Chiusura raster
    G_close_cell(outfd)

    # Chiude la copertura
    Vect_close(map_info)
    grass.message("Elaborazione terminata")

    return 0
    #End main

if __name__ == "__main__":
    options, flags = grass.parser()
    sys.exit(main())

Visualizza il file sorgente di tin.to.raster.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: