handle 2D, no CRS, and read CRS from S3

Author: sclaw

external_caller.py

@@ -1,15 +1,32 @@
+import json
+import logging
 import sys
 import warnings
 
+from botocore.exceptions import NoSuchKey
+
 from ras_stac.ras1d.converter import process_in_place_s3
+from ras_stac.ras1d.utils.s3_utils import str_from_s3
 
 warnings.filterwarnings("ignore")
 
 from papipyplug import parse_input, plugin_logger, print_results
 
 PLUGIN_PARAMS = {"required": ["in_prefix", "crs", "out_prefix"], "optional": []}
 
-if __name__ == "__main__":
+
+def owp_wrapper(in_prefix: str, out_prefix: str) -> dict:
+    """Temporary wrapper to meet OWP MIP 30% deliverable deadline"""
+    crs_path = in_prefix.replace("source_models", "source_crs") + "crs_inference.json"
+    try:
+        crs_dict = json.loads(str_from_s3(crs_path))
+        crs = crs_dict["best_crs"]
+    except NoSuchKey:
+        crs = None
+    return process_in_place_s3(in_prefix, crs, out_prefix)
+
+
+def main():
     plugin_logger()
 
     input_params = parse_input(sys.argv, PLUGIN_PARAMS)
@@ -18,8 +35,19 @@
     crs = input_params.get("crs")
     out_prefix = input_params.get("out_prefix")
 
+    # Debugging option
     if in_prefix == "TEST":
         results = f"{crs} | {out_prefix}"
-    else:
-        results = process_in_place_s3(in_prefix, crs, out_prefix)
-    print_results(results)
+        print_results(results)
+        return
+
+    # Process
+    try:
+        results = owp_wrapper(in_prefix, out_prefix)
+        print_results(results)
+    except Exception as e:
+        logging.warning(e)
+
+
+if __name__ == "__main__":
+    main()

ras_stac/ras1d/converter.py

@@ -1,7 +1,8 @@
-import datetime
 import io
 import json
+import logging
 import sys
+from datetime import datetime
 from pathlib import Path
 
 import pandas as pd
@@ -14,6 +15,7 @@
 from ras_stac.ras1d.utils.classes import (
     GenericAsset,
     GeometryAsset,
+    NullGeometryAsset,
     PlanAsset,
     SteadyFlowAsset,
     ThumbAsset,
@@ -48,6 +50,8 @@ def export_stac(self, output_path: str) -> None:
     def export_thumbnail(self, thumb_path: str) -> None:
         """Generate STAC thumbnail, save to S3, and log path."""
         gdfs = self.primary_geometry.gdfs
+        if len(gdfs) == 0 or "null" in gdfs:
+            return
         thumb = make_thumbnail(gdfs)
         if file_location(thumb_path) == "local":
             thumb.savefig(thumb_path, dpi=80)
@@ -72,15 +76,16 @@ def stac_item(self) -> dict:
         )
         stor_ext = StorageExtension.ext(stac, add_if_missing=True)
         stor_ext.apply(platform="AWS", region="us-east-1")
-        prj_ext = AssetProjectionExtension.ext(stac, add_if_missing=True)
-        og_crs = CRS(self.crs)
-        prj_ext.apply(
-            epsg=og_crs.to_epsg(),
-            wkt2=og_crs.to_wkt(),
-            geometry=self.get_footprint(),
-            bbox=self.get_bbox(),
-            centroid=to_geojson(self.get_centroid()),
-        )
+        if self.crs:
+            prj_ext = AssetProjectionExtension.ext(stac, add_if_missing=True)
+            og_crs = CRS(self.crs)
+            prj_ext.apply(
+                epsg=og_crs.to_epsg(),
+                wkt2=og_crs.to_wkt(),
+                geometry=self.get_footprint(),
+                bbox=self.get_bbox(),
+                centroid=to_geojson(self.get_centroid()),
+            )
         return stac
 
     @property
@@ -139,6 +144,7 @@ def stac_properties(self):
             "datetime_source": "processing_time" if self.primary_geometry.last_update is None else "model_geometry",
             "assigned_HUC8": self.huc8,
             "has_2d": any([a.has_2d for a in self.assets if isinstance(a, GeometryAsset)]),
+            "has_1d": any([a.has_1d for a in self.assets if isinstance(a, GeometryAsset)]),
         }
         for p in self.custom_properties:
             properties[p] = self.custom_properties[p]
@@ -177,7 +183,16 @@ def primary_plan(self) -> PlanAsset:
     @property
     def primary_geometry(self) -> GeometryAsset:
         """The geometry file listed in the primary plan"""
