Creazione di una mappa topografica del Nepal utilizzando Python

Mappa topografica Nepal con Python

Introduzione

Ti sei mai chiesto come la topografia del tuo paese influisce sullo sviluppo economico e politico? Le mappe topografiche, che sono mappe della superficie terrestre che utilizzano linee di contorno per la visualizzazione, possono aiutare a rispondere a queste domande! Utilizzeremo Python per creare una mappa topografica per il Nepal, un paese con un interessante ambiente topografico. Imparerai come leggere dati geospaziali che descrivono la topografia di un paese, come interpretare questi dati e come visualizzarli. La mappa risultante può essere combinata con altri dati di interesse a livelli subnazionali molto disaggregati per capire come la topografia di un paese influisce sul suo sviluppo economico e/o politico. Questo articolo ti insegnerà come generare uno strumento molto interessante che può informare le politiche e lo sviluppo del settore privato!

Obiettivi di Apprendimento

  • Acquisire competenze nelle tecniche di analisi dei dati per i dati di elevazione digitale.
  • Imparare come utilizzare i dati geospaziali e gli strumenti di analisi correlati in Python.
  • Acquisire conoscenze sulle tecniche di mappatura.
  • Sviluppare competenze nella visualizzazione efficace dei dati per la comunicazione.
  • Comprendere l’importanza dell’elevazione per l’ineguaglianza e la povertà.

Questo articolo è stato pubblicato come parte del Data Science Blogathon.

Cosa sono le Mappe Topografiche?

Le mappe topografiche sono mappe della superficie terrestre che utilizzano linee di contorno per la visualizzazione. Le mappe topografiche sono uno strumento prezioso per navigare in terreni sconosciuti e servono come input per la pianificazione urbana e la gestione delle catastrofi. Sono spesso utilizzate per comprendere il contesto spaziale delle politiche o dei progetti del settore privato relativi allo sviluppo delle infrastrutture, per identificare aree vulnerabili ai disastri naturali o con accesso limitato ai servizi essenziali, come l’istruzione, l’assistenza sanitaria e le infrastrutture, o per la gestione delle risorse naturali. In ultima analisi, queste mappe possono fungere da input per la presa di decisioni basata sulle evidenze. In questo articolo, utilizzeremo Python per creare una mappa topografica per il Nepal, un paese con un ambiente topografico davvero interessante.

Descrizione dei Dati

Per generare la nostra mappa, ci affideremo ai dati pubblicati dal United States Geological Survey (USGS). USGS è un’agenzia scientifica del governo federale degli Stati Uniti che genera dati e ricerche su risorse naturali, geologia, geografia, risorse idriche e rischi naturali. Per accedere alla loro pagina dei dati, digita “USGS Data” su Google o clicca sul link che ti indirizza al loro Earth Explorer. Earth Explorer è uno strumento online e un portale dati che ti consente di cercare, accedere e scaricare una vasta gamma di dati scientifici sulla Terra. È necessario creare un account e accedere per utilizzare pienamente i dati.

Download dei Dati

Questo articolo utilizzerà il Nepal come esempio a causa delle sue caratteristiche topografiche uniche. Il Nepal ha una delle topografie più sfidanti e interessanti al mondo. 8 delle 14 montagne oltre gli 8.000 metri si trovano in Nepal (Trekking Trail Nepal), e il paese è diviso in tre regioni topografiche molto diverse: le Montagne, le Colline e il Terai (o pianure) (DHS). Sebbene queste caratteristiche rendano il paese unico e interessante, alcune ricerche mostrano che la topografia del Nepal rende difficile collegare il paese, fornire servizi essenziali alla sua popolazione e imporre rischi e ostacoli a un percorso di sviluppo sostenibile.

A tal fine, effettueremo una ricerca per il Nepal nei Criteri di Ricerca, come indicato nell’immagine sottostante. Una volta selezionato il Nepal, abbiamo selezionato il nostro set di dati di interesse. Per farlo, clicca sulla scheda Data Sets e scegli Digital Elevation. Ci sono diverse opzioni per i dati di elevazione digitale e, sebbene tu possa utilizzare diversi di questi set di dati, utilizzeremo i dati Global Multi-resolution Terrain Elevation Data 2010 GMTED2010. Questi dati forniscono una copertura globale del terreno terrestre a diverse risoluzioni (da 7,5 secondi di arco (circa 250 metri) a 30 secondi di arco (circa 1 chilometro)). Sono generati da dati di rilevamento remoto satellitare e aereo, compresi l’altimetria satellitare, l’immagine stereo e le mappe topografiche.

