Rasterio: Numpy for Geospatial Data

Ciao, sono Alex di DEV STORIES DOT EU. Rasterio: Numpy per i dati geospaziali, episodio 1 di 5. Vuoi estrarre i valori dei pixel da un'immagine satellitare, ma le librerie di immagini standard si bloccano sui formati geospaziali e i binding ufficiali sembrano codice C scritto in Python. È proprio per questo che esiste Rasterio. Rasterio è una libreria che legge e scrive dati geospaziali basati su griglia, come i file GeoTIFF. È molto probabile che tu abbia già sentito parlare di GDAL, la libreria geospaziale standard del settore. Un malinteso comune è pensare che Rasterio sostituisca GDAL. Non è così. In realtà Rasterio usa GDAL sotto il cofano. La differenza sta nell'interfaccia. Se hai mai usato i binding nativi di GDAL per Python, sai che espongono una API C macchinosa, piena di gestione manuale delle risorse e con una sintassi che risulta completamente estranea a Python. Rasterio si lascia tutto questo alle spalle. Offre un'interfaccia moderna, Python-first. Tratta i dataset raster come normali file e gestisce i loro pixel come array Numpy standard. Supponi di avere un GeoTIFF di immagini Landsat sul disco. Apri questo file usando la funzione rasterio dot open. Dato che Rasterio abbraccia i pattern standard di Python, lo fai usando uno statement with. Questo crea un context manager. Proprio come quando apri un semplice file di testo, il blocco with garantisce che il dataset venga chiuso in modo pulito e che la memoria venga liberata non appena il tuo codice termina l'esecuzione. Non devi mai distruggere manualmente gli oggetti o preoccuparti dei memory leak causati da file non chiusi. All'interno di quel blocco with, la funzione open ti restituisce un oggetto dataset reader. Prima di leggere i pixel veri e propri, puoi ispezionare il file. Puoi chiedere al dataset il suo count, che ti dice il numero totale di bande, o layer, nell'immagine. Per un'immagine Landsat, questo count potrebbe essere qualcosa come sette o undici, che rappresentano diverse lunghezze d'onda della luce. Puoi leggere le property width e height per ottenere le dimensioni spaziali in pixel. Puoi anche controllare la property dtypes. Questo ti fornisce il data type dei pixel per ciascuna banda, come ad esempio interi senza segno a 16 bit o float a 32 bit. Tutti questi metadati sono disponibili all'istante, senza dover caricare i pesanti dati dell'immagine nella memoria di sistema. Ecco il punto chiave. Quando sei pronto a guardare i pixel veri e propri, chiami il metodo read sull'oggetto dataset. Puoi passare un indice specifico per caricare un solo layer. Fai attenzione a questo dettaglio: le bande di Rasterio sono one-indexed, il che significa che chiedi la banda uno, non la banda zero. Se chiami read e passi il numero uno, il metodo restituisce quei pixel come un array Numpy bidimensionale standard. Se chiami read senza alcun argomento, legge tutto, restituendo un array Numpy tridimensionale contenente tutte le bande. Non c'è nessuno speciale oggetto pixel proprietario. Ottieni semplicemente un array di numeri. Una volta ottenuto quell'array Numpy, sei di nuovo in un territorio familiare. Puoi fare slicing, eseguire operazioni statistiche o passarlo direttamente a modelli di machine learning usando l'esatto stesso codice che usi per qualsiasi altra matrice. Rasterio fa da ponte. Gestisce i complicati formati di file e i metadati su disco, e ti consegna un oggetto matematico pulito in memoria. Se ti piace lo show e vuoi supportare quello che facciamo, puoi cercare DevStoriesEU su Patreon: è un aiuto enorme. Il vero potere di Rasterio non è che reinventa il processing geospaziale, ma che lo rende noioso; trasformando immagini satellitari complesse in array Numpy standard, ti permette di smettere di combattere con i formati GIS e iniziare a fare vera data science. Grazie per l'ascolto, happy coding a tutti!