-        return self.extension_dict[self.primary_plan.geometry]
+        if not self.crs:
+            return NullGeometryAsset()
+        try:
+            geom = self.extension_dict[self.primary_plan.geometry]
+        except Exception:
+            return NullGeometryAsset()
+        if not geom.has_1d:
+            return NullGeometryAsset()
+        else:
+            return geom
 
     def check_for_mip(self) -> None:
         mip_data = [a for a in self.assets if a.name == "mip_package_geolocation_metadata.json"]
@@ -214,18 +229,24 @@ def ras_to_stac(ras_dir: str, crs: str):
 
 def process_in_place_s3(in_prefix: str, crs: str, out_prefix: str):
     """Convert a HEC-RAS model to a STAC item and save to same directory."""
+    logging.info(f"Processing model with crs {crs} at prefix {in_prefix}")
+    logging.info("Discovering model contents")
     converter = from_directory(in_prefix, crs)
     converter.check_for_mip()
     thumb_path = out_prefix + "Thumbnail.png"
+    logging.info(f"Generating thumbnail at {thumb_path}")
     converter.export_thumbnail(thumb_path)
     stac_path = out_prefix + f"{converter.idx}.json"
+    logging.info(f"Generating STAC item at {thumb_path}")
     converter.export_stac(stac_path)
     return {"in_path": in_prefix, "crs": crs, "thumb_path": thumb_path, "stac_path": stac_path}
 
 
 if __name__ == "__main__":
     ras_dir = sys.argv[1]
     crs = sys.argv[2]
+    if crs == "None":
+        crs = None
     out_dir = sys.argv[3]
-    process_in_place_s3(ras_dir, crs, out_dir)
-    # ras_to_stac(ras_dir, crs)
+    # process_in_place_s3(ras_dir, crs, out_dir)
+    ras_to_stac(ras_dir, crs)

ras_stac/ras1d/utils/classes.py

@@ -10,13 +10,9 @@
 import pystac
 from pyproj import CRS
 from shapely import make_valid, union_all
-from shapely.geometry import (
-    LineString,
-    MultiPolygon,
-    Point,
-    Polygon,
-)
+from shapely.geometry import LineString, MultiPolygon, Point, Polygon, shape
 
+from ras_stac.ras1d.data.us_geom import us_bounds
 from ras_stac.ras1d.utils.common import file_location
 from ras_stac.ras1d.utils.ras_utils import (
     data_pairs_from_text_block,
@@ -47,7 +43,7 @@ def check_crs(func):
 
     def wrapper(self, *args, **kwargs):
         if self.crs is None:
-            raise ValueError("Projection cannot be None")
+            return None
         return func(self, *args, **kwargs)
 
     return wrapper
@@ -221,6 +217,24 @@ def _extra_fields(self):
         return ex
 
 
+class NullGeometryAsset(GenericAsset):
+
+    def __init__(self):
+        features = us_bounds["features"]
+        geometries = [shape(feature["geometry"]) for feature in features]
+        properties = [feature["properties"] for feature in features]
+        gdf = gpd.GeoDataFrame(properties, geometry=geometries, crs="epsg:4326")
+        self.gdfs = {"null": gdf}
+        self.concave_hull = gdf
+
+    def __getattr__(self, name):
+        # Return None if the attribute is not found
+        return None
+
+    def get_river_miles(self):
+        return None
+
+
 class GeometryAsset(GenericAsset):
 
     def __init__(self, url: str):
@@ -253,7 +267,7 @@ def version(self):
         return search_contents(self.contents, "Program Version", expect_one=False)
 
     @property
-    @check_crs
+    # @check_crs
     def reaches(self) -> dict:
         """A dictionary of the reaches contained in the HEC-RAS geometry file."""
         river_reaches = search_contents(self.contents, "River Reach", expect_one=False)
@@ -263,7 +277,7 @@ def reaches(self) -> dict:
         return reaches
 
     @property
-    @check_crs
+    # @check_crs
     def rivers(self) -> dict:
         """A nested river-reach dictionary of the rivers/reaches contained in the HEC-RAS geometry file."""
         rivers = {}
@@ -370,7 +384,7 @@ def n_junctions(self):
         return len(self.junctions)
 
     @property
-    @check_crs
+    # @check_crs
     def n_rivers(self):
         """Number of rivers in the HEC-RAS geometry file."""
         return len(self.rivers)
@@ -505,6 +519,11 @@ def has_2d(self):
                 return True
         return False
 
+    @property
+    def has_1d(self):
+        """Check if RAS geometry has any 1D components"""
+        return self.n_rivers > 0
+
 
 class XS:
     """HEC-RAS Cross Section."""