Use Projected CRS for Mask Buffering (Instead of Silencing geopandas Warning)

NextGen_Forcings_Engine_BMI/NextGen_Forcings_Engine/core/geoMod.py

@@ -250,6 +250,14 @@ def y_coords(self) -> np.ndarray:
                     y_coords[:, :] = np.flipud(y_coords[:, :])
             return y_coords
 
+    @property
+    def approx_centroid_global_xy(self) -> tuple[float, float]:
+        """The approximate centroid in global coordinates (a tuple of 2 floats: (x, y))"""
+        x_mean = np.mean([c[0] for c in self.elementcoords_global])
+        y_mean = np.mean([c[1] for c in self.elementcoords_global])
+        LOG.debug(f"Approximate centroid: ({x_mean}, {y_mean})")
+        return (x_mean, y_mean)
+
 
 class GriddedGeoMeta(GeoMeta):
     """Class for handling information about the gridded domains for forcing."""

NextGen_Forcings_Engine_BMI/NextGen_Forcings_Engine/general_utils.py

@@ -7,6 +7,7 @@
 from collections import OrderedDict
 
 import numpy as np
+from pyproj import CRS
 
 JSON_NOT_SERIALIZABLE_SENTINEL = "ERR_NOT_JSON_SERIALIZABLE"
 JSON_NOT_SERIALIZABLE_FORMAT = JSON_NOT_SERIALIZABLE_SENTINEL + ":TYPE:{typ}"
@@ -239,3 +240,13 @@ def rand_str(length: int) -> str:
             f"length requested was {length}, but this function only supports length 1 through 32"
         )
     return str(uuid.uuid4()).replace("-", "")[:length]
+
+
+def crs_assert_projected_horizontal_meters(crs: CRS) -> None:
+    """Assert that the CRS is projected and has horizontal units of meters."""
+    if not crs.is_projected:
+        raise ValueError(f"CRS is not projected: {crs}")
+    if crs.axis_info[0].unit_conversion_factor != 1:
+        raise ValueError(
+            f"Expected crs.axis_info[0].unit_conversion_factor == 1, but got: {crs.axis_info[0].unit_conversion_factor}"
+        )

NextGen_Forcings_Engine_BMI/NextGen_Forcings_Engine/historical_forcing.py

@@ -28,9 +28,11 @@
     ConfigOptions,
 )
 from NextGen_Forcings_Engine_BMI.NextGen_Forcings_Engine.core.parallel import MpiConfig
-from NextGen_Forcings_Engine_BMI.NextGen_Forcings_Engine.general_utils import rand_str
+from NextGen_Forcings_Engine_BMI.NextGen_Forcings_Engine.general_utils import (
+    crs_assert_projected_horizontal_meters,
+    rand_str,
+)
 
-warnings.filterwarnings("ignore", module="geopandas")
 LOG = logging.getLogger("FORCING")
 
 zarr.config.set({"async.concurrency": 100})