Ciao, sono Alex di DEV STORIES DOT EU. Rasterio: Numpy per i dati geospaziali, episodio 2 di 5. Un array Numpy di base è semplicemente una griglia piatta di numeri. Se lo visualizzi, ottieni semplicemente un rettangolo colorato sul tuo monitor. Cosa trasforma effettivamente questi numeri astratti in una mappa precisa dell'altopiano del Colorado? La risposta è l'Affine Transform. Spesso si pensa che i dati geospaziali siano semplicemente un'immagine in cui gli angoli sono vagamente ancorati alle coordinate di una mappa. Non è così. La transform è un collegamento matematico rigoroso che calcola la posizione reale di ogni singolo pixel nel tuo array. Prima di poter mappare i pixel sulla Terra, devi sapere come stai misurando la Terra. Questo è il Coordinate Reference System. Puoi leggerlo dal dataset accedendo a dataset punto crs. Questa property indica a Rasterio se le tue coordinate spaziali sono misurate in gradi di longitudine e latitudine, oppure in metri rispetto a uno specifico equatore e primo meridiano. Una volta impostato il sistema di riferimento, ti serve la logica di mappatura. La trovi controllando dataset punto transform. Questo restituisce una matrice di affine transform. Sembra algebra lineare pesante, ma in pratica è un set di sei numeri che descrivono l'impronta fisica della tua griglia. La matrice definisce le esatte coordinate spaziali x e y dell'angolo in alto a sinistra del pixel in alto a sinistra. Definisce la larghezza fisica di un pixel, ad esempio, esattamente trenta metri. Definisce l'altezza fisica di un pixel, che di solito è negativa perché le righe dell'immagine contano verso il basso, mentre le coordinate della mappa contano verso l'alto. La matrice contiene anche i valori di rotazione nel caso in cui il satellite fosse inclinato quando è stata catturata l'immagine. Ecco il punto chiave. Non devi mai moltiplicare a mano i valori di questa matrice. Rasterio ti offre dei metodi diretti per saltare avanti e indietro tra lo spazio dei pixel e lo spazio delle coordinate. Supponi di analizzare un'immagine e di trovare un'anomalia in una posizione specifica, diciamo riga cinquecento e colonna mille. Per scoprire dove si trova sul pianeta reale, passi quei due interi al metodo dataset punto xy. Rasterio passa la riga e la colonna attraverso la matrice affine e restituisce le esatte coordinate spaziali x e y per il punto centrale di quel pixel specifico. Puoi anche far girare questa macchina al contrario. Diciamo che hai le esatte coordinate spaziali di una posizione fisica. Magari è un canyon situato esattamente a cento chilometri a est e cinquanta chilometri a sud dell'origine della tua immagine. Devi estrarre i dati dei pixel sottostanti per quel canyon. Prendi i tuoi valori spaziali x e y e li passi al metodo dataset punto index. Rasterio inverte la matrice affine, processa le tue coordinate e ti restituisce gli interi esatti di riga e colonna. Usi quindi quegli interi per fare uno slice nel tuo array Numpy e leggere l'esatta misurazione registrata dal satellite. L'Affine Transform isola il tuo array di dati grezzi dalla tua geometria spaziale, garantendo che, non importa come fai il crop o leggi i tuoi dati, non perdi mai la tua posizione precisa sul globo. Per oggi è tutto. Grazie per l'ascolto — vai a costruire qualcosa di fantastico.

Ciao, sono Alex di DEV STORIES DOT EU. Rasterio: Numpy per i dati geospaziali, episodio 3 di 5. Se provi a leggere un GeoTIFF da dieci gigabyte direttamente nella memoria del tuo portatile, il tuo script andrà quasi certamente in crash. Magari ti serve solo una piccola porzione di pixel per testare un algoritmo, ma il sistema cerca di inghiottire l'intero continente in una volta sola. La soluzione a questo limite di memoria è il Windowed Processing. Il Windowed Processing ti permette di leggere o scrivere specifici sottoinsiemi rettangolari di un file raster, mentre il resto del file rimane al sicuro sul tuo hard drive. È un processo interamente basato sugli indici. Oggi non ci occupiamo di coordinate geografiche. Lavoriamo esclusivamente con righe e colonne di pixel. Per leggere un sottoinsieme, ti serve un modo per dire a Rasterio quali pixel prendere. Lo fai importando l'oggetto Window dal modulo windows di Rasterio. Il modo standard per definire una Window è fornire quattro numeri: column offset, row offset, width e height. Fai attenzione a questa parte. L'ordine degli argomenti è importante. Richiede prima il column offset, e poi il row offset. Il column offset è la tua posizione orizzontale di partenza, misurata dal bordo sinistro. Il row offset è la tua posizione verticale di partenza, misurata scendendo dal bordo superiore. Dopo gli offset, specifichi quanti pixel di larghezza e quanti di altezza deve avere la tua selezione. Immagina un enorme file di test da un gigabyte. Vuoi estrarre una piccola tile larga 256 pixel e alta 512 pixel. Vuoi che questa tile inizi a 1000 pixel dal bordo sinistro e a 2000 pixel dal bordo superiore. Istanzi la tua Window passando 1000 per il column offset, 2000 per il row offset, 256 per la width e 512 per la height. Dopodiché, apri il tuo dataset usando la funzione open standard. Quando chiami il metodo read, non lasciare gli argomenti vuoti. Gli argomenti vuoti dicono a Rasterio di caricare l'intero file in memoria. Invece, assegni il tuo nuovo oggetto Window al keyword argument window all'interno del metodo read. Rasterio traduce quella window in byte offset sul disco, recupera solo quello specifico chunk da 256 per 512, e lo restituisce come un array numpy. Il consumo della tua memoria si muove a malapena. Se lavori con numpy tutto il giorno, passare offset e dimensioni potrebbe sembrarti un po' innaturale. Probabilmente sei abituato a definire delle slice con gli indici di start e stop per le tue righe e colonne. Rasterio supporta questo workflow con un metodo alternativo chiamato Window dot from slices. Questo metodo prende due input. Il primo input è una tupla che definisce il row start index e il row stop index. Il secondo input è una tupla che definisce il column start index e il column stop index. Nota che qui l'ordine è prima le righe, poi le colonne. Questo rispecchia esattamente come fai lo slice di un array numpy bidimensionale. Dietro le quinte, questo metodo calcola gli offset e le dimensioni per te, e restituisce un oggetto Window standard che puoi passare al metodo read proprio come prima. Forzando delle windowed read nelle tue data pipeline, disaccoppi in modo permanente la dimensione fisica dei tuoi file di input dai limiti hardware della tua macchina. Questo è tutto per questo episodio. Grazie per l'ascolto, e continua a sviluppare!

