Cómo agregar rásteres de diferentes tamaños en GDAL para que el resultado sea solo en la región intersectada

11

Estoy escribiendo un método de Python que agrega dos rásteres y genera una única salida de ráster. Por razones que escapan a mi control, la extensión de los rásteres de entrada es diferente, pero se superponen.

¿Es posible operar exclusivamente en los píxeles de ráster de entrada que se superponen en las 2 regiones superpuestas para generar mi salida de modo que la extensión del ráster de salida sea solo la región de intersección de las dos entradas?

python raster gdal Rico
fuente

Respuestas:

12

Lo primero que debe hacer es determinar el rectángulo superpuesto en coordenadas geoespaciales. Para hacer esto, obtienes la geotransformación para cada imagen de origen:

gt1 = ds1.GetGeoTransform()
# r1 has left, top, right, bottom of dataset's bounds in geospatial coordinates.
r1 = [gt1[0], gt1[3], gt1[0] + (gt1[1] * ds1.RasterXSize), gt1[3] + (gt1[5] * ds1.RasterYSize)]

# Do the same for dataset 2 ...

intersection = [max(r1[0], r2[0]), min(r1[1], r2[1]), min(r1[2], r2[2]), max(r1[3], r2[3])]

Luego convierta ese rectángulo en píxeles para cada imagen restando las coordenadas superior e izquierda y dividiendo por el tamaño del píxel, redondeando hacia arriba.

Desde aquí puede llamar ReadRaster()a cada imagen, dándole las extensiones de píxeles que acaba de calcular:

band.ReadRaster(px1[0], px1[1], px1[2] - px1[0], px1[3] - px1[1], px1[2] - px1[0], px1[3] - px1[1],
                # <band's datatype here>
)

Estoy un poco cansado, así que si esto no tiene mucho sentido, ¡házmelo saber!

MerseyViking
fuente
¿Esto también funciona para rásteres con diferentes proyecciones (también conocidos como sistemas de referencia de coordenadas / sistemas de referencia espacial)? E incluso si las proyecciones son las mismas: ¿Esto también funciona si gt1[1]y gt2[1](o gt1[5]y gt2[5]) tienen signos opuestos? (Creo que sería voltear uno de los rásteres vertical u horizontalmente). ¿O si abs(gt1[2])y abs(gt1[4])son mayores que abs(gt1[1])y abs(gt1[5])pero abs(gt2[2])y abs(gt2[4])son más pequeños que abs(gt2[1])y abs(gt2[5])(lo que probablemente voltearía uno de los rásteres en diagonal)?
das-g
6

El tercer elemento de intersección debe ser min (r1 [2], r2 [2]):

intersection = [max(r1[0], r2[0]), min(r1[1], r2[1]), min(r1[2], r2[2]), max(r1[3], r2[3])]

Además, recomendaría algo de lógica para verificar que los conjuntos de datos realmente se crucen.

Este es un enfoque:

if (intersection[2] < intersection[0]) or (intersection[1] < intersection[3]):
    intersection = None
David Shean
fuente