Creación de una trama de porcentaje de píxeles ocupados a partir de un shapefile

Question

Partiendo de un shapefile poligonal, ¿cómo puedo producir un raster de 1 km, donde los valores de las celdas son el porcentaje del píxel ocupado por el shapefile? Por ejemplo, los valores irían de 0 a 100%. Me gustaría utilizar herramientas de python disponibles gratuitamente como osgeo, shapely, rasterstats, etc. para producir el raster de salida y tiene que ser una solución eficiente porque se trata de un cálculo global de 1 km.

Answer 1

Aquí hay una función que escribí que hace exactamente esto, excepto que la cuadrícula está definida por una trama de entrada y mi polígono no era un MultiPolígono, y no tenía agujeros. En tu caso, "geom" sería tu shapefile azul, y necesitas especificar los xs y los ys de tu cuadrícula de 1km. Actualmente, esos valores se leen del objeto raster subyacente (fsrc). Al final del código, tendrás que reemplazar gt[0]*gt[4] por la resolución de tu celda. Creo que el objeto "geom" es una geometría de forma, pero no estoy seguro. Es lo que rasterstats devuelve en su cuerpo de código principal en la siguiente línea: geom = shape(feat['geometry']) .

Tenga en cuenta que es mejor que su cuadrícula y su polígono estén en la misma proyección.

def fractional_pixel_weights(fsrc, geom):

""" 
Returns a grid of the same size as fsrc, where each value represents the
fraction of the cell that is filled by the polygon in geom.

fsrc is a rasterstats-created object (I think rasterio).
"""

gt = fsrc.affine
xs = np.arange(gt[2], gt[2] +  gt[0]* (1 + fsrc.shape[1]), gt[0])
ys = np.arange(gt[5], gt[5] +  gt[4]* (1 + fsrc.shape[0]), gt[4])

# Convert geom into ogr geometry
geom_ogr = ogr.CreateGeometryFromWkt(geom.to_wkt())

# Loop through each grid cell, compute the intersecting area
overlapping_areas = np.empty((len(ys)-1, len(xs)-1))
for ix in range(len(xs)-1):
    xmin = xs[ix]
    xmax = xs[ix + 1]
    for iy in range(len(ys)-1):
        ymax = ys[iy]
        ymin = ys[iy + 1]

        # Intersecting area
        coords_wkt = "POLYGON ((" + str(xmin) + ' ' + str(ymax) + ', ' + str(xmax) + ' ' + str(ymax) + ', ' + str(xmax) + ' ' + str(ymin) + ', ' + str(xmin) + ' ' + str(ymin) + ', ' + str(xmin) + ' ' + str(ymax) + "))"
        polycell = ogr.CreateGeometryFromWkt(coords_wkt)

        overlapping_areas[iy, ix] = polycell.Intersection(geom_ogr).Area()

# Ratio of overlapped area to pixel area
frac_intersected = overlapping_areas / (abs(gt[0] * gt[4]))

return frac_intersected

Esta no es una respuesta completa porque, como se ha dicho, utiliza una trama subyacente para crear la cuadrícula, pero espero que sea un comienzo lo suficientemente bueno para que puedas implementar lo que quieres.

También me acabo de dar cuenta de que tienes lo que parece un shapefile MultiPolígono. Yo seguiría el consejo del otro contestador y lo disolvería primero en uno.

Answer 2

Depende de lo que quieras optimizar (velocidad, cpu, memoria, etc)

Una opción sería

1) rasterizar sus polígonos - las celdas completamente dentro/fuera de los polígonos se establecen con precisión, las celdas en el borde no

2) disolver los polígonos para eliminar los límites entre polígonos (no tiene que ser un solo polígono, sino los menos posibles)

3) en el código, crear un polígono cuadrado que tenga el mismo tamaño que una celda raster.

4) Haga un bucle a través de cada uno de sus polígonos disueltos, mueva el polígono cuadrado a lo largo del borde del polígono (alineado en la coordenada de las celdas) y calcule la superposición.