.il primo giorno di GIS, tra ombre e rilievi

Come faccio a dimenticare il giorno in cui per la prima volta mi avvicinai alla Geomatica, era l’Aprile del 2002, il mese dopo sarei diventato uno Scienziato della Terra. Ricordo che la prima cosa che ho imparato è stato creare un DTM (Digital Terrain Model) a partire da un piano quotato e da curve di livello. Rimasi affascinato di come da questo semplice dato raster si potevano ottenere molteplici dati derivati: pendenze, esposizioni ed anche una mappa delle ombre, più comunemente conosciuto come shaded relief o hillshade.

In modo particolare mi è piaciuto l’hillshade, sarà stato l’effetto 3D generato dalle ombre o forse il fatto che si riusciva a dare una lettura più chiara sulla morfologia del terreno, non saprei, ma fra tutti i dati derivati da un DTM, il modello delle ombre ha ricevuto sempre un’attenzione particolare da parte mia. Crearlo era ed è abbastanza semplice, oltre al DTM come input il modello delle ombre richiede alcuni parametri necessari, parametri che dipendono dalla posizione del sole ed in modo specifico dalla sua altezza rispetto all’orizzonte.

In QGIS, per esempio, esistono attualmente almeno tre modi per creare un hillshade:

  1. attraverso il plugin Terrain Analysis;
  2. usando il plugin Processing e GDAL come backend;
  3. usando il plugin Processing e SAGA come backend.

In tutte e tre i casi è necessario conoscere ed inserire l’azimuth e l’altitudine del sole espressi in angoli. Ma se volessimo creare il nostro hillshade sulla base della posizione del sole in un determinato giorno e ad una determinata ora? Basta conoscere azimuth e altitudine del sole per quella data! Come?

Schermata 2016-04-08 alle 22.20.15

In letteratura si trovano diversi modelli per il calcolo della posizione del sole, io mi sono limitato ad utilizzarne uno già bello e pronto in una libreria python. La libreria si chiama PySolar. Tra le diverse funzioni presenti nel modulo, ce ne sono due in particolare che consentono di ricavare l’azimuth e l’altitudine del sole sulla base di una data temporale precisa. Quindi è sufficiente inserire questi due valori in uno dei tool presenti in QGIS.

Vediamo come fare questa semplice procedura attraverso l’ausilio di python. Scarichiamo la libreria PySolar e, come al solito, copiamo il modulo in una directory sulla nostra macchina o direttamente nella cartella {utente}/.qgis2/python.

Apriamo la Console Python di QGIS ed incolliamo il codice seguente per importare le librerie necessarie.

import sys
sys.path.append('/Users/larosa/dev/ext-libs')

import processing
from pysolar.solar import GetAltitude, GetAzimuth
import datetime

Ora definiamo i parametri di input necessari per ottenere la posizione del sole in un determinato giorno dell’anno e ad una specifica posizione geografica

DEM = '/Users/larosa/test/dtm.tif'

# date
DAY = 10
MONTH = 4
YEAR = 2002

# time
H = 11
M = 00

NOW = False

Adesso utilizziamo la libreria PySolar per recuperare i valori di azimuth e altitude del sole rispetto ai dati di input impostati precedentemente. Tali valori vengono ottenuti utilizzando come location le coordinate geografiche del centro della mappa

def get_solar_params():
    mc = iface.mapCanvas()
    s_srs = mc.mapSettings().destinationCrs()
    d_srs = QgsCoordinateReferenceSystem()
    d_srs.createFromSrid(4326)

    # initialize qgis class to convert the coordinates
    tr = QgsCoordinateTransform(s_srs, d_srs)

    # get center coords of the map canvas
    center = mc.center()
    point = QgsPoint(center)

    # convert coords from s_srs to d_srs
    point_tr = tr.transform(point)

    if NOW:
        d = datetime.datetime.utcnow()
    else:
        d = datetime.datetime(YEAR, MONTH, DAY, H, M)

    params = {'altitude': GetAltitude(point_tr.y(), point_tr.x(), d),
              'azimuth': GetAzimuth(point_tr.y(), point_tr.x(), d) - 180}

    # returns a dictionary with azimuth and altitude values
    return params

Se la variabile NOW è impostata a True, il calcolo dell’hillshade viene fatto alla data odierna.

Il passo successivo consiste nel memorizzare in due variabili separate il valore dell’azimuth e dell’altitudine

altitude = get_solar_params()['altitude']
azimuth = get_solar_params()['azimuth']

Ora non resta che eseguire l’algoritmo per la creazione dell’hillshade. Il plugin Processing mette a disposizione anche un’interfaccia Python per eseguire gli algoritmi utilizzando la Console Python in un modo molto semplice. Per prima cosa ricerchiamo il nostro tool attraverso la funzione alglist(“nome_del_tool”)

>>> processing.alglist("hillshade")
 Hillshade-------------------------------------------->gdalogr:hillshade

Per sapere invece quali parametri sono necessari per eseguire il tool usiamo la funzione alghelp(“nome_del_tool”)

>>> processing.alghelp("gdalogr:hillshade")
ALGORITHM: Hillshade
 INPUT <ParameterRaster>
 BAND <ParameterNumber>
 COMPUTE_EDGES <ParameterBoolean>
 ZEVENBERGEN <ParameterBoolean>
 Z_FACTOR <ParameterNumber>
 SCALE <ParameterNumber>
 AZIMUTH <ParameterNumber>
 ALTITUDE <ParameterNumber>
 OUTPUT <OutputRaster>

Ed infine, eseguiamo l’algoritmo gdalogr:hillshade, utilizzando la funzione runanload(), per la creazione dell’hillshade, passando i valori ottenuti con PySolar sulla posizione della nostra stella .

if altitude > 0:
    processing.runandload('gdalogr:hillshade',
                          DEM, 1, False, False, 1, 1,
                          abs(azimuth),
                          abs(altitude), None)
else:
    print('no shadow')

Il risultato finale sarà caricato automaticamente in QGIS.

Schermata 2016-04-07 alle 17.10.15

 

Advertisements

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s