Rasterio: Numpy for Geospatial Data

Hola, soy Alex de DEV STORIES DOT EU. Rasterio: Numpy para datos geoespaciales, episodio 1 de 5. Quieres extraer valores de píxeles de una imagen por satélite, pero las librerías de imágenes estándar se atragantan con los formatos geoespaciales, y los bindings oficiales son como escribir código C en Python. Esa tensión es exactamente la razón por la que existe Rasterio. Rasterio es una librería que lee y escribe datos geoespaciales basados en grid, como los archivos GeoTIFF. Es muy probable que hayas oído hablar de GDAL, la librería geoespacial estándar de la industria. Un error común es pensar que Rasterio reemplaza a GDAL. No es así. De hecho, Rasterio utiliza GDAL por debajo. La diferencia es la interfaz. Si alguna vez has usado los bindings nativos de GDAL para Python, sabrás que exponen una API de C engorrosa, llena de gestión manual de recursos y una sintaxis que resulta completamente ajena a Python. Rasterio deja eso atrás. Proporciona una interfaz moderna, Python-first. Trata los datasets ráster como archivos normales y maneja sus píxeles como arrays de Numpy estándar. Supongamos que tienes un archivo GeoTIFF de imágenes Landsat en tu disco. Abres este archivo usando la función rasterio punto open. Como Rasterio adopta patrones estándar de Python, lo haces usando un statement with. Esto crea un context manager. Al igual que al abrir un archivo de texto básico, el bloque with garantiza que el dataset se cierre limpiamente y que la memoria se libere en el momento en que tu código termina de ejecutarse. Nunca tendrás que destruir objetos manualmente ni preocuparte por memory leaks de archivos sin cerrar. Dentro de ese bloque with, la función open te da un objeto reader del dataset. Antes de leer los píxeles reales, puedes inspeccionar el archivo. Puedes pedirle al dataset su count, que te indica el número total de bandas, o capas, en la imagen. Para una imagen Landsat, este count podría ser algo así como siete u once, representando diferentes longitudes de onda de luz. Puedes leer las propiedades width y height para obtener las dimensiones espaciales en píxeles. También puedes comprobar la propiedad dtypes. Esto te da el data type de los píxeles para cada banda, como enteros sin signo de 16 bits o floats de 32 bits. Todos estos metadatos están disponibles al instante, sin cargar los pesados datos de la imagen en la memoria de tu sistema. Aquí está la clave. Cuando estés listo para ver los píxeles reales, llamas al método read en el objeto dataset. Puedes pasar un índice específico para cargar solo una capa. Presta atención a esta parte: las bandas de Rasterio son one-indexed, lo que significa que pides la banda uno, no la banda cero. Si llamas a read y pasas el número uno, el método devuelve esos píxeles como un array de Numpy bidimensional estándar. Si llamas a read sin ningún argumento, lo lee todo, devolviendo un array de Numpy tridimensional que contiene todas las bandas. No hay ningún objeto pixel especial y propietario. Simplemente obtienes un array de números. Una vez que tienes ese array de Numpy, vuelves a estar en territorio familiar. Puedes hacerle slice, ejecutar operaciones estadísticas o pasarlo directamente a modelos de machine learning usando exactamente el mismo código que usas para cualquier otra matriz. Rasterio actúa como un puente. Gestiona los formatos de archivo complejos y los metadatos en disco, y te entrega un objeto matemático limpio en memoria. Si disfrutas del programa y quieres apoyar lo que hacemos, puedes buscar DevStoriesEU en Patreon; es de gran ayuda. El verdadero poder de Rasterio no es que reinvente el procesamiento geoespacial, sino que lo hace aburrido; al convertir imágenes por satélite complejas en arrays de Numpy estándar, te permite dejar de pelear con los formatos GIS y empezar a hacer data science de verdad. ¡Gracias por escuchar, happy coding a todos!

Hola, soy Alex de DEV STORIES DOT EU. Rasterio: Numpy para datos geoespaciales, episodio 2 de 5. Un array básico de Numpy es simplemente una cuadrícula plana de números. Si lo visualizas, solo verás un rectángulo de color en tu monitor. ¿Qué convierte esos números abstractos en un mapa preciso de la Meseta de Colorado? La respuesta es la transformación afín. La gente suele pensar que los datos geoespaciales son solo una imagen donde las esquinas están fijadas de forma aproximada a las coordenadas del mapa. Eso no es así. La transformación es un enlace matemático estricto que calcula la posición exacta en el mundo real de cada píxel de tu array. Antes de poder mapear píxeles a la Tierra, tienes que saber cómo estás midiendo la Tierra. Este es el Sistema de Referencia de Coordenadas. Lo lees del dataset accediendo a dataset punto crs. Esta propiedad le dice a Rasterio si tus coordenadas espaciales se miden en grados de longitud y latitud, o en metros en relación a un ecuador y un meridiano de referencia específicos. Una vez configurado el sistema de referencia, necesitas la lógica de mapeo. La encuentras comprobando dataset punto transform. Esto te da una matriz de transformación afín. Suena a álgebra lineal dura, pero en realidad es un conjunto de seis números que describen la huella física de tu cuadrícula. La matriz define las coordenadas espaciales exactas x e y de la esquina superior izquierda del píxel de arriba a la izquierda. Define el ancho físico de un píxel, por ejemplo, exactamente treinta metros. Define la altura física de un píxel, que normalmente es negativa porque las filas de la imagen cuentan hacia abajo mientras que las coordenadas del mapa cuentan hacia arriba. La matriz también guarda valores de rotación por si el satélite estaba inclinado cuando se capturó la imagen. Aquí está la clave. Nunca tienes que multiplicar estos valores de la matriz a mano. Rasterio te da métodos directos para saltar de un lado a otro entre el espacio de píxeles y el espacio de coordenadas. Imagina que estás analizando una imagen y encuentras una anomalía en una ubicación específica, digamos la fila quinientos y la columna mil. Para averiguar dónde está eso en el planeta real, le pasas esos dos enteros al método dataset punto xy. Rasterio pasa la fila y la columna por la matriz afín y devuelve las coordenadas espaciales exactas x e y del punto central de ese píxel en concreto. También puedes hacer funcionar esta máquina a la inversa. Digamos que tienes las coordenadas espaciales exactas de una ubicación física. Quizás sea un cañón situado exactamente a cien kilómetros al este y cincuenta kilómetros al sur del origen de tu imagen. Necesitas extraer los datos de los píxeles subyacentes para ese cañón. Coges tus valores espaciales x e y, y los pasas al método dataset punto index. Rasterio invierte la matriz afín, procesa tus coordenadas y te devuelve los enteros exactos de la fila y la columna. Luego usas esos enteros para hacer un slice en tu array de Numpy y leer la medición exacta registrada por el satélite. La transformación afín aísla tu array de datos en crudo de tu geometría espacial, garantizando que, sin importar cómo recortes o leas tus datos, nunca pierdas tu posición precisa en el globo. Eso es todo por hoy. Gracias por escuchar, ve a crear algo genial.

Hola, soy Alex de DEV STORIES DOT EU. Rasterio: Numpy para datos geoespaciales, episodio 3 de 5. Si intentas leer un archivo GeoTIFF de diez gigabytes directamente en la memoria de tu portátil, es casi seguro que tu script crasheará. Puede que solo necesites un pequeño fragmento de píxeles para probar un algoritmo, pero el sistema intenta tragarse todo el continente a la vez. La solución a este límite de memoria es el Windowed Processing. El Windowed Processing te permite leer o escribir subconjuntos rectangulares específicos de un archivo raster mientras el resto del archivo se queda a salvo en tu disco duro. Esto está basado completamente en índices. Hoy no vamos a tratar con coordenadas geográficas. Vamos a tratar estrictamente con filas y columnas de píxeles. Para leer un subconjunto, necesitas una forma de decirle a Rasterio qué píxeles coger. Esto lo haces importando el objeto Window del módulo windows de rasterio. La forma estándar de definir un Window es proporcionando cuatro números: column offset, row offset, width y height. Presta atención a esta parte. El orden de los argumentos importa. Requiere el column offset primero, y luego el row offset. El column offset es tu posición horizontal inicial, midiendo desde el borde izquierdo. El row offset es tu posición vertical inicial, midiendo hacia abajo desde el borde superior. Después de los offsets, especificas cuántos píxeles de width y cuántos píxeles de height debería tener tu selección. Imagina un archivo de prueba masivo del tamaño de un gigabyte. Quieres extraer un pequeño tile que tiene 256 píxeles de width y 512 píxeles de height. Quieres que este tile empiece mil píxeles hacia adentro desde la izquierda y dos mil píxeles hacia abajo desde arriba. Instancias tu Window pasando mil para el column offset, dos mil para el row offset, 256 para el width y 512 para el height. A continuación, abres tu dataset usando la función open estándar. Cuando llamas al método read, no dejas los argumentos vacíos. Los argumentos vacíos le dicen a Rasterio que cargue todo el archivo en memoria. En su lugar, asignas tu nuevo objeto Window al keyword argument window dentro del método read. Rasterio traduce ese Window en byte offsets en el disco, recupera solo ese chunk específico de 256 por 512, y lo devuelve como un numpy array. Tu uso de memoria apenas se mueve. Si trabajas con numpy todo el día, pasar offsets y dimensiones puede resultarte un poco antinatural. Probablemente estés acostumbrado a definir slices con índices start y stop para tus filas y columnas. Rasterio soporta este workflow con un método alternativo llamado Window punto from slices. Este método recibe dos inputs. El primer input es una tuple que define el row start index y el row stop index. El segundo input es una tuple que define el column start index y el column stop index. Fíjate en que el orden aquí es filas primero, luego columnas. Esto refleja exactamente cómo haces un slice de un numpy array bidimensional. Por debajo, este método calcula los offsets y dimensiones por ti y devuelve un objeto Window estándar que puedes pasar al método read igual que antes. Al forzar windowed reads en todos tus data pipelines, desacoplas permanentemente el tamaño físico de tus archivos de entrada de los límites de hardware de tu máquina. Eso es todo por este episodio. Gracias por escuchar, y ¡sigue construyendo!