Dopo aver scelto i dati, clicca sulla scheda Risultati. Ora puoi scaricare i dati cliccando sul simbolo con le opzioni di download. Puoi anche visualizzare i dati tramite l’icona della traccia. Scarichiamo i dati nella loro massima risoluzione (7,5 secondi di arco). È importante notare che, per coprire tutto il Nepal, è necessario scaricare due mosaici (parti) diverse dei dati sottostanti e combinarli successivamente. Vedrai che il set di dati risultante è in formato tif, che indica dati raster.

Prepara i Dati con Strumenti di Analisi Geospaziale in Python

Python fornisce diversi strumenti per l’analisi geospaziale. In questo post del blog, ci basiamo sulla libreria Rasterio che rende possibile leggere e scrivere dati raster geospaziali (dati a griglia). Iniziamo e leggiamo il primo mosaico (parte) dei dati che abbiamo precedentemente scaricato nel nostro Jupyter Notebook:

#importa le librerie pertinenti (dopo averle installate)
import rasterio
import matplotlib.pyplot as plt
import numpy as np

#Leggi i dati e mostra la forma del dataset
file = rasterio.open(r'percorso\10n060e_20101117_gmted_mea075.tif')
dataset = file.read()
print(dataset.shape)

Carichiamo anche il secondo mosaico e combiniamoli unendoli. A tal fine, seguiamo le tecniche standard di lettura e manipolazione dei dati raster in Python come segue:

#Carica il secondo dataset e mostra la forma del dataset
file2 = rasterio.open(r'percorso\30n060e_20101117_gmted_mea075.tif')
dataset2 = file2.read()
print(dataset2.shape)


#Combina entrambi i dataset
from rasterio.merge import merge
from rasterio.plot import show

#Crea una lista vuota
src_files_to_mosaic = []

#Aggiungi la lista con entrambi i file
src_files_to_mosaic.append(file)
src_files_to_mosaic.append(file2)
src_files_to_mosaic

#Unisci entrambi i file
mosaic, out_trans = merge(src_files_to_mosaic)

#Copia i metadati
output_meta = file.meta.copy()

#Aggiorna i metadati
output_meta.update(
    {"driver": "GTiff",
        "height": mosaic.shape[1],
        "width": mosaic.shape[2],
        "transform": out_trans,
    }
)

#Scrivi sulla destinazione
# Scrivi il raster del mosaico su disco
out_fp = r"percorso\Nepal_Mosaic.tif"

with rasterio.open(out_fp, "w", **output_meta) as dest:
        dest.write(mosaic)

#Apri i dati raster combinati
file_mosaic = rasterio.open(out_fp)

#Leggi i dati
dataset_mosaic = file_mosaic.read()
print(file_mosaic.shape)

#Mostra i dati
plt.imshow(dataset_mosaic[0], cmap='Spectral')
plt.show()

Dati Globali di Elevazione del Terreno a Risoluzione Multipla

Ora abbiamo un set di dati combinato Global Multi-resolution Terrain Elevation Data 2010 GMTED2010 per tutto il Nepal, ma il file copre anche ampie parti dell’area circostante che non fanno parte del Nepal. Limitiamo l’area al Nepal utilizzando uno shapefile del Nepal. Utilizzeremo uno shapefile con i confini dei paesi del mondo. Puoi scaricare questo dataset qui. Quindi ritagliamo i dati raster e lo shapefile utilizzando la funzione “mask”. Utilizzeremo solo la prima riga dello shapefile e la colonna di geometria. Il risultato di questa operazione viene memorizzato in “clipped_array”, che è il dato raster ritagliato, e “clipped_transform”, che rappresenta le informazioni di trasformazione del dato raster ritagliato.

import geopandas as gpd
from shapely.geometry import mapping
from rasterio import mask as msk

#Carica lo shapefile con i confini dei paesi del mondo
df = gpd.read_file(r'percorso/world-administrative-boundaries.shp')

#Limita al Nepal
nepal = df.loc[df.name=="Nepal"]
nepal.head()

#Ritaglia i dati
clipped_array, clipped_transform = msk.mask(file_mosaic, [mapping(nepal.iloc[0].geometry)], crop=True)

#

C’è un altro problema. I valori di “no data” nei dati raster sono molto dannosi. Pertanto, distorcono la visualizzazione della nostra mappa, poiché fanno parte dell’intervallo di valori.

Capire il Problema

