hexbin_population.py

from h3 import h3
import json

from rasterstats import gen_zonal_stats


def hexbin_features(bounding_polygon, zoom=8):
    """Takes a GeoJSON-like geometry dictionary with type Polygon
    and yields GeoJSON-like Features representing the hexbins
    that cover the polygon at a given zoom level
    """
    assert bounding_polygon["type"] == "Polygon"

    # Reverse coordinate order, h3 expects lat, lng on input
    bounding_coords = bounding_polygon["coordinates"][0]
    bounding_geom_latlng = {
        "type": "Polygon",
        "coordinates": [[(lat, lng) for lng, lat in bounding_coords]],
    }

    for hexbin in h3.polyfill(bounding_geom_latlng, zoom):
        coords = h3.h3_set_to_multi_polygon([hexbin], geo_json=True)
        yield {
            "type": "Feature",
            "properties": {"hexid": hexbin},
            "geometry": {"type": "Polygon", "coordinates": coords[0]},
        }


def bbox_to_polygon(w, s, e, n):
    return {
        "type": "Polygon",
        "coordinates": [[[w, s], [w, n], [e, n], [e, s], [w, s]]],
    }


if __name__ == "__main__":

    # Raster data is the 30 arc-second gridded population of the world, v4
    raster_path = "/Users/mperry/Downloads/gpw-v4-population-count-rev11_2020_30_sec_tif/gpw_v4_population_count_rev11_2020_30_sec.tif"

    # Vector data is a generator of hexbins as GeoJSON-like features
    # For state of Colorado at hexbin level 7
    co = (-109.060253, 36.992426, -102.041524, 41.003444)
    zoom = 5
    poly = bbox_to_polygon(*co)
    hexbin_generator = hexbin_features(poly, zoom)

    # Outputs hexbin feautres with additional properties: "population_sum": <float>
    for feature in gen_zonal_stats(
        hexbin_generator,
        raster_path,
        stats="sum std",
        prefix="population_",
        geojson_out=True,
    ):
        print(json.dumps(feature))