Hola, soy Alex de DEV STORIES DOT EU. Rasterio: Numpy para datos geoespaciales, episodio 4 de 5. Acabas de descargar una imagen de satélite enorme que cubre todo un estado, pero solo te interesan los píxeles dentro de los límites de una granja específica. Cargar todo el grid desperdicia memoria, y leer por pixel offsets es inútil cuando tu límite tiene forma de polígono irregular. La solución es el Masking con vectores. En un episodio anterior, vimos el windowed reading. Ese enfoque requería pixel indices exactos, como pedir la fila cincuenta y la columna cien. El Masking es completamente diferente. El Masking se basa en coordenadas geográficas del mundo real. En lugar de contar píxeles, le pasas una geometría vectorial, como un polígono definido por longitudes y latitudes reales. Rasterio averigua qué píxeles caen dentro de esa forma y se encarga de la compleja traducción entre el espacio de coordenadas vectoriales continuo y el grid raster discreto. La herramienta principal para este workflow es la función mask, ubicada en el módulo rasterio dot mask. Le pasas dos argumentos principales. Primero, un dataset raster abierto. Segundo, un iterable de formas geométricas. Estas formas deben estar formateadas como diccionarios que cumplan con la especificación GeoJSON. Normalmente, lees un shapefile que contiene el polígono de tu granja, extraes el diccionario de geometría para esa feature específica, lo metes en una lista estándar de Python y se lo pasas a la función mask. Aquí está la clave. La función mask tiene dos comportamientos distintos controlados por un único parámetro booleano llamado crop. Tienes que entender la diferencia entre simplemente hacer masking a un dataset y realmente hacerle un crop. Si ejecutas la función con crop en false, Rasterio mira el polígono de tu granja y asigna un valor nodata a cada píxel que quede fuera de ese límite. Lee el valor nodata por defecto directamente de los metadatos del raster original. El array de salida que recibes tiene exactamente las mismas dimensiones que tu enorme imagen original de todo el estado. Los píxeles de dentro de la granja conservan sus valores reales, pero todo lo de fuera se queda en blanco. La extensión espacial total del archivo no cambia. Si pones crop en true, el comportamiento cambia. Rasterio sigue asignando el valor nodata a los píxeles de fuera de tu polígono. Sin embargo, luego descarta las enormes áreas de espacio vacío. Calcula el bounding box rectangular del polígono irregular de tu granja y reduce el array devuelto para que encaje exactamente dentro de esa caja específica. Tu array gigante de todo el estado se reduce físicamente a un grid diminuto y manejable que contiene solo la granja y un poco de padding de nodata alrededor de sus bordes irregulares. Como hacer cropping cambia físicamente las dimensiones del grid, también cambia el punto de partida geográfico. La esquina superior izquierda ya no es la esquina superior izquierda del estado. Ahora es la esquina superior izquierda del bounding box de la granja. Por lo tanto, la función mask devuelve una tupla que contiene dos elementos. El primer elemento es tu nuevo array de numpy recortado. El segundo elemento es un objeto affine transform totalmente nuevo. Este transform actualizado contiene las nuevas coordenadas de origen. Si quieres guardar este array recién recortado en un archivo nuevo, debes pasar este transform actualizado a tu write profile. Si reutilizas por error el transform original de todo el estado, la imagen de tu pequeña granja se estirará por todo el estado. El verdadero valor de la función mask es que se encarga de todo este recálculo espacial por ti automáticamente. Espero que te haya sido útil. Gracias por escuchar, y disfruta del resto de tu día.

