ABI products now support latitude/longitude bounding boxes

- Enable lat/lon bounding boxes for all ABI imagery, projection, navigation, and ancillary objects - Breaking API change: ABI navigation objects now take `lat_bounds` and `lon_bounds` arguments to constrain to individual pixels or pixel ranges. The deprecated arguments `lat_deg` and `lon_deg` are no longer supported and their functionality is maintained in `lat_bounds` and `lon_bounds`. - Fix bug where upscaled and subsetted ABI RGB images could be misaligned between channels - Add new function util.align_idx() to align ABI Fixed Grid indices and continuous slices between different spatial resolutions
heregoesradio · Feb 29, 2024 · caeea2e · caeea2e
1 parent f92eec2
commit caeea2e
Show file tree

Hide file tree

Showing 9 changed files with 701 additions and 350 deletions.
diff --git a/README.md b/README.md
@@ -21,7 +21,7 @@ This is the core processing library used by [Here GOES Radiotelescope (2020)](ht
 
 ## Setup
 
-Clone this repository and install the Conda environment. For Intel machines, use heregoes-env-intel.yml which uses MKL for acceleration. For other architectures, including ARM64 (e.g. Raspberry Pi 4), use heregoes-env-other.yml which installs with OpenBLAS.
+Clone this repository and install the Conda environment. For Intel machines, use heregoes-env-intel.yml which uses MKL for acceleration. For other architectures, including ARM64 (e.g. Raspberry Pi 5), use heregoes-env-other.yml which installs with OpenBLAS.
 ```
 conda env create -f heregoes-env-intel.yml
 conda activate heregoes-env
@@ -91,7 +91,7 @@ Which returns:
 
 ## Imagery Examples
 
-The heregoes library renders standard imagery for the GOES-R ABI and SUVI instruments following available literature wherever possible.
+The heregoes library renders standard imagery for the GOES-R ABI and SUVI instruments following available literature wherever possible. To reduce memory utilization, subsets of ABI images can be processed by specifying either an array index or continuous slice with the `index` argument, or with a latitude and longitude bounding box: `lat_bounds=[ul_lat, lr_lat]` and `lon_bounds=[ul_lon, lr_lon]`.
 
 ### Render a single-channel SUVI image
 
@@ -180,16 +180,16 @@ abi_navigation.sat_za #satellite zenith angle
 abi_navigation.sat_az #satellite azimuth in North-clockwise convention
 ```
 
-### Find a pixel in an ABI image from latitude and longitude
+### Find pixels in an ABI image from latitude and longitude
 
 With an ABI L1b radiance netCDF file `abi_nc`:
 
 ```python
 from heregoes import load, navigation
 
 abi_data = load(abi_nc)
-abi_navigation = navigation.ABINavigation(abi_data, lat_deg=44.72609499, lon_deg=-93.02279070)
-abi_navigation.index #nearest (y, x) index of the ABI image at provided latitude and longitude
+abi_navigation = navigation.ABINavigation(abi_data, lat_bounds=[44.731495, 44.602596], lon_bounds=[-93.01798, -92.85648])
+abi_navigation.index #nearest (y, x) index or slice of the ABI image at provided latitude and longitude bounds
 ```
 
 ---

diff --git a/heregoes/ancillary.py b/heregoes/ancillary.py
@@ -104,21 +104,26 @@ class IREMIS(AncillaryDataset):
 
     Arguments:
         - `abi_data`: The ABIObject formed on a GOES-R ABI L1b Radiance netCDF file as returned by `heregoes.load()`
-        - `index`: Optionally constrains the IREMIS dataset to an index or continuous slice on the ABI Fixed Grid matching the resolution of the provided `abi_data` object
+        - `index`: Optionally constrains the IREMIS dataset to an array index or continuous slice on the ABI Fixed Grid matching the resolution of the provided `abi_data` object
+        - `lat_bounds`, `lon_bounds`: Optionally constrains the IREMIS dataset to a latitude and longitude bounding box defined by the upper left and lower right points, e.g. `lat_bounds=[ul_lat, lr_lat]`, `lon_bounds=[ul_lon, lr_lon]`
         - `iremis_dir`: Location of IREMIS netCDF files. Defaults to the directory set by the HEREGOES_ENV_IREMIS_DIR environmental variable
     """
 
-    def __init__(self, abi_data, index=None, iremis_dir=IREMIS_DIR):
+    def __init__(
+        self,
+        abi_data,
+        index=None,
+        lat_bounds=None,
+        lon_bounds=None,
+        iremis_dir=IREMIS_DIR,
+    ):
         super().__init__()
 
         self.abi_data = abi_data
         self.index = index
         month = self.abi_data.time_coverage_start.month
         self.dataset_name = "iremis_month" + str(month).zfill(2)
 
-        if self.index is None:
-            self.index = np.s_[:, :]
-
         try:
             iremis_dir = Path(iremis_dir)
         except Exception as e:
@@ -128,6 +133,19 @@ def __init__(self, abi_data, index=None, iremis_dir=IREMIS_DIR):
                 exception=e,
             )
 
+        if lat_bounds is not None and lon_bounds is not None:
+            self.abi_nav = navigation.ABINavigation(
+                self.abi_data,
+                lat_bounds=np.atleast_1d(lat_bounds),
+                lon_bounds=np.atleast_1d(lon_bounds),
+            )
+            self.index = self.abi_nav.index
+            self.lat_deg = self.abi_nav.lat_deg
+            self.lon_deg = self.abi_nav.lon_deg
+
+        elif self.index is None:
+            self.index = np.s_[:, :]
+
         iremis_locations_nc = iremis_dir.joinpath("global_emis_inf10_location.nc")
         iremis_months = [
             "global_emis_inf10_monthFilled_MYD11C3.A2016001.041.nc",
@@ -189,12 +207,14 @@ def __init__(self, abi_data, index=None, iremis_dir=IREMIS_DIR):
         abi_projection = projection.ABIProjection(self.abi_data, index=self.index)
         self.data["c07_land_emissivity"] = abi_projection.resample2abi(
             self.data["c07_land_emissivity"],
-            latlon_bounds=[iremis_ul_lon, iremis_ul_lat, iremis_lr_lon, iremis_lr_lat],
+            lat_bounds=[iremis_ul_lat, iremis_lr_lat],
+            lon_bounds=[iremis_ul_lon, iremis_lr_lon],
             interpolation="bilinear",
         )
         self.data["c14_land_emissivity"] = abi_projection.resample2abi(
             self.data["c14_land_emissivity"],
-            latlon_bounds=[iremis_ul_lon, iremis_ul_lat, iremis_lr_lon, iremis_lr_lat],
+            lat_bounds=[iremis_ul_lat, iremis_lr_lat],
+            lon_bounds=[iremis_ul_lon, iremis_lr_lon],
             interpolation="bilinear",
         )
 
@@ -209,7 +229,8 @@ class WaterMask(AncillaryDataset):
 
     Arguments:
         - `abi_data`: The ABIObject formed on a GOES-R ABI L1b Radiance netCDF file as returned by `heregoes.load()`
-        - `index`: Optionally constrains the GSHHS dataset to an index or continuous slice on the ABI Fixed Grid matching the resolution of the provided `abi_data` object
+        - `index`: Optionally constrains the GSHHS dataset to an array index or continuous slice on the ABI Fixed Grid matching the resolution of the provided `abi_data` object
+        - `lat_bounds`, `lon_bounds`: Optionally constrains the GSHHS dataset to a latitude and longitude bounding box defined by the upper left and lower right points, e.g. `lat_bounds=[ul_lat, lr_lat]`, `lon_bounds=[ul_lon, lr_lon]`
         - `gshhs_scale`: 'auto', 'coarse', 'low', 'intermediate', 'high, or 'full' (https://scitools.org.uk/cartopy/docs/latest/reference/generated/cartopy.feature.GSHHSFeature.html)
         - `rivers`: Default `False`
     """
@@ -218,6 +239,8 @@ def __init__(
         self,
         abi_data,
         index=None,
+        lat_bounds=None,
+        lon_bounds=None,
         gshhs_scale="intermediate",
         rivers=False,
     ):
@@ -227,7 +250,18 @@ def __init__(
         self.index = index
         self.dataset_name = "gshhs_" + gshhs_scale
 
-        if self.index is None:
+        if lat_bounds is not None and lon_bounds is not None:
+            self.lat_deg = np.atleast_1d(lat_bounds)
+            self.lon_deg = np.atleast_1d(lon_bounds)
+
+            self.abi_nav = navigation.ABINavigation(
+                self.abi_data,
+                lat_bounds=self.lat_deg,
+                lon_bounds=self.lon_deg,
+            )
+            self.index = self.abi_nav.index
+
+        elif self.index is None:
             self.index = np.s_[:, :]
 
         y_rad = self.abi_data["y"][self.index[0]]

diff --git a/heregoes/exceptions.py b/heregoes/exceptions.py
@@ -34,43 +34,52 @@
 
 
 class HereGOESException(Exception):
-    def __init__(self, caller, msg, exception=None):
+    def __init__(self, msg="", caller=None, exception=None):
         if exception is not None:
             msg += f" Exception: {exception}"
-        logger.critical(msg, extra={"caller": caller})
+        if len(msg) > 0:
+            logger.critical(msg, extra={"caller": caller})
         super().__init__(msg)
 
+    def __str__(self):
+        return self.__class__.__qualname__
+
+    def __repr__(self):
+        return self.__str__()
+
 
 class HereGOESIOError(HereGOESException):
-    pass
+    def __init__(self, *args, **kwargs):
+        super().__init__(*args, **kwargs)
 
 
 class HereGOESIOReadException(HereGOESIOError):
     def __init__(self, caller, filepath, exception=None):
         super().__init__(
-            caller=caller,
             msg=f"Could not read from file at {filepath}.",
+            caller=caller,
             exception=exception,
         )
 
 
 class HereGOESIOWriteException(HereGOESIOError):
     def __init__(self, caller, filepath, exception=None):
         super().__init__(
-            caller=caller,
             msg=f"Could not write to file at {filepath}.",
+            caller=caller,
             exception=exception,
         )
 
 
 class HereGOESValueError(HereGOESException):
-    pass
+    def __init__(self, *args, **kwargs):
+        super().__init__(*args, **kwargs)
 
 
 class HereGOESUnsupportedProductException(HereGOESValueError):
     def __init__(self, caller, filepath, exception=None):
         super().__init__(
-            caller=caller,
             msg=f"Product type is not supported: {filepath}.",
+            caller=caller,
             exception=exception,
         )
diff --git a/heregoes/image.py b/heregoes/image.py
@@ -29,9 +29,9 @@
 from astropy.convolution import Gaussian2DKernel, interpolate_replace_nans
 from scipy import ndimage
 
-from heregoes import exceptions, load
+from heregoes import exceptions, load, navigation
 from heregoes.goesr import abi, suvi
-from heregoes.util import make_8bit, scale_idx
+from heregoes.util import align_idx, make_8bit, scale_idx
 
 logger = logging.getLogger("heregoes-logger")
 safe_time_format = "%Y-%m-%dT%H%M%SZ"
@@ -84,6 +84,8 @@ def __init__(
         self,
         abi_nc,
         index=None,
+        lat_bounds=None,
+        lon_bounds=None,
         gamma=1.0,
         black_space=False,
     ):
@@ -101,7 +103,8 @@ def __init__(
 
         Arguments:
             - `abi_nc`: String or Path object pointing to a GOES-R ABI L1b Radiance netCDF file
-            - `index`: Optionally constrains the ABI image to an index or continuous slice on the ABI Fixed Grid matching the resolution in the provided `abi_nc` file
+            - `index`: Optionally constrains the ABI image to an array index or continuous slice on the ABI Fixed Grid matching the resolution in the provided `abi_nc` file
+            - `lat_bounds`, `lon_bounds`: Optionally constrains the ABI image to a latitude and longitude bounding box defined by the upper left and lower right points, e.g. `lat_bounds=[ul_lat, lr_lat]`, `lon_bounds=[ul_lon, lr_lon]`
             - `gamma`: Optional gamma correction for reflective ABI brightness value. Defaults to no correction
             - `black_space`: Optionally overwrites the masked pixels in the final ABI image (nominally the "space" background) to be black. Defaults to no overwriting, or white pixels for reflective imagery and black pixels for emissive imagery. Default `True`
         """
@@ -115,7 +118,17 @@ def __init__(
 
         self.abi_data = load(abi_nc)
 
-        if self.index is None:
+        if lat_bounds is not None and lon_bounds is not None:
+            self.abi_nav = navigation.ABINavigation(
+                self.abi_data,
+                lat_bounds=np.atleast_1d(lat_bounds),
+                lon_bounds=np.atleast_1d(lon_bounds),
+            )
+            self.index = self.abi_nav.index
+            self.lat_deg = self.abi_nav.lat_deg
+            self.lon_deg = self.abi_nav.lon_deg
+
+        elif self.index is None:
             self.index = np.s_[:, :]
 
         self.rad = self.abi_data["Rad"][self.index]
@@ -198,6 +211,8 @@ def __init__(
         green_nc,
         blue_nc,
         index=None,
+        lat_bounds=None,
+        lon_bounds=None,
         r_coeff=0.45,
         g_coeff=0.1,
         b_coeff=0.45,
@@ -211,7 +226,8 @@ def __init__(
 
         Arguments:
             - `red_nc`, `green_nc`, `blue_nc`: Strings or Path objects pointing to GOES-R ABI L1b Radiance netCDF files for the red (0.64 μm), green (0.86 μm), and blue (0.47 μm) components
-            - `index`: Optionally constrains the ABI imagery to an index or continuous slice on the ABI Fixed Grid. If `upscale` is `False` (default), this is an index on the 1 km Fixed Grid. Otherwise, if `upscale` is `True`, a 500 m Fixed Grid is used
+            - `index`: Optionally constrains the ABI imagery to an array index or continuous slice on the ABI Fixed Grid. If `upscale` is `False` (default), this is an index on the 1 km Fixed Grid. Otherwise, if `upscale` is `True`, a 500 m Fixed Grid is used
+            - `lat_bounds`, `lon_bounds`: Optionally constrains the ABI imagery to a latitude and longitude bounding box defined by the upper left and lower right points, e.g. `lat_bounds=[ul_lat, lr_lat]`, `lon_bounds=[ul_lon, lr_lon]`
             - `r_coeff`, `g_coeff`, `b_coeff`: Coefficients for the fractional combination "green" band method described in Bah et. al (2018)
             - `upscale`: Whether to scale up green and blue images (1 km) to match the red image (500 m) (`True`) or vice versa (`False`, Default)
             - `upscale_algo`: The OpenCV interpolation algorithm used for upscaling green and blue images. See https://docs.opencv.org/4.6.0/da/d54/group__imgproc__transform.html#ga5bb5a1fea74ea38e1a5445ca803ff121. Default `cv2.INTER_CUBIC`
@@ -222,17 +238,42 @@ def __init__(
         super().__init__()
         self.index = index
 
-        if self.index is None:
-            self.index = np.s_[:, :]
-
         if upscale:
+            nav_nc = red_nc
             scaler_500m = 1.0
             scaler_1km = 2.0
 
         else:
+            nav_nc = green_nc
             scaler_500m = 0.5
             scaler_1km = 1.0
 
+        if lat_bounds is not None and lon_bounds is not None:
+            self.abi_nav = navigation.ABINavigation(
+                load(nav_nc),
+                lat_bounds=np.atleast_1d(lat_bounds),
+                lon_bounds=np.atleast_1d(lon_bounds),
+            )
+            self.index = self.abi_nav.index
+            self.lat_deg = self.abi_nav.lat_deg
+            self.lon_deg = self.abi_nav.lon_deg
+
+        elif self.index is None:
+            self.index = np.s_[:, :]
+
+        if upscale:
+            aligned_idx = align_idx(self.index, 2)
+            if self.index != aligned_idx:
+                # only warn for index alignment if index was provided
+                if lat_bounds is None or lon_bounds is None:
+                    logger.warning(
+                        "Adjusting provided index %s to %s to align to the 1 km ABI Fixed Grid.",
+                        str(self.index),
+                        str(aligned_idx),
+                        extra={"caller": f"{__name__}.{self.__class__.__name__}"},
+                    )
+                self.index = aligned_idx
+
         # make images
         red_image = ABIImage(
             red_nc,