Efficient AI4EO OpenSource framework
You can not select more than 25 topics Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.

83 lines
3.0KB

  1. from rasterio.crs import CRS
  2. from rasterio.warp import transform_geom
  3. from rasterio.features import rasterize
  4. from rasterio.transform import from_bounds
  5. import mercantile
  6. from supermercado import burntiles
  7. from shapely.geometry import shape, mapping
  8. from neat_eo.tiles import tile_bbox
  9. def geojson_parse_feature(zoom, srid, feature_map, feature, buffer=0):
  10. def geojson_parse_polygon(zoom, srid, feature_map, polygon):
  11. if isinstance(polygon["coordinates"], list): # https://github.com/Toblerity/Shapely/issues/245
  12. for i, ring in enumerate(polygon["coordinates"]): # GeoJSON coordinates could be N dimensionals
  13. polygon["coordinates"][i] = [[x, y] for point in ring for x, y in zip([point[0]], [point[1]])]
  14. if srid != 4326:
  15. try:
  16. polygon = transform_geom(CRS.from_epsg(srid), CRS.from_epsg(4326), polygon)
  17. except: # negative buffer could lead to empty/invalid geom
  18. return feature_map
  19. try:
  20. for tile in burntiles.burn([{"type": "feature", "geometry": polygon}], zoom=zoom):
  21. feature_map[mercantile.Tile(*tile)].append({"type": "feature", "geometry": polygon})
  22. except:
  23. pass
  24. return feature_map
  25. def geojson_parse_geometry(zoom, srid, feature_map, geometry, buffer):
  26. if buffer:
  27. geometry = transform_geom(CRS.from_epsg(srid), CRS.from_epsg(3857), geometry) # be sure to be planar
  28. geometry = mapping(shape(geometry).buffer(buffer))
  29. srid = 3857
  30. if geometry["type"] == "Polygon":
  31. feature_map = geojson_parse_polygon(zoom, srid, feature_map, geometry)
  32. elif geometry["type"] == "MultiPolygon":
  33. for polygon in geometry["coordinates"]:
  34. feature_map = geojson_parse_polygon(zoom, srid, feature_map, {"type": "Polygon", "coordinates": polygon})
  35. return feature_map
  36. if not feature or not feature["geometry"]:
  37. return feature_map
  38. if feature["geometry"]["type"] == "GeometryCollection":
  39. for geometry in feature["geometry"]["geometries"]:
  40. feature_map = geojson_parse_geometry(zoom, srid, feature_map, geometry, buffer)
  41. else:
  42. feature_map = geojson_parse_geometry(zoom, srid, feature_map, feature["geometry"], buffer)
  43. return feature_map
  44. def geojson_srid(feature_collection):
  45. try:
  46. crs_mapping = {"CRS84": "4326", "900913": "3857"}
  47. srid = feature_collection["crs"]["properties"]["name"].split(":")[-1]
  48. srid = int(srid) if srid not in crs_mapping else int(crs_mapping[srid])
  49. except:
  50. srid = int(4326)
  51. return srid
  52. def geojson_tile_burn(tile, features, srid, ts, burn_value=1):
  53. """Burn tile with GeoJSON features."""
  54. crs = (CRS.from_epsg(srid), CRS.from_epsg(3857))
  55. shapes = ((transform_geom(*crs, feature["geometry"]), burn_value) for feature in features)
  56. try:
  57. return rasterize(shapes, out_shape=ts, transform=from_bounds(*tile_bbox(tile, mercator=True), *ts))
  58. except:
  59. return None