Más

Extracción de la proporción de clases en píxeles de un ráster más fino en uno más grueso


Tengo dos rásteres, uno tiene una resolución de 1 km. El otro tiene una resolución espacial de 30 m. Este último pixelwise es una clasificación. Quiero superponer el ráster de 1 km en el ráster de clasificación y extraer para cada píxel de 1 km la distribución de clases correspondientes a la clasificación de 30 m.

La clasificación de 30 m incluye 7 clases, por lo que me gustaría que este proceso produjera 7 nuevas capas a 1 km. Uno para cada clase. Cada píxel de 1 km ahora tendrá asociado un% correspondiente de la presencia de cada clase presente en la clasificación de 30 m. Si uno hiciera una tabla donde cada fila es un píxel, entonces tendría una tabla como la siguiente:

Una dificultad adicional es que los rásteres son bastante grandes. Realmente me gustaría lograr esto usando Python.

EDITAR:

He agregado subconjuntos de mis dos conjuntos de datos en este enlace (a resoluciones de 30 my 1000 m respectivamente), en caso de que alguien quiera proponer una prueba de concepto; sería increíble.


Paso 1

Cree rásteres de bits para cada una de las clases únicas. Puede ser un ráster de 1 banda para cada clase o un ráster único con una banda para cada clase (por ejemplo, GeoTIFF). Si usa GTiff, puede usar la opción de creaciónNBITS = 1para conservar espacio. También puede considerar rásteres de dos bits para almacenar lógica de tres valores donde el tercero (por ejemplo, 2) es NODATA, lo que requeriríaNBITS = 2.

Para cada banda / clase, los píxeles serán 0 donde esa clase no existe, o 1 donde existe esa clase.

import numpy as np from osgeo import gdal gdal.UseExceptions () # Obtener datos de ráster con clasificaciones ds = gdal.Open ('cropped_30m.tif') band = ds.GetRasterBand (1) class_ar = band.ReadAsArray () gt = ds .GetGeoTransform () pj = ds.GetProjection () ds = band = None # close # Define los valores ráster para cada clase, para relacionarlos con cada banda class_ids = (np.arange (31) + 1) .tolist () # Make un nuevo bit ráster drv = gdal.GetDriverByName ('GTiff') ds = drv.Create ('bit_raster.tif', class_ar.shape [1], class_ar.shape [0], len (class_ids), gdal.GDT_Byte, [ 'NBITS = 1', 'COMPRESS = LZW', 'INTERLEAVE = BAND']) ds.SetGeoTransform (gt) ds.SetProjection (pj) para bidx en rango (ds.RasterCount): band = ds.GetRasterBand (bidx + 1 ) # crear selección booleana = (class_ar == class_ids [bidx]) band.WriteArray (selection.astype ('B')) ds = band = None # guardar, cerrar

P.ej. para Rojo = Banda 10, Verde = Banda 13 y Azul = Banda 27, se puede visualizar:

El negro indica que es otra clase / banda.

Paso 2

Utilice gdalwarp para obtener el promedio de valores 0 y 1 de una grilla de muestra más gruesa. Este método requiere GDAL> = 1.10.0. Una forma realmente sencilla de hacer esto es especificar una resolución de destino en la línea de comando, p. Ej. 1 km:

$ gdalwarp -ot Float32 -tr 1000 1000 -r promedio bit_raster.tif promedio_1km.tif

El resultado tendrá tantas bandas como clases, y cada banda describirá la fracción de esa clase para cada píxel, entre 0 y 1. Si prefiere el porcentaje, multiplíquelo por 100 (p. Ej., Consulte gdal_calc.py).

También hay una interfaz Python para gdalwarp (por ejemplo, esta respuesta), que puede usar un ráster de plantilla para la forma y donde se debe usar NODATA. Por ejemplo:

# Abrir ráster del paso 1 src_ds = gdal.Open ('bit_raster.tif') # Abrir una plantilla o copiar matriz, para dimensiones y máscara NODATA cpy_ds = gdal.Open ('cropped_1000m.tif') band = cpy_ds.GetRasterBand (1 ) cpy_mask = (band.ReadAsArray () == band.GetNoDataValue ()) # Ráster de resultado, con la misma resolución y posición que el ráster de copia dst_ds = drv.Create ('average2_1km.tif', cpy_ds.RasterXSize, cpy_ds.RasterYSize, len (class_ids), gdal.GDT_Float32, ['INTERLEAVE = BAND']) dst_ds.SetGeoTransform (cpy_ds.GetGeoTransform ()) dst_ds.SetProjection (cpy_ds.GetProjection ()) # Haga lo mismo que gdalwarp -r promedio; esto puede tardar un poco en finalizar gdal.ReprojectImage (src_ds, dst_ds, None, None, gdal.GRA_Average) # Convierta todas las fracciones a porcentaje y aplique la misma máscara NODATA del raster de copia NODATA = -9999 para bidx en rango (dst_ds .RasterCount): band = dst_ds.GetRasterBand (bidx + 1) ar = band.ReadAsArray () * 100.0 ar [cpy_mask] = NODATA band.WriteArray (ar) band.SetNoDataValue (NODATA) # Guarde y cierre todos los rásteres src_ds = cpy_ds = dst_ds = band = Ninguno

P.ej. las mismas tres bandas (clases) se pueden visualizar como RGB:

O introdúzcalo en QGIS y use la herramienta Identificar características para inspeccionar la fracción de cada banda / clasificación en una sola ubicación:

Estos números suman 100%, lo cual es un buen control.


Ver el vídeo: ArcGIS Raster Calculator (Septiembre 2021).