Skip to content

Instantly share code, notes, and snippets.

@perrygeo
Created January 12, 2020 20:16
Show Gist options
  • Save perrygeo/0781c4724d3c5b4b1afbfba2fcace811 to your computer and use it in GitHub Desktop.
Save perrygeo/0781c4724d3c5b4b1afbfba2fcace811 to your computer and use it in GitHub Desktop.
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))
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment