Merge pull request #60 from developmentseed/feature/only-epsg

Feature/only epsg
developmentseed · Oct 29, 2024 · cff3fcc · cff3fcc
2 parents b4c7285 + a3915ae
commit cff3fcc
Show file tree

Hide file tree

Showing 5 changed files with 38 additions and 26 deletions.
diff --git a/CHANGES.md b/CHANGES.md
@@ -2,6 +2,8 @@
 ## 0.10.0 (2024-10-29)
 
 * handle `TIFFTAG_DATETIME` metadata for STAC datetime
+* only set `proj:epsg` if present else use `proj:wkt2` or `proj:projjson`
+* add `geographic_crs` parameter to `create_stac_item` function to enable non-earth raster dataset
 
 ## 0.9.0 (2023-12-08)
 

diff --git 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
diff --git a/tests/test_cli.py b/tests/test_cli.py
@@ -129,6 +129,6 @@ def test_rio_stac_cli(runner):
             ]
             assert "datetime" in stac_item["properties"]
             assert "proj:epsg" in stac_item["properties"]
-            assert "proj:projjson" in stac_item["properties"]
+            assert "proj:projjson" not in stac_item["properties"]
             assert "raster:bands" in stac_item["assets"]["asset"]
             assert "eo:bands" in stac_item["assets"]["asset"]
diff --git 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,19 @@ 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" in item_dict["properties"]
+        or "proj:wkt2" in item_dict["properties"]
+    )