@@ -49,17 +51,19 @@ def __init__(
         self.config_options = config_options
         self.mpi_config = mpi_config
         self.wrf_hydro_geo_meta = wrf_hydro_geo_meta
-        self.dest_crs = CRS(4326)
-        self.buffer = 0.02  # degree buffer around bounding box
 
     @cached_property
     def bounds(self) -> tuple[float, float, float, float]:
         """Get bounding box from geospatial dataframe.
 
-        Apply buffer in known crs/units (degrees) and then convert back to src_crs.
+        Apply buffer in known crs/units (m) and then convert back to src_crs.
         """
+        LOG.debug(
+            f"Temporary CRS for creating a mask (will buffer AOI by {self.buffer} in this CRS): {self._temp_crs}"
+        )
+        crs_assert_projected_horizontal_meters(self._temp_crs)
         return (
-            self.gdf.to_crs(self.dest_crs)
+            self.gdf.to_crs(self._temp_crs)
             .buffer(self.buffer)
             .to_crs(self.src_crs)
             .total_bounds
@@ -396,6 +400,8 @@ def __init__(
         self.x_label = "longitude"
         self.y_label = "latitude"
         self.time_label = "time"
+        self.buffer = 6000  # m buffer around bounding box. Use 6km buffer in case someone applies this to legacy 4km AORC data instead of the newer 1km AORC data.
+        self._temp_crs = CRS(5070)
 
     @cached_property
     def src_crs(self) -> CRS:
@@ -466,6 +472,8 @@ def __init__(
         self.x_label = "longitude"
         self.y_label = "latitude"
         self.time_label = "time"
+        self.buffer = 6000  # m buffer around bounding box. Use 6km buffer in case someone applies this to legacy 4km AORC data instead of the newer 1km AORC data.
+        self._temp_crs = CRS(3338)
 
     @cached_property
     def src_crs(self):
@@ -536,6 +544,7 @@ def __init__(
         self.x_label = "x"
         self.y_label = "y"
         self.time_label = "time"
+        self.buffer = 6000  # m buffer around bounding box
 
     @property
     def vars(
@@ -565,6 +574,7 @@ def __init__(
     ):
         """Initialize NWM CONUS processor."""
         super().__init__(config_options, mpi_config, wrf_hydro_geo_meta)
+        self._temp_crs = CRS(5070)
 
     def url(self, var: str) -> str:
         """Generate NWM S3 zarr URL for current variable.
@@ -685,17 +695,54 @@ def s3_lazy_ds(self) -> xr.Dataset:
         return xr.open_zarr(ObjectStore(object_store))
 
 
+class NWMV3PuertoRicoProcessor(NWMV3OConusProcessor):
+    """Processor for NWM Puerto Rico data."""
+
+    def __init__(
+        self,
+        config_options: ConfigOptions,
+        mpi_config: MpiConfig,
+        wrf_hydro_geo_meta: dict,
+    ):
+        """Initialize NWM Puerto Rico processor."""
+        super().__init__(config_options, mpi_config, wrf_hydro_geo_meta)
+        self._temp_crs = CRS(32161)
+
+
+class NWMV3HawaiiProcessor(NWMV3OConusProcessor):
+    """Processor for NWM Hawaii data."""
+
+    def __init__(
+        self,
+        config_options: ConfigOptions,
+        mpi_config: MpiConfig,
+        wrf_hydro_geo_meta: dict,
+    ):
+        """Initialize NWM Hawaii processor."""
+        super().__init__(config_options, mpi_config, wrf_hydro_geo_meta)
+        lon, lat = wrf_hydro_geo_meta.approx_centroid_global_xy
+        if not -180 < lon < 180:
+            raise ValueError(f"Unexpected (lon, lat) = ({lon}, {lat})")
+        utm_zone_number = int((lon + 180) / 6) + 1
+        if utm_zone_number not in (1, 2, 3, 4, 5):
+            raise ValueError(
+                f"Unexpected UTM zone {utm_zone_number} for Hawaii. Expected zone 1 through 5. (lon, lat) = ({lon}, {lat})"
+            )
+        self._temp_crs = CRS(f"EPSG:3260{utm_zone_number}")
+
+
 class NWMV3AlaskaProcessor(NWMV3Processor):
-    """Processor for NWM OCONUS data."""
+    """Processor for NWM Alaska data."""
 
     def __init__(
         self,
         config_options: ConfigOptions,
         mpi_config: MpiConfig,
         wrf_hydro_geo_meta: dict,
     ):
-        """Initialize NWM OCONUS processor."""
+        """Initialize NWM Alaska processor."""
         super().__init__(config_options, mpi_config, wrf_hydro_geo_meta)
+        self._temp_crs = CRS(3338)
 
     @cached_property
     def url(self) -> str:

NextGen_Forcings_Engine_BMI/NextGen_Forcings_Engine/model.py

@@ -1,8 +1,9 @@
 import datetime
+import logging
 import os
 from contextlib import contextmanager
 from time import time
-import logging
+
 import numpy as np
 import pandas as pd
 from ewts import Payload as Pld
@@ -29,7 +30,8 @@
     AORCConusProcessor,
     NWMV3AlaskaProcessor,
     NWMV3ConusProcessor,
-    NWMV3OConusProcessor,
+    NWMV3HawaiiProcessor,
+    NWMV3PuertoRicoProcessor,
 )
 
 LOG = logging.getLogger("FORCING")
@@ -227,8 +229,8 @@ def determine_forecast(
             # If we're in an AnA configuration, then must offset the BMI future
             # timestamp to account for the "lookback" period being properly iterated
             # over between 3-28 hour look back time period and operation configuration
-            
-            #if config_options.input_forcings[0] in [20, 22]:
+
+            # if config_options.input_forcings[0] in [20, 22]:
             #     config_options.current_fcst_cycle = (
             #         config_options.b_date_proc
             #         + pd.TimedeltaIndex(
@@ -243,19 +245,19 @@ def determine_forecast(
             #     )
             #     config_options.future_time = future_time
             # else:
-             
+
             # Puerto Rico / Hawaii AnA: 1-hour lookback (based on 6-hourly forecast cycles)
             config_options.current_fcst_cycle = (
                 config_options.b_date_proc
-                + pd.TimedeltaIndex(
-                    np.array([future_time - 3600.0], dtype=float), "s"
-                )[0]
+                + pd.TimedeltaIndex(np.array([future_time - 3600.0], dtype=float), "s")[
+                    0
+                ]
             )
             config_options.current_time = (
                 config_options.b_date_proc
-                + pd.TimedeltaIndex(
-                    np.array([future_time - 3600.0], dtype=float), "s"
-                )[0]
+                + pd.TimedeltaIndex(np.array([future_time - 3600.0], dtype=float), "s")[
+                    0
+                ]
             )
         else:
             # Forecast-only mode — use BMI timestamp as-is
@@ -467,13 +469,19 @@ def loop_through_forcing_products(
                                     self.source_data_processor = NWMV3ConusProcessor(
                                         config_options, mpi_config, wrf_hydro_geo_meta
                                     )
-                                elif config_options.nwm_domain in [
-                                    "Hawaii",
-                                    "PR",
-                                ]:
-                                    self.source_data_processor = NWMV3OConusProcessor(
+                                elif config_options.nwm_domain == "Hawaii":
+                                    self.source_data_processor = NWMV3HawaiiProcessor(
                                         config_options, mpi_config, wrf_hydro_geo_meta
                                     )
+
+                                elif config_options.nwm_domain == "PR":
+                                    self.source_data_processor = (
+                                        NWMV3PuertoRicoProcessor(
+                                            config_options,
+                                            mpi_config,
+                                            wrf_hydro_geo_meta,
+                                        )
+                                    )
                                 elif config_options.nwm_domain == "Alaska":
                                     self.source_data_processor = NWMV3AlaskaProcessor(
                                         config_options, mpi_config, wrf_hydro_geo_meta

tests/test_data/expected_results/test_expected_bmi_model_after_update_n1_rank0__step_1.json

@@ -118,6 +118,10 @@
     "time_step_seconds": 3600
   },
   "geo_meta": {
+    "approx_centroid_global_xy": [
+      -72.05503025291092,
+      41.7594933573475
+    ],
     "centerCoords": null,
     "config_options": {
       "ExactExtract": null,
@@ -166,7 +170,7 @@
               "long_name": "longitude",
               "units": "degrees_east"
             },
-            "data": "hash_-7248319024934117860_len_14",
+            "data": "hash_-5399622063019475031_len_28",
             "dims": [
               "x"
             ]
@@ -176,7 +180,7 @@
               "long_name": "latitude",
               "units": "degrees_north"
             },
-            "data": "hash_5170267599511687273_len_29",
+            "data": "hash_-1457289664996315734_len_37",
             "dims": [
               "y"
             ]
@@ -313,8 +317,8 @@
           }
         },
         "dims": {
-          "x": 14,
-          "y": 29
+          "x": 28,
+          "y": 37
         }
       },
       "aws_time": null,
@@ -530,7 +534,7 @@
       "_coords": [
         [
           "hash_2904886930620536845_len_582",
-          "hash_2117205719057770460_len_582"
+          "hash_902249663378761532_len_582"
         ],
         [
           [
@@ -685,7 +689,7 @@
     "height_elem": null,
     "heights_global": null,
     "inds": null,
-    "lat_bounds": "hash_2117205719057770460_len_582",
+    "lat_bounds": "hash_902249663378761532_len_582",
     "latitude_grid": [
       41.67904491423198,
       41.73961258049449,
@@ -766,7 +770,7 @@
                 "long_name": "longitude",
                 "units": "degrees_east"
               },
-              "data": "hash_-7248319024934117860_len_14",
+              "data": "hash_-5399622063019475031_len_28",
               "dims": [
                 "x"
               ]
@@ -776,7 +780,7 @@
                 "long_name": "latitude",
                 "units": "degrees_north"
               },
-              "data": "hash_5170267599511687273_len_29",
+              "data": "hash_-1457289664996315734_len_37",
               "dims": [
                 "y"
               ]
@@ -913,8 +917,8 @@
             }
           },
           "dims": {
-            "x": 14,
-            "y": 29
+            "x": 28,
+            "y": 37
           }
         },
         "aws_time": null,