Sumando rásteres que se solapan sólo en cierta extensión con rasterio

Cuando los rásteres no se solapan adecuadamente en toda su extensión no es posible hacer álgebra de mapas porque al estar almacenados los values en arrays de tipo matricial estas operaciones no pueden realizarse si las matrices no son conformes. Sin embargo, rasterio tiene disponible el método 'Window' que permite la lectura y escritura por ventanas de datos. Si se determina cuáles son las áreas coincidentes es posible hacer la lectura de los arrays en esas zonas específicas y aplicar sólo allí el álgebra de mapas.

Una forma de hacer esto es determinando la bounding box para cada ráster e intersectándolas empleando el módulo pyhon shapely. La bounding box de intersección es usada para determinar los índices de columnna, fila para el primero y el último píxel de cada ráster. Por tanto, estos valores son usados para la lectura por ventanas de cada ráster como arrays antes de calcular la suma respectiva. Finalmente, este array es escrito como ráster empleando el método 'write' de rasterio con las opciones de configuración adecuadas. El código desarrollado se presenta a continuación:

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66

import rasterio from rasterio.windows import Window from rasterio.transform import Affine from shapely.geometry import box raster1 = rasterio.open("/home/zeito/pyqgis_data/c1.tif") raster2 = rasterio.open("/home/zeito/pyqgis_data/c2.tif") bb_raster1 = box(raster1.bounds[0], raster1.bounds[1], raster1.bounds[2], raster1.bounds[3]) bb_raster2 = box(raster2.bounds[0], raster2.bounds[1], raster2.bounds[2], raster2.bounds[3]) xminR1, yminR1, xmaxR1, ymaxR1 = raster1.bounds xminR2, yminR2, xmaxR2, ymaxR2 = raster2.bounds intersection = bb_raster1.intersection(bb_raster2) transform = Affine(raster1.res[0], 0.0, intersection.bounds[0], 0.0, -raster1.res[1], intersection.bounds[3]) p1Y = intersection.bounds[3] - raster1.res[1]/2 p1X = intersection.bounds[0] + raster1.res[0]/2 p2Y = intersection.bounds[1] + raster1.res[1]/2 p2X = intersection.bounds[2] - raster1.res[0]/2 #row index raster1 row1R1 = int((ymaxR1 - p1Y)/raster1.res[1]) #row index raster2 row1R2 = int((ymaxR2 - p1Y)/raster2.res[1]) #column index raster1 col1R1 = int((p1X - xminR1)/raster1.res[0]) #column index raster2 col1R2 = int((p1X - xminR2)/raster1.res[0]) #row index raster1 row2R1 = int((ymaxR1 - p2Y)/raster1.res[1]) #row index raster2 row2R2 = int((ymaxR2 - p2Y)/raster2.res[1]) #column index raster1 col2R1 = int((p2X - xminR1)/raster1.res[0]) #column index raster2 col2R2 = int((p2X - xminR2)/raster1.res[0]) width1 = col2R1 - col1R1 + 1 width2 = col2R2 - col1R2 + 1 height1 = row2R1 - row1R1 + 1 height2 = row2R2 - row1R2 + 1 arr_raster1 = raster1.read(1, window=Window(col1R1, row1R1, width1, height1)) arr_raster2 = raster2.read(1, window=Window(col1R2, row1R2, width2, height2)) arr_sum = arr_raster1 + arr_raster2 # Register GDAL format drivers and configuration options with a # context manager. with rasterio.open('/home/zeito/pyqgis_data/sum3.tif', 'w', driver='GTiff', height=arr_sum.shape[0], width=arr_sum.shape[1], count=1, dtype=arr_sum.dtype, crs=raster1.crs, transform=transform) as dst: dst.write(arr_sum, 1) dst.close()

El código anterior se probó en PyCharm empleando los dos rásteres relativamente grandes que se visualizan en la imagen siguiente:

Después de ejecutar el script, el ráster resultante se cargó también en QGIS y, con la ayuda del plugin Value Tool, se corroboró en mucho puntos (incluyendo el primero y último píxel) que la suma de rásteres se efectuó tal como se esperaba (ver imagen siguiente).