dask-geomodeling AggregateRaster example

aggregate_raster.py

### INSTALLATION
# follow https://dask-geomodeling.readthedocs.io/en/docs/installation.html#local-setup-for-development-ubuntu
# additionally, constrain geopandas to 0.6  (pip install geopandas==0.6)

from dask_geomodeling.geometry import GeometryFileSource, AggregateRaster, GeometryTiler
from dask_geomodeling.raster import RasterFileSource
from shapely.geometry import box

# construct the View
rs = RasterFileSource("path/to/tif")
gs = GeometryFileSource(
    "path/to/shp",
    id_field="FID"  # specify ID field if it is non-standard like this
)
agg = AggregateRaster(gs, rs, statistic="max")

# Build the data request using the complete raster extent in its native projection
x1, x2, y1, y2 = rs.geometry.GetEnvelope()
projection = rs.projection
request = dict(geometry=box(x1, y1, x2, y2), projection=projection)

# get the data directly
result = agg.get_data(**request)
print("Feature count:", len(result["features"]))

# tile it, get a graph and compute it remotely
tiled = GeometryTiler(agg, size=500, projection=projection)
graph, name = tiled.get_compute_graph(**request)
# -> see dask docs

# output to a file with tiling (convenience function)
agg.to_file("output2.shp", tile_size=500, **request)