Prendiamoci cura di questo problema come segue, come descritto in questo post del blog:

  • Costruiamo una funzione che si occupa dei valori di “no data”. Costruiamo un parametro “no-data” per specificare il valore considerato “no data” nell’array ritagliato. In questo caso, viene impostato su (np.amax(clipped_array[0]) + 1), il che significa che è uguale al valore massimo nell’array ritagliato più uno. Questo valore sarà considerato come il valore “no data”.
  • Regoliamo l’array ritagliato aggiungendo il valore assoluto del valore minimo nell’array ritagliato al primo canale (indice 0) dell’array ritagliato. Questo passaggio garantisce che tutti i valori nell’array ritagliato diventino non negativi.
  • Calcoliamo anche l’intervallo di valori dell’array ritagliato. Viene aggiunto il valore massimo e il valore assoluto del valore minimo nell’array ritagliato. La variabile “value_range” conterrà l’intervallo di valori calcolato.
  • Utilizziamo un dizionario di valori di colore costruito manualmente basato su uno già esistente (il dizionario “seismic”) e definiamo il nostro colore di sfondo per i valori di “no data”.
  • Nell’ultimo passaggio, plottiamo la mappa con il nuovo intervallo di colori chiamato “new_seismic”.
#Indaghiamo i valori di dati mancanti

nodata_value = file_mosaic.nodata 
print("Valore di dati mancanti:", nodata_value)
#Valore di dati mancanti: -32768.0

#Cambia il valore di dati mancanti con uno superiore alla massima elevazione
def clip_raster(gdf, img):
    clipped_array, clipped_transform = msk.mask(img, [mapping(gdf.iloc[0].geometry)], crop=True)
    clipped_array, clipped_transform = msk.mask(img, [mapping(gdf.iloc[0].geometry)],
                                                          crop=True, nodata=(np.amax(clipped_array[0]) + 1))
    clipped_array[0] = clipped_array[0] + abs(np.amin(clipped_array))
    value_range = np.amax(clipped_array) + abs(np.amin(clipped_array))
    return clipped_array, value_range

nepal_topography, value_range = clip_raster(nepal, file_mosaic)


#Verifica che ciò abbia funzionato
print(value_range)


#Diamo al valore di dati mancanti un nuovo colore di sfondo
from matplotlib import cm
from matplotlib.colors import ListedColormap,LinearSegmentedColormap

#Seismic
new_seismic = cm.get_cmap('seismic', 8828)

#Definisci il colore di sfondo
background_color = np.array([0.9882352941176471, 0.9647058823529412, 0.9607843137254902, 1.0])

#Usa la mappa dei colori
newcolors = new_seismic(np.linspace(0, 1, 8828))

# Aggiungi il colore di sfondo come ultima riga all'array newcolors.
newcolors = np.vstack((newcolors, background_color))

#Usa la nuova mappa dei colori per l'Italia
new_seismic = ListedColormap(newcolors)

#Crea la mappa finale e salvala
plt.figure(figsize=(10,10))
c = plt.imshow(nepal_topography[0], cmap = new_seismic)
clb = plt.colorbar(c, shrink=0.4)
clb.ax.set_title('Elevazione (metri)',fontsize=10)

plt.savefig(r'percorso\Topographic_Map_Nepal.png', bbox_inches='tight')
plt.show()

Voilá! Abbiamo una mappa topografica del Nepal che indica chiaramente le diverse elevazioni nel paese e le tre zone topografiche.

Conclusioni

Hai imparato a generare una mappa topografica in Python utilizzando dati geospaziali dell’United States Geological Survey (USGS). Hai anche appreso l’importanza di prendersi cura dei valori mancanti nell’insieme di dati finale per la visualizzazione.

Politici o professionisti possono ora utilizzare questa mappa per ulteriori analisi, come combinarla con altre mappe, come mappe della povertà o disastri naturali, per analizzare se esiste qualche connessione. Abbiamo generato uno strumento prezioso che può contribuire a decisioni politiche basate su evidenze!

Punti chiave

  • Le mappe topografiche sono strumenti utili per decisioni basate su evidenze.
  • Topografia ed elevazione svolgono un ruolo cruciale nella pianificazione urbana, nella fornitura di servizi e nell’ineguaglianza.
  • Python dispone di strumenti utili per analizzare dati geospaziali.
  • Prestare attenzione ai valori di dati mancanti in questo tipo di dati è cruciale per la visualizzazione.
  • La visualizzazione di dati geospaziali può generare informazioni preziose a livelli disaggregati.

Spero che tu abbia trovato questo articolo informativo. Non esitare a contattarmi su LinkedIn. Connettiamoci e lavoriamo insieme per sfruttare i dati per un cambiamento positivo.

Domande frequenti

I media mostrati in questo articolo non sono di proprietà di Analytics Vidhya e vengono utilizzati a discrezione dell’autore.