Hola, soy Alex de DEV STORIES DOT EU. Rasterio: Numpy para datos geoespaciales, episodio 5 de 5. A veces tienes imágenes de satélite preciosas y de alta resolución, pero tu modelo climático solo acepta inputs de baja resolución. No puedes simplemente hacer un slice del array, porque los píxeles resultantes todavía tienen que representar exactamente la misma área física, solo que en un grid de menor densidad. Para hacer eso correctamente, necesitas Resampling y Writing. Esto no va de reproyección. El coordinate reference system se queda exactamente igual. Solo estamos cambiando la resolución. Digamos que necesitas hacer un downsample de un raster a la mitad. En Rasterio, haces el resampling en el momento exacto en el que lees los datos en memoria. Cuando llamas al método read en tu dataset, pasas un parámetro llamado out shape. En lugar de leer el array completo, Rasterio lo lee y lo escala inmediatamente para que encaje en el shape que has pedido. Para un raster de una sola banda, tu out shape sería una tupla especificando una banda, luego la mitad del height original, y la mitad del width original. También tienes que decirle a Rasterio cómo calcular los valores de estos nuevos píxeles más grandes. Esto lo haces pasando un parámetro resampling. Rasterio proporciona un enum Resampling para esto. Si estás trabajando con datos continuos como la temperatura, podrías pasar Resampling punto bilinear, que calcula una media ponderada de los píxeles de alrededor. Si estás trabajando con datos categóricos como la cobertura del suelo, usarías Resampling punto nearest para preservar los valores de clase originales exactos. Aquí está la clave. Ahora tienes un nuevo array Numpy más pequeño, pero sigue representando exactamente la misma área geográfica en la Tierra. Como el array tiene menos píxeles, cada píxel individual cubre más espacio físico. Esto significa que tu Affine transform original ahora está completamente mal. Si escribes el nuevo array a disco usando el transform antiguo, tu imagen encogerá en el mapa para cubrir solo una cuarta parte de su área original. Tienes que escalar el transform. Esto lo haces cogiendo el transform del dataset original y multiplicándolo por una scaling matrix de Affine. Encuentras los scaling ratios dividiendo el width original entre tu nuevo width, y el height original entre tu nuevo height. Como hemos reducido las dimensiones a la mitad, ambos ratios son exactamente dos. Cuando multiplicas el transform original por estos scaling factors, obtienes un nuevo transform donde las dimensiones de los píxeles se duplican, pero el punto de origen arriba a la izquierda se queda perfectamente anclado. Ahora tienes tu array con el resampling y tu transform escalado. El paso final es escribir un GeoTIFF válido a disco. Para hacer esto, necesitas metadatos. El enfoque más seguro es copiar la property profile de tu dataset de origen. Esto te da un diccionario que contiene el file driver, el coordinate reference system, el data type, y el valor nodata. Luego actualizas este diccionario profile con tu nuevo height, tu nuevo width, y tu transform recién calculado. Con el profile actualizado listo, llamas a rasterio punto open. Le pasas el nuevo filename, pones el mode en write, y desempaquetas tu diccionario profile como keyword arguments. Esto crea un dataset vacío en disco formateado con tus especificaciones exactas. Finalmente, llamas al método write en este nuevo objeto dataset y le pasas tu array Numpy con el downsample. El array se escribe en el archivo, y tu nuevo GeoTIFF de menor resolución está listo para el modelo climático. Recuerda, modificar datos raster es siempre un contrato de dos partes: si cambias el shape del array Numpy, tienes que cambiar explícitamente el Affine transform para que coincida, o tu spatial footprint se romperá. Eso nos lleva al final de nuestra serie. Te animo a leer la documentación oficial de Rasterio para explorar otros algoritmos de resampling y a probar a construir estas pipelines de forma práctica. Si quieres sugerir temas para una futura serie, visita devstories punto eu. Me gustaría tomarme un momento para darte las gracias por escucharnos; nos ayuda mucho. ¡Que vaya genial!