Ciao, sono Alex di DEV STORIES DOT EU. Rasterio: Numpy per i dati geospaziali, episodio 4 di 5. Hai appena scaricato un'enorme immagine satellitare che copre un intero stato, ma ti interessano solo i pixel all'interno del confine di una specifica fattoria. Caricare l'intera griglia spreca memoria, e leggere tramite offset dei pixel è inutile quando il tuo confine è un poligono irregolare. La soluzione è il masking con i vettori. In un episodio precedente, abbiamo visto il windowed reading. Questo approccio richiedeva indici esatti dei pixel, come chiedere la riga cinquanta e la colonna cento. Il masking è completamente diverso. Il masking si basa su coordinate geografiche reali. Invece di contare i pixel, fornisci una geometria vettoriale, come un poligono definito da vere longitudini e latitudini. Rasterio capisce quali pixel rientrano in quella forma e gestisce la complessa conversione tra lo spazio continuo delle coordinate vettoriali e la griglia raster discreta. Lo strumento principale per questo workflow è la funzione mask, che si trova nel modulo rasterio punto mask. Le passi due argomenti principali. Primo, un dataset raster aperto. Secondo, un iterabile di forme geometriche. Queste forme devono essere formattate come dizionari che rispettano le specifiche GeoJSON. In genere, leggi uno shapefile che contiene il poligono della tua fattoria, estrai il dizionario della geometry per quella specifica feature, lo metti in una lista Python standard e lo passi alla funzione mask. Ecco il punto chiave. La funzione mask ha due comportamenti distinti, controllati da un singolo parametro booleano chiamato crop. Devi capire la differenza tra fare semplicemente il masking di un dataset e fare un vero e proprio crop. Se esegui la funzione con crop impostato su false, Rasterio guarda il poligono della tua fattoria e imposta ogni pixel fuori da quel confine a un valore nodata. Legge il valore nodata di default direttamente dai metadata del raster originale. L'array di output che ricevi ha esattamente le stesse dimensioni della tua enorme immagine originale dell'intero stato. I pixel all'interno della fattoria mantengono i loro valori reali, ma tutto ciò che è all'esterno viene oscurato. L'estensione spaziale complessiva del file non cambia. Se imposti crop su true, il comportamento cambia. Rasterio continua a impostare i pixel esterni al tuo poligono a nodata. Tuttavia, poi scarta le vaste aree di spazio vuoto. Calcola il bounding box rettangolare del poligono irregolare della tua fattoria e rimpicciolisce l'array restituito per adattarlo esattamente a quel box specifico. Il tuo gigantesco array dell'intero stato viene fisicamente ridotto a una piccola griglia gestibile, che contiene solo la fattoria e un po' di padding di nodata attorno ai suoi bordi irregolari. Dato che il crop modifica fisicamente le dimensioni della griglia, cambia anche il punto di partenza geografico. L'angolo in alto a sinistra non è più l'angolo in alto a sinistra dello stato. Ora è l'angolo in alto a sinistra del bounding box della fattoria. Di conseguenza, la funzione mask restituisce una tupla contenente due elementi. Il primo elemento è il tuo nuovo array numpy appena croppato. Il secondo elemento è un nuovissimo oggetto affine transform. Questa transform aggiornata contiene le nuove coordinate di origine. Se vuoi salvare questo array appena croppato in un nuovo file, devi passare questa transform aggiornata al tuo write profile. Se per errore riutilizzi la transform originale dell'intero stato, la tua piccola immagine della fattoria verrà stirata di nuovo su tutto lo stato. Il vero valore della funzione mask è che gestisce automaticamente tutti questi ricalcoli spaziali per te. Spero sia stato utile. Grazie per aver ascoltato, e buona continuazione di giornata.

Ciao, sono Alex di DEV STORIES DOT EU. Rasterio: Numpy per i dati geospaziali, episodio 5 di 5. A volte hai a disposizione bellissime immagini satellitari ad alta risoluzione, ma il tuo modello climatico accetta solo input grezzi e a bassa risoluzione. Non puoi semplicemente fare lo slice dell'array, perché i pixel risultanti devono comunque rappresentare l'esatta stessa area fisica, solo in una griglia a densità inferiore. Per farlo correttamente, hai bisogno di Resampling e Writing. Non si tratta di reprojection. Il coordinate reference system rimane esattamente lo stesso. Stiamo solo cambiando la risoluzione. Mettiamo che tu debba fare il downsample di un raster della metà. In Rasterio, esegui il resampling nell'esatto momento in cui leggi i dati in memoria. Quando chiami il metodo read sul tuo dataset, passi un parametro chiamato out shape. Invece di leggere l'intero array, Rasterio lo legge e lo scala immediatamente per adattarlo alla shape che hai richiesto. Per un raster a banda singola, la tua out shape sarebbe una tuple che specifica una banda, poi metà dell'altezza originale e metà della larghezza originale. Devi anche dire a Rasterio come calcolare i valori di questi nuovi pixel più grandi. Lo fai passando un parametro resampling. Rasterio fornisce un enum Resampling per questo. Se stai lavorando con dati continui come la temperatura, potresti passare Resampling dot bilinear, che calcola una media ponderata dei pixel circostanti. Se stai lavorando con dati categorici come la copertura del suolo, useresti Resampling dot nearest per preservare gli esatti valori di classe originali. Ecco il punto chiave. Ora hai un nuovo array Numpy più piccolo, ma rappresenta ancora l'esatta stessa area geografica sulla Terra. Poiché l'array ha meno pixel, ogni singolo pixel copre più spazio fisico. Questo significa che la tua Affine transform originale ora è completamente sbagliata. Se scrivi il nuovo array su disco usando la vecchia transform, la tua immagine si restringerà sulla mappa fino a coprire solo un quarto della sua area originale. Devi scalare la transform. Lo fai prendendo la transform del dataset originale e moltiplicandola per una Affine scaling matrix. Trovi gli scaling ratio dividendo la larghezza originale per la tua nuova larghezza, e l'altezza originale per la tua nuova altezza. Dato che abbiamo dimezzato le dimensioni, entrambi questi ratio sono esattamente due. Quando moltiplichi la transform originale per questi scaling factor, ottieni una nuova transform in cui le dimensioni dei pixel sono raddoppiate, ma il punto di origine in alto a sinistra rimane perfettamente ancorato. Ora hai il tuo array resampled e la tua transform scalata. Il passaggio finale è scrivere un GeoTIFF valido su disco. Per farlo, hai bisogno dei metadati. L'approccio più sicuro è copiare la property profile dal tuo dataset di origine. Questo ti dà un dictionary contenente il file driver, il coordinate reference system, il data type e il valore nodata. Quindi aggiorni questo dictionary profile con la tua nuova altezza, la tua nuova larghezza e la tua transform appena calcolata. Con il profile aggiornato pronto, chiami rasterio dot open. Fornisci il nuovo filename, imposti la modalità su write e fai l'unpack del tuo dictionary profile come keyword arguments. Questo crea un dataset vuoto su disco formattato secondo le tue esatte specifiche. Infine, chiami il metodo write su questo nuovo oggetto dataset e passi il tuo array Numpy downsampled. L'array viene scritto nel file, e il tuo nuovo GeoTIFF a bassa risoluzione è pronto per il modello climatico. Ricorda, modificare i dati raster è sempre un contratto in due parti: se cambi la shape dell'array Numpy, devi cambiare esplicitamente l'Affine transform per farla corrispondere, o la tua impronta spaziale sarà compromessa. Questo ci porta alla fine della nostra serie. Ti incoraggio a leggere la documentazione ufficiale di Rasterio per esplorare altri algoritmi di resampling e provare a costruire queste pipeline hands-on. Se vuoi suggerire argomenti per una serie futura, visita devstories dot eu. Vorrei prendermi un momento per ringraziarti per l'ascolto — ci aiuta tantissimo. Buona giornata!