diff --git a/rio_stac/stac.py b/rio_stac/stac.py index b8b4f21..51a1387 100644 --- a/rio_stac/stac.py +++ b/rio_stac/stac.py @@ -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: @@ -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]) @@ -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: @@ -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), @@ -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 @@ -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. @@ -352,6 +345,7 @@ def create_stac_item( src_dst, densify_pts=geom_densify_pts, precision=geom_precision, + geographic_crs=geographic_crs, ) media_type = ( diff --git a/tests/fixtures/dataset_mars.tif b/tests/fixtures/dataset_mars.tif new file mode 100644 index 0000000..1c624f6 Binary files /dev/null and b/tests/fixtures/dataset_mars.tif differ diff --git a/tests/test_create_item.py b/tests/test_create_item.py index a62fad9..bd99d78 100644 --- a/tests/test_create_item.py +++ b/tests/test_create_item.py @@ -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 @@ -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"]