Skip to content

Commit

Permalink
only forward epsg if present
Browse files Browse the repository at this point in the history
  • Loading branch information
vincentsarago committed Oct 29, 2024
1 parent b4c7285 commit 1fdfabe
Show file tree
Hide file tree
Showing 3 changed files with 33 additions and 25 deletions.
40 changes: 17 additions & 23 deletions rio_stac/stac.py
Original file line number Diff line number Diff line change
Expand Up @@ -43,6 +43,7 @@ def get_dataset_geom(
src_dst: Union[DatasetReader, DatasetWriter, WarpedVRT, MemoryFile],
densify_pts: int = 0,
precision: int = -1,
geographic_crs: rasterio.crs.CRS = EPSG_4326,
) -> Dict:
"""Get Raster Footprint."""
if densify_pts < 0:
Expand All @@ -53,7 +54,7 @@ def get_dataset_geom(
geom = bbox_to_geom(src_dst.bounds)

# 2. Densify the Polygon geometry
if src_dst.crs != EPSG_4326 and densify_pts:
if src_dst.crs != geographic_crs and densify_pts:
# Derived from code found at
# https://stackoverflow.com/questions/64995977/generating-equidistance-points-along-the-boundary-of-a-polygon-but-cw-ccw
coordinates = numpy.asarray(geom["coordinates"][0])
Expand All @@ -69,7 +70,7 @@ def get_dataset_geom(
}

# 3. Reproject the geometry to "epsg:4326"
geom = warp.transform_geom(src_dst.crs, EPSG_4326, geom, precision=precision)
geom = warp.transform_geom(src_dst.crs, geographic_crs, geom, precision=precision)
bbox = feature_bounds(geom)

else:
Expand Down Expand Up @@ -99,27 +100,12 @@ def get_projection_info(
see: https://github.com/stac-extensions/projection
"""
projjson = None
wkt2 = None

epsg = None
if src_dst.crs is not None:
# EPSG
epsg = src_dst.crs.to_epsg() if src_dst.crs.is_epsg_code else None

# PROJJSON
try:
projjson = src_dst.crs.to_dict(projjson=True)
except (AttributeError, TypeError) as ex:
warnings.warn(f"Could not get PROJJSON from dataset : {ex}")
pass

# WKT2
try:
wkt2 = src_dst.crs.to_wkt()
except Exception as ex:
warnings.warn(f"Could not get WKT2 from dataset : {ex}")
pass

meta = {
"epsg": epsg,
"geometry": bbox_to_geom(src_dst.bounds),
Expand All @@ -128,11 +114,17 @@ def get_projection_info(
"transform": list(src_dst.transform),
}

if projjson is not None:
meta["projjson"] = projjson

if wkt2 is not None:
meta["wkt2"] = wkt2
if not epsg and src_dst.crs:
# WKT2
try:
meta["wkt2"] = src_dst.crs.to_wkt()
except Exception as ex:
warnings.warn(f"Could not get WKT2 from dataset : {ex}")
# PROJJSON
try:
meta["projjson"] = src_dst.crs.to_dict(projjson=True)
except (AttributeError, TypeError) as ex:
warnings.warn(f"Could not get PROJJSON from dataset : {ex}")

return meta

Expand Down Expand Up @@ -300,6 +292,7 @@ def create_stac_item(
raster_max_size: int = 1024,
geom_densify_pts: int = 0,
geom_precision: int = -1,
geographic_crs: rasterio.crs.CRS = EPSG_4326,
) -> pystac.Item:
"""Create a Stac Item.
Expand Down Expand Up @@ -352,6 +345,7 @@ def create_stac_item(
src_dst,
densify_pts=geom_densify_pts,
precision=geom_precision,
geographic_crs=geographic_crs,
)

media_type = (
Expand Down
Binary file added tests/fixtures/dataset_mars.tif
Binary file not shown.
18 changes: 16 additions & 2 deletions tests/test_create_item.py
Original file line number Diff line number Diff line change
Expand Up @@ -115,8 +115,8 @@ def test_create_item_options():
]
assert "datetime" in item_dict["properties"]
assert "proj:epsg" in item_dict["properties"]
assert "proj:wkt2" in item_dict["properties"]
assert "proj:projjson" in item_dict["properties"]
assert "proj:wkt2" not in item_dict["properties"]
assert "proj:projjson" not in item_dict["properties"]
assert "sci:citation" in item_dict["properties"]

# external assets
Expand Down Expand Up @@ -330,3 +330,17 @@ def test_densify_geom():
item_dens_dict = item_dens.to_dict()

assert item_dict["bbox"] != item_dens_dict["bbox"]


def test_mars_dataset():
"""Test with Mars Dataset."""
MARS2000_SPHERE = rasterio.crs.CRS.from_proj4("+proj=longlat +R=3396190 +no_defs")
src_path = os.path.join(PREFIX, "dataset_mars.tif")

item = create_stac_item(src_path, geographic_crs=MARS2000_SPHERE, with_proj=True)
assert item.validate()
item_dict = item.to_dict()

assert not item_dict["properties"].get("proj:epsg")
assert "proj:projjson" not in item_dict["properties"]
assert "proj:wkt2" in item_dict["properties"]

0 comments on commit 1fdfabe

Please sign in to comment.