diff --git a/README.md b/README.md
index c18b1b0..62b60c2 100644
--- 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,7 +180,7 @@ 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`:
@@ -188,8 +188,8 @@ With an ABI L1b radiance netCDF file `abi_nc`:
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
index b7a9600..7803ed0 100755
--- a/heregoes/ancillary.py
+++ b/heregoes/ancillary.py
@@ -104,11 +104,19 @@ 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
@@ -116,9 +124,6 @@ def __init__(self, abi_data, index=None, iremis_dir=IREMIS_DIR):
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
index 3b44bc9..b1665e7 100644
--- a/heregoes/exceptions.py
+++ b/heregoes/exceptions.py
@@ -34,22 +34,30 @@
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,
)
@@ -57,20 +65,21 @@ def __init__(self, caller, filepath, exception=None):
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
index 41f95a8..a7425dd 100755
--- 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,
diff --git a/heregoes/navigation.py b/heregoes/navigation.py
index 8d743c8..7867d12 100755
--- a/heregoes/navigation.py
+++ b/heregoes/navigation.py
@@ -18,19 +18,24 @@
# You should have received a copy of the GNU General Public License
# along with this program. If not, see .
+import logging
+
import astropy.units as u
import numpy as np
from astropy import coordinates
from astropy.time import Time
+from heregoes import exceptions
from heregoes.util import nearest_2d, njit, orbital, rad2deg
+logger = logging.getLogger("heregoes-logger")
+
class ABINavigation:
"""
This is a class for GOES-R ABI navigation routines using an ABIObject to access netCDF data.
- The routines may be constrained to a single index or continuous slice of the ABI Fixed Grid with `index`, or to a single area by providing `lat_deg` and `lon_deg` (degrees).
- All calculations return 32-bit floating point NumPy arrays which should be accurate enough for most applications at this scale.
+ The routines may be constrained to a single index or continuous slice of the ABI Fixed Grid with `index`, or to a single area by providing `lat_bounds` and `lon_bounds` (degrees).
+ All calculations return 32-bit floating point NumPy arrays which should be accurate enough for most applications at ABI scale.
The following quantities are always generated upon instantiation:
- Instrument scanning angle (`y_rad`, `x_rad`)
@@ -44,6 +49,8 @@ class ABINavigation:
Arguments:
- `abi_data`: The ABIObject formed on a GOES-R ABI L1b Radiance netCDF file as returned by `heregoes.load()`
+ - `index`: Optionally constrains the navigation routines to an array index or continuous slice on the ABI Fixed Grid
+ - `lat_bounds`, `lon_bounds`: Optionally constrains the navigation routines 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]`
- `hae_m`: The Height Above Ellipsoid (HAE) in meters of the ABI array to correct for terrain height. Default 0.0 (no correction)
- `time`: The time for which the Sun position is valid. The product midpoint time is used if not provided
- `precise_sun`: Whether to calculate solar position using Equation of Time with Pyorbital (`False`, default) or real ephemeris with Astropy (`True`)
@@ -54,8 +61,8 @@ def __init__(
self,
abi_data,
index=None,
- lat_deg=None,
- lon_deg=None,
+ lat_bounds=None,
+ lon_bounds=None,
hae_m=0.0,
time=None,
precise_sun=False,
@@ -77,12 +84,15 @@ def __init__(
if self.index is None:
self.index = np.s_[:, :]
+ if self.time is None:
+ self.time = self.abi_data.midpoint_time
+
self.x_rad, self.y_rad = np.meshgrid(
self.abi_data["x"][self.index[1]],
self.abi_data["y"][self.index[0]],
)
- if lat_deg is None or lon_deg is None:
+ if lat_bounds is None or lon_bounds is None:
self.lat_deg, self.lon_deg = self.navigate(
self.y_rad,
self.x_rad,
@@ -96,9 +106,6 @@ def __init__(
].perspective_point_height,
)
- if self.hae_m.shape != self.lat_deg.shape:
- self.hae_m = np.full(self.lat_deg.shape, self.hae_m, dtype=np.float32)
-
# correct for terrain parallax if HAE is provided
if (self.hae_m != 0.0).any():
self.y_rad, self.x_rad = self.reverse_navigate(
@@ -112,7 +119,7 @@ def __init__(
sat_height=self.abi_data[
"goes_imager_projection"
].perspective_point_height,
- feature_height=self.hae_m,
+ feature_height=np.broadcast_to(self.hae_m, self.lat_deg.shape),
)
self.lat_deg, self.lon_deg = self.navigate(
self.y_rad,
@@ -128,12 +135,32 @@ def __init__(
)
else:
- self.lat_deg = np.atleast_1d(lat_deg)
- self.lon_deg = np.atleast_1d(lon_deg)
+ lat_bounds = np.atleast_1d(lat_bounds).astype(np.float32)
+ lon_bounds = np.atleast_1d(lon_bounds).astype(np.float32)
+
+ if lat_bounds.shape != lon_bounds.shape:
+ raise exceptions.HereGOESValueError(
+ msg="`lat_bounds` and `lon_bounds` must be the same shape.",
+ caller=f"{__name__}.{self.__class__.__name__}",
+ )
+
+ if np.isnan(lat_bounds).any() | np.isnan(lon_bounds).any():
+ logger.warning(
+ "Provided lat/lon bounds contain NaN; check bounds for off-Earth pixels.",
+ extra={"caller": f"{__name__}.{self.__class__.__name__}"},
+ )
+
+ # this is tricky to deal with because of Numba broadcasting issues, and the desired shape of hae_m being unknown at this step when using lat/lon bounds
+ # related: https://github.com/numba/numba/issues/4632
+ if (self.hae_m != 0.0).any():
+ logger.warning(
+ "Terrain height correction is not currently supported when `lat_bounds` and `lon_bounds` are specified; ignoring argument `hae_m`.",
+ extra={"caller": f"{__name__}.{self.__class__.__name__}"},
+ )
derived_y_rad, derived_x_rad = self.reverse_navigate(
- self.lat_deg,
- self.lon_deg,
+ lat_bounds,
+ lon_bounds,
lon_origin=self.abi_data[
"goes_imager_projection"
].longitude_of_projection_origin,
@@ -142,27 +169,33 @@ def __init__(
sat_height=self.abi_data[
"goes_imager_projection"
].perspective_point_height,
- feature_height=self.hae_m,
+ feature_height=np.broadcast_to(
+ np.atleast_1d(0.0).astype(np.float32), lat_bounds.shape
+ ),
)
+
self.index = nearest_2d(
self.y_rad, self.x_rad, derived_y_rad, derived_x_rad
)
- # form a tuple of slices to encompass a continuous range
- if len(self.index[0]) > 1 or len(self.index[1]) > 1:
- self.index = tuple(
- slice(idx.T.min(), idx.T.max() + 1) for idx in self.index
- )
-
- # or form a single tuple index
- else:
- self.index = (self.index[0].item(), self.index[1].item())
+ self.y_rad = np.atleast_1d(self.y_rad[self.index]).astype(np.float32)
+ self.x_rad = np.atleast_1d(self.x_rad[self.index]).astype(np.float32)
- self.y_rad = derived_y_rad
- self.x_rad = derived_x_rad
+ self.lat_deg, self.lon_deg = self.navigate(
+ self.y_rad,
+ self.x_rad,
+ lon_origin=self.abi_data[
+ "goes_imager_projection"
+ ].longitude_of_projection_origin,
+ r_eq=self.abi_data["goes_imager_projection"].semi_major_axis,
+ r_pol=self.abi_data["goes_imager_projection"].semi_minor_axis,
+ sat_height=self.abi_data[
+ "goes_imager_projection"
+ ].perspective_point_height,
+ )
- if self.time is None:
- self.time = self.abi_data.midpoint_time
+ if self.hae_m.shape != self.lat_deg.shape:
+ self.hae_m = np.broadcast_to(self.hae_m, self.lat_deg.shape)
@property
def sat_za(self):
diff --git a/heregoes/projection.py b/heregoes/projection.py
index d27a38b..1429814 100755
--- a/heregoes/projection.py
+++ b/heregoes/projection.py
@@ -27,7 +27,7 @@
gdal.UseExceptions()
-from heregoes import NUM_CPUS
+from heregoes import NUM_CPUS, navigation
GDAL_PARALLEL = False
if os.getenv("HEREGOES_ENV_PARALLEL", "False").lower() == "true":
@@ -40,7 +40,9 @@ class ABIProjection:
Arguments:
- `abi_data`: The ABIObject formed on a GOES-R ABI L1b Radiance netCDF file as returned by `heregoes.load()`
- - `index`: Optionally constrains the projection to an index or continuous slice on the ABI Fixed Grid matching the resolution of the provided `abi_data` object
+ - `index`: Optionally constrains the projection 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 projection 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]`
+ - When projecting to or from a subset of an ABI image, the `index` of the image is preferred to using `lat_bounds` and `lon_bounds` to ensure the image bounds match exactly on the ABI Fixed Grid
Class methods:
- `resample2abi(latlon_array)` resamples an array with WGS84 lat/lon projection to the ABI Fixed Grid domain. Returns the resampled array if convert_np is `True` (default), otherwise returns the GDAL dataset
@@ -48,11 +50,21 @@ class ABIProjection:
- `resample2cog(abi_array, cog_filepath)` resamples an ABI array from the ABI Fixed Grid domain to WGS84 lat/lon projection and saves to a Cloud Optimized GeoTIFF (COG) at the filepath `cog_filepath`
"""
- def __init__(self, abi_data, index=None):
+ def __init__(self, abi_data, index=None, lat_bounds=None, lon_bounds=None):
self.abi_data = abi_data
self.index = index
- 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_[:, :]
h = self.abi_data["goes_imager_projection"].perspective_point_height
@@ -183,15 +195,19 @@ def _resample(self, ds, translate_options, warp_options):
def resample2abi(
self,
latlon_array,
- latlon_bounds=[-180.0, 90.0, 180.0, -90.0],
+ lat_bounds=[90.0, -90.0],
+ lon_bounds=[-180.0, 180.0],
interpolation="nearest",
convert_np=True,
):
+ ul_y, lr_y = lat_bounds
+ ul_x, lr_x = lon_bounds
+
ds = self._make_dataset(latlon_array)
translate_options = gdal.TranslateOptions(
outputSRS=self.latlon_srs,
- outputBounds=latlon_bounds,
+ outputBounds=[ul_x, ul_y, lr_x, lr_y],
format=self._intermediate_format,
resampleAlg=interpolation.lower(),
creationOptions=self._intermediate_gdal_options,
diff --git a/heregoes/util/__init__.py b/heregoes/util/__init__.py
index 96ca316..4d1b5e2 100644
--- a/heregoes/util/__init__.py
+++ b/heregoes/util/__init__.py
@@ -76,10 +76,49 @@ def x4(arr):
return nearest_scale(arr, k=4)
+def align_idx(idx, modulus):
+ """
+ Aligns an array index or continuous slice `idx` to the closest index modulo `modulus`.
+ Useful for aligning indices between different spatial resolutions, e.g. aligning a 500 m index to the 1 km ABI Fixed Grid with a `modulus` of 2, or to the 2 km Fixed Grid with a `modulus` of 4.
+ """
+
+ def safe_align(value, modulus):
+ if value is None or value is Ellipsis:
+ return value
+
+ else:
+ i = value - (value % modulus)
+ j = (value + modulus) - (value % modulus)
+
+ if j - value < value - i:
+ return int(j)
+
+ else:
+ return int(i)
+
+ def slice_align(slc, modulus):
+ start = safe_align(slc.start, modulus)
+ stop = safe_align(slc.stop, modulus)
+
+ return np.s_[start:stop:None]
+
+ aligned_idx = []
+ for i in tuple(idx):
+ if isinstance(i, slice):
+ slc = slice_align(i, modulus)
+
+ else:
+ slc = safe_align(i, modulus)
+
+ aligned_idx.append(slc)
+
+ return tuple(aligned_idx)
+
+
def scale_idx(idx, scale_factor):
"""
- Adjusts an array index or slice by `scale_factor`.
- Useful for referencing the same point between different spatial resolutions, e.g. converting a 500 m index for use with a 2 km product using a scale factor of 0.25.
+ Adjusts an array index or continous slice `idx` by `scale_factor`.
+ Useful for referencing the same point between different spatial resolutions, e.g. converting a 500 m index for use with a 2 km product using a `scale_factor` of 0.25.
"""
def safe_floor(value):
@@ -92,19 +131,20 @@ def safe_floor(value):
def slice_floor(slc):
start = safe_floor(slc.start)
stop = safe_floor(slc.stop)
- step = safe_floor(slc.step)
- return np.s_[start:stop:step]
+ return np.s_[start:stop:None]
idx = np.s_[idx]
scaled_idx = []
for i in tuple(idx):
if isinstance(i, slice):
- scaled_idx.append(slice_floor(i))
+ slc = slice_floor(i)
else:
- scaled_idx.append(safe_floor(i))
+ slc = safe_floor(i)
+
+ scaled_idx.append(slc)
return tuple(scaled_idx)
@@ -240,7 +280,6 @@ def unravel_index(index, shape):
return to_fixed_tuple(result, len(shape))
-@njit.heregoes_njit_noparallel
def nearest_2d(y_arr, x_arr, target_y, target_x):
"""
Finds the nearest index of a value simultaneously for two arrays `y_arr` and `x_arr` of the same shape.
@@ -248,35 +287,51 @@ def nearest_2d(y_arr, x_arr, target_y, target_x):
Adapted from: https://github.com/blaylockbk/pyBKB_v3/blob/master/demo/Nearest_lat-lon_Grid.ipynb (MIT)
"""
- if y_arr.shape != x_arr.shape:
- return
+ @njit.heregoes_njit_noparallel
+ def _nearest_2d(y_arr, x_arr, target_y, target_x):
+ if y_arr.shape != x_arr.shape:
+ return
+
+ if np.isnan(y_arr).any() | np.isnan(x_arr).any():
+ y_arr = y_arr.copy()
+ x_arr = x_arr.copy()
+
+ y_arr.ravel()[np.nonzero(np.isnan(y_arr.ravel()))] = np.inf
+ x_arr.ravel()[np.nonzero(np.isnan(x_arr.ravel()))] = np.inf
+
+ nearest_ys = np.zeros(target_y.ravel().shape, dtype=np.int64)
+ nearest_xs = np.zeros(target_x.ravel().shape, dtype=np.int64)
+
+ flatidx = 0
+ for idx, val in np.ndenumerate(target_y):
+ nearest_y, nearest_x = unravel_index(
+ int(
+ np.argmin(
+ np.maximum(
+ np.abs(y_arr - target_y[idx]), np.abs(x_arr - target_x[idx])
+ )
+ )
+ ),
+ y_arr.shape,
+ )
- if np.isnan(y_arr).any() | np.isnan(x_arr).any():
- y_arr = y_arr.copy()
- x_arr = x_arr.copy()
+ nearest_ys[flatidx] = nearest_y
+ nearest_xs[flatidx] = nearest_x
- y_arr.ravel()[np.nonzero(np.isnan(y_arr.ravel()))] = np.inf
- x_arr.ravel()[np.nonzero(np.isnan(x_arr.ravel()))] = np.inf
+ flatidx += 1
- nearest_ys = np.zeros(target_y.ravel().shape, dtype=np.int64)
- nearest_xs = np.zeros(target_x.ravel().shape, dtype=np.int64)
+ return (nearest_ys, nearest_xs)
- flatidx = 0
- for idx, val in np.ndenumerate(target_y):
- nearest_y, nearest_x = unravel_index(
- int(
- np.argmin(
- np.maximum(
- np.abs(y_arr - target_y[idx]), np.abs(x_arr - target_x[idx])
- )
- )
- ),
- y_arr.shape,
- )
+ nearest_indices = _nearest_2d(y_arr, x_arr, target_y, target_x)
- nearest_ys[flatidx] = nearest_y
- nearest_xs[flatidx] = nearest_x
+ # form a tuple of slices to encompass a continuous range
+ if len(nearest_indices[0]) > 1 or len(nearest_indices[1]) > 1:
+ nearest_indices = tuple(
+ slice(idx.T.min(), idx.T.max() + 1) for idx in nearest_indices
+ )
- flatidx += 1
+ # or form a single tuple index
+ else:
+ nearest_indices = (nearest_indices[0].item(), nearest_indices[1].item())
- return (nearest_ys, nearest_xs)
+ return nearest_indices
diff --git a/test/heregoes_test.py b/test/heregoes_test.py
index 650866d..3c4b085 100755
--- a/test/heregoes_test.py
+++ b/test/heregoes_test.py
@@ -18,141 +18,28 @@
# You should have received a copy of the GNU General Public License
# along with this program. If not, see .
-from pathlib import Path
+import gc
import cv2
import numpy as np
from heregoes import ancillary, exceptions, image, load, navigation, projection
-from heregoes.util import minmax
+from heregoes.util import minmax, scale_idx
-SCRIPT_PATH = Path(__file__).parent.resolve()
+from .test_resources import *
-input_dir = SCRIPT_PATH.joinpath("input")
-input_dir.mkdir(parents=True, exist_ok=True)
output_dir = SCRIPT_PATH.joinpath("output")
output_dir.mkdir(parents=True, exist_ok=True)
for output_file in output_dir.glob("*"):
output_file.unlink()
-# abi
-abi_mc01_nc = input_dir.joinpath(
- "abi/OR_ABI-L1b-RadM1-M6C01_G16_s20211691942252_e20211691942310_c20211691942342.nc"
-)
-abi_mc02_nc = input_dir.joinpath(
- "abi/OR_ABI-L1b-RadM1-M6C02_G16_s20211691942252_e20211691942310_c20211691942334.nc"
-)
-abi_mc03_nc = input_dir.joinpath(
- "abi/OR_ABI-L1b-RadM1-M6C03_G16_s20211691942252_e20211691942310_c20211691942351.nc"
-)
-abi_mc04_nc = input_dir.joinpath(
- "abi/OR_ABI-L1b-RadM1-M6C04_G16_s20211691942252_e20211691942310_c20211691942340.nc"
-)
-abi_mc05_nc = input_dir.joinpath(
- "abi/OR_ABI-L1b-RadM1-M6C05_G16_s20211691942252_e20211691942310_c20211691942347.nc"
-)
-abi_mc06_nc = input_dir.joinpath(
- "abi/OR_ABI-L1b-RadM1-M6C06_G16_s20211691942252_e20211691942315_c20211691942345.nc"
-)
-abi_mc07_nc = input_dir.joinpath(
- "abi/OR_ABI-L1b-RadM1-M6C07_G16_s20211691942252_e20211691942321_c20211691942355.nc"
-)
-abi_mc08_nc = input_dir.joinpath(
- "abi/OR_ABI-L1b-RadM1-M6C08_G16_s20211691942252_e20211691942310_c20211691942357.nc"
-)
-abi_mc09_nc = input_dir.joinpath(
- "abi/OR_ABI-L1b-RadM1-M6C09_G16_s20211691942252_e20211691942315_c20211691942368.nc"
-)
-abi_mc10_nc = input_dir.joinpath(
- "abi/OR_ABI-L1b-RadM1-M6C10_G16_s20211691942252_e20211691942322_c20211691942353.nc"
-)
-abi_mc11_nc = input_dir.joinpath(
- "abi/OR_ABI-L1b-RadM1-M6C11_G16_s20211691942252_e20211691942310_c20211691942348.nc"
-)
-abi_mc12_nc = input_dir.joinpath(
- "abi/OR_ABI-L1b-RadM1-M6C12_G16_s20211691942252_e20211691942316_c20211691942356.nc"
-)
-abi_mc13_nc = input_dir.joinpath(
- "abi/OR_ABI-L1b-RadM1-M6C13_G16_s20211691942252_e20211691942322_c20211691942361.nc"
-)
-abi_mc14_nc = input_dir.joinpath(
- "abi/OR_ABI-L1b-RadM1-M6C14_G16_s20211691942252_e20211691942310_c20211691942364.nc"
-)
-abi_mc15_nc = input_dir.joinpath(
- "abi/OR_ABI-L1b-RadM1-M6C15_G16_s20211691942252_e20211691942316_c20211691942358.nc"
-)
-abi_mc16_nc = input_dir.joinpath(
- "abi/OR_ABI-L1b-RadM1-M6C16_G16_s20211691942252_e20211691942322_c20211691942366.nc"
-)
-
-# add a few g16 conuses to test off-earth pixels
-abi_cc01_nc = input_dir.joinpath(
- "abi/OR_ABI-L1b-RadC-M6C01_G16_s20211691941174_e20211691943547_c20211691943589.nc"
-)
-abi_cc02_nc = input_dir.joinpath(
- "abi/OR_ABI-L1b-RadC-M6C02_G16_s20211691941174_e20211691943547_c20211691943571.nc"
-)
-abi_cc03_nc = input_dir.joinpath(
- "abi/OR_ABI-L1b-RadC-M6C03_G16_s20211691941174_e20211691943547_c20211691943587.nc"
-)
-abi_cc07_nc = input_dir.joinpath(
- "abi/OR_ABI-L1b-RadC-M6C07_G16_s20211691941174_e20211691943558_c20211691944002.nc"
-)
-
-abi_ncs = [
- abi_cc01_nc,
- abi_cc02_nc,
- abi_cc03_nc,
- abi_cc07_nc,
- abi_mc01_nc,
- abi_mc02_nc,
- abi_mc03_nc,
- abi_mc04_nc,
- abi_mc05_nc,
- abi_mc06_nc,
- abi_mc07_nc,
- abi_mc08_nc,
- abi_mc09_nc,
- abi_mc10_nc,
- abi_mc11_nc,
- abi_mc12_nc,
- abi_mc13_nc,
- abi_mc14_nc,
- abi_mc15_nc,
- abi_mc16_nc,
-]
-
-# suvi
-suvi_094_nc = input_dir.joinpath(
- "suvi/OR_SUVI-L1b-Fe093_G16_s20203160623501_e20203160623511_c20203160624091.nc"
-)
-suvi_131_nc = input_dir.joinpath(
- "suvi/OR_SUVI-L1b-Fe131_G16_s20203160623001_e20203160623011_c20203160623196.nc"
-)
-suvi_171_nc = input_dir.joinpath(
- "suvi/OR_SUVI-L1b-Fe171_G16_s20203160624201_e20203160624211_c20203160624396.nc"
-)
-suvi_195_nc = input_dir.joinpath(
- "suvi/OR_SUVI-L1b-Fe195_G16_s20203160623301_e20203160623311_c20203160623491.nc"
-)
-suvi_284_nc = input_dir.joinpath(
- "suvi/OR_SUVI-L1b-Fe284_G16_s20203160624501_e20203160624511_c20203160625090.nc"
-)
-suvi_304_nc = input_dir.joinpath(
- "suvi/OR_SUVI-L1b-He303_G16_s20203160622501_e20203160622511_c20203160623090.nc"
-)
-suvi_ncs = [
- suvi_094_nc,
- suvi_131_nc,
- suvi_171_nc,
- suvi_195_nc,
- suvi_284_nc,
- suvi_304_nc,
-]
+epsilon_degrees = 1e-4
+epsilon_meters = 1
def test_abi_image():
+ # test single-channel images
for abi_nc in abi_ncs:
abi_image = image.ABIImage(abi_nc, gamma=0.75)
@@ -162,26 +49,100 @@ def test_abi_image():
abi_image.save(file_path=output_dir, file_ext=".jpg")
- abi_rgb = image.ABINaturalRGB(
- abi_mc02_nc,
- abi_mc03_nc,
- abi_mc01_nc,
- upscale=True,
- gamma=0.75,
- black_space=True,
- )
-
- assert abi_rgb.bv.dtype == np.uint8
-
- abi_rgb.save(file_path=output_dir, file_ext=".jpeg")
-
- abi_rgb = image.ABINaturalRGB(
- abi_cc02_nc, abi_cc03_nc, abi_cc01_nc, gamma=0.75, black_space=True
- )
-
- assert abi_rgb.bv.dtype == np.uint8
-
- abi_rgb.save(file_path=output_dir, file_ext=".jpeg")
+ # test index alignment for subsetted RGB images
+ lat_bounds_500m = (46.0225830078125, 43.89013671875)
+ lon_bounds_500m = [-94.68467712402344, -91.75820922851562]
+ lat_bounds_1km = [46.02677536010742, 43.90188217163086]
+ lon_bounds_1km = (-94.6901626586914, -91.77256774902344)
+
+ for scene in ["meso", "conus"]:
+ if scene == "meso":
+ slc_500m = np.s_[213:474, 11:307]
+ r_nc = abi_mc02_nc
+ g_nc = abi_mc03_nc
+ b_nc = abi_mc01_nc
+
+ elif scene == "conus":
+ slc_500m = np.s_[613:875, 4451:4747]
+ r_nc = abi_cc02_nc
+ g_nc = abi_cc03_nc
+ b_nc = abi_cc01_nc
+
+ slc_1km = scale_idx(slc_500m, 0.5)
+
+ for upscale in [True, False]:
+ if upscale:
+ slc = slc_500m
+ lat_bounds = lat_bounds_500m
+ lon_bounds = lon_bounds_500m
+ else:
+ slc = slc_1km
+ lat_bounds = lat_bounds_1km
+ lon_bounds = lon_bounds_1km
+
+ # full RGB
+ abi_rgb_full = image.ABINaturalRGB(
+ r_nc,
+ g_nc,
+ b_nc,
+ upscale=upscale,
+ gamma=0.75,
+ black_space=True,
+ )
+
+ assert abi_rgb_full.bv.dtype == np.uint8
+
+ abi_rgb_full.save(file_path=output_dir, file_ext=".jpeg")
+
+ # indexed RGB
+ abi_rgb_indexed_bounds = image.ABINaturalRGB(
+ r_nc,
+ g_nc,
+ b_nc,
+ index=slc,
+ upscale=upscale,
+ gamma=0.75,
+ black_space=True,
+ )
+
+ # latlon RGB
+ abi_rgb_latlon_bounds = image.ABINaturalRGB(
+ r_nc,
+ g_nc,
+ b_nc,
+ lat_bounds=lat_bounds,
+ lon_bounds=lon_bounds,
+ upscale=upscale,
+ gamma=0.75,
+ black_space=True,
+ )
+
+ # get the original index of the brightest pixel within the slice
+ brightest_idx_500m = np.unravel_index(
+ np.nanargmax(np.sum(abi_rgb_full.bv[slc], axis=2)),
+ abi_rgb_full.bv[slc].shape[0:2],
+ )
+
+ # if the RGB image is upscaled, then the 500m slice of the below subsetted images will have been aligned +1,+1 pixels to the 1 km FGF
+ brightest_idx_500m_aligned = tuple(
+ [i + upscale for i in brightest_idx_500m]
+ )
+ assert (
+ brightest_idx_500m_aligned
+ == np.unravel_index(
+ np.nanargmax(np.sum(abi_rgb_indexed_bounds.bv, axis=2)),
+ abi_rgb_indexed_bounds.bv.shape[0:2],
+ )
+ == np.unravel_index(
+ np.nanargmax(np.sum(abi_rgb_latlon_bounds.bv, axis=2)),
+ abi_rgb_latlon_bounds.bv.shape[0:2],
+ )
+ )
+
+ del abi_rgb_full
+ del abi_rgb_indexed_bounds
+ del abi_rgb_latlon_bounds
+ _ = gc.collect()
def test_suvi_image():
@@ -199,10 +160,12 @@ def test_projection():
abi_mc02_nc,
abi_mc03_nc,
abi_mc01_nc,
+ lat_bounds=[45.4524291240206, 44.567740562044825],
+ lon_bounds=[-93.86441304735817, -92.69697639796475],
upscale=True,
gamma=1.0,
)
- abi_projection = projection.ABIProjection(abi_rgb.abi_data)
+ abi_projection = projection.ABIProjection(abi_rgb.abi_data, index=abi_rgb.index)
abi_rgb = abi_projection.resample2cog(
abi_rgb.bv, output_dir.joinpath(abi_rgb.default_filename + ".tiff")
)
@@ -211,173 +174,252 @@ def test_projection():
abi_rgb = image.ABINaturalRGB(
abi_cc02_nc, abi_cc03_nc, abi_cc01_nc, index=slc, upscale=True, gamma=0.75
)
- abi_projection = projection.ABIProjection(abi_rgb.abi_data, index=slc)
+ abi_projection = projection.ABIProjection(abi_rgb.abi_data, index=abi_rgb.index)
abi_rgb = abi_projection.resample2cog(
abi_rgb.bv, output_dir.joinpath(abi_rgb.default_filename + ".tiff")
)
-def test_navigation():
+def test_point_navigation():
idx = (92, 42)
- abi_data = load(abi_mc07_nc)
+ meso_data = load(abi_mc07_nc)
- abi_nav = navigation.ABINavigation(abi_data, precise_sun=False)
+ meso_nav = navigation.ABINavigation(meso_data, precise_sun=False)
- assert abi_nav.lat_deg.dtype == np.float32
- assert abi_nav.lon_deg.dtype == np.float32
- assert abi_nav.lat_deg[idx] == 44.73149490356445
- assert abi_nav.lon_deg[idx] == -93.01798248291016
+ assert meso_nav.lat_deg.dtype == np.float32
+ assert meso_nav.lon_deg.dtype == np.float32
+ assert meso_nav.lat_deg[idx] == 44.73149490356445
+ assert meso_nav.lon_deg[idx] == -93.01798248291016
- assert abi_nav.area_m.dtype == np.float32
- assert abi_nav.area_m[idx] == 4398247.5
+ assert meso_nav.area_m.dtype == np.float32
+ assert meso_nav.area_m[idx] == 4398247.5
- assert abi_nav.sat_za.dtype == np.float32
- assert abi_nav.sat_az.dtype == np.float32
- assert abi_nav.sat_za[idx] == 0.9509874582290649
- assert abi_nav.sat_az[idx] == 2.7129034996032715
+ assert meso_nav.sat_za.dtype == np.float32
+ assert meso_nav.sat_az.dtype == np.float32
+ assert meso_nav.sat_za[idx] == 0.9509874582290649
+ assert meso_nav.sat_az[idx] == 2.7129034996032715
- assert abi_nav.sun_za.dtype == np.float32
- assert abi_nav.sun_az.dtype == np.float32
- assert abi_nav.sun_za[idx] == 0.4886804521083832
- assert abi_nav.sun_az[idx] == 3.9760777950286865
+ assert meso_nav.sun_za.dtype == np.float32
+ assert meso_nav.sun_az.dtype == np.float32
+ assert meso_nav.sun_za[idx] == 0.4886804521083832
+ assert meso_nav.sun_az[idx] == 3.9760777950286865
# test on astropy sun
- abi_nav = navigation.ABINavigation(abi_data, index=idx, precise_sun=True)
+ meso_nav = navigation.ABINavigation(meso_data, index=idx, precise_sun=True)
- assert abi_nav.sun_za.dtype == np.float32
- assert abi_nav.sun_az.dtype == np.float32
- assert abi_nav.sun_za[0] == 0.4887488782405853
- assert abi_nav.sun_az[0] == 3.976273775100708
+ assert meso_nav.sun_za.dtype == np.float32
+ assert meso_nav.sun_az.dtype == np.float32
+ assert meso_nav.sun_za[0] == 0.4887488782405853
+ assert meso_nav.sun_az[0] == 3.976273775100708
# test with height correction
- abi_nav = navigation.ABINavigation(abi_data, precise_sun=False, hae_m=1.2345678)
+ meso_nav = navigation.ABINavigation(meso_data, precise_sun=False, hae_m=1.2345678)
- assert abi_nav.lat_deg.dtype == np.float32
- assert abi_nav.lon_deg.dtype == np.float32
- assert abi_nav.lat_deg[idx] == 44.73153305053711
- assert abi_nav.lon_deg[idx] == -93.01801300048828
+ assert meso_nav.lat_deg.dtype == np.float32
+ assert meso_nav.lon_deg.dtype == np.float32
+ assert meso_nav.lat_deg[idx] == 44.73153305053711
+ assert meso_nav.lon_deg[idx] == -93.01801300048828
- assert abi_nav.area_m.dtype == np.float32
- assert abi_nav.area_m[idx] == 4398248.0
+ assert meso_nav.area_m.dtype == np.float32
+ assert meso_nav.area_m[idx] == 4398248.0
- assert abi_nav.sat_za.dtype == np.float32
- assert abi_nav.sat_az.dtype == np.float32
- assert abi_nav.sat_za[idx] == 0.9509882926940918
- assert abi_nav.sat_az[idx] == 2.7129032611846924
+ assert meso_nav.sat_za.dtype == np.float32
+ assert meso_nav.sat_az.dtype == np.float32
+ assert meso_nav.sat_za[idx] == 0.9509882926940918
+ assert meso_nav.sat_az[idx] == 2.7129032611846924
- assert abi_nav.sun_za.dtype == np.float32
- assert abi_nav.sun_az.dtype == np.float32
- assert abi_nav.sun_za[idx] == 0.4886806607246399
- assert abi_nav.sun_az[idx] == 3.976076126098633
+ assert meso_nav.sun_za.dtype == np.float32
+ assert meso_nav.sun_az.dtype == np.float32
+ assert meso_nav.sun_za[idx] == 0.4886806607246399
+ assert meso_nav.sun_az[idx] == 3.976076126098633
# test with astropy sun and height correction
- abi_nav = navigation.ABINavigation(
- abi_data, index=idx, precise_sun=True, hae_m=1.2345678
+ meso_nav = navigation.ABINavigation(
+ meso_data, index=idx, precise_sun=True, hae_m=1.2345678
)
- assert abi_nav.sun_za.dtype == np.float32
- assert abi_nav.sun_az.dtype == np.float32
- assert abi_nav.sun_za[0] == 0.48874911665916443
- assert abi_nav.sun_az[0] == 3.9762721061706543
+ assert meso_nav.sun_za.dtype == np.float32
+ assert meso_nav.sun_az.dtype == np.float32
+ assert meso_nav.sun_za[0] == 0.48874911665916443
+ assert meso_nav.sun_az[0] == 3.9762721061706543
- # test on a navigation dataset that does not contain NaNs (G16 meso)
- abi_nav = navigation.ABINavigation(
- abi_data, lat_deg=44.72609499, lon_deg=-93.02279070
+ # test single-point indexing on a navigation dataset that does not contain NaNs (G16 meso)
+ meso_nav = navigation.ABINavigation(
+ meso_data, lat_bounds=44.72609499, lon_bounds=-93.02279070
)
- assert abi_nav.index == idx
+ assert meso_nav.index == idx
- # test retrieving indices fron multiple points
- abi_nav = navigation.ABINavigation(
- abi_data,
- lat_deg=np.array(
- [
- [44.73149490356445, 44.73029708862305],
- [44.70037841796875, 44.699180603027344],
- ],
- dtype=np.float32,
- ),
- lon_deg=np.array(
- [
- [-93.01798248291016, -92.98883819580078],
- [-93.00658416748047, -92.97746276855469],
- ],
- dtype=np.float32,
- ),
+ # test single-point indexing on a navigation dataset that contains NaNs (G16 CONUS)
+ conus_data = load(abi_cc07_nc)
+ conus_nav = navigation.ABINavigation(
+ conus_data, lat_bounds=44.72609499, lon_bounds=-93.02279070
)
- assert abi_nav.index == np.s_[92:94, 42:44]
+ assert conus_nav.index == (192, 1152)
- # test that the indices work
- abi_nav2 = navigation.ABINavigation(
- abi_data,
- index=abi_nav.index,
+
+def test_range_navigation():
+ # test meso
+ meso_data = load(abi_mc07_nc)
+
+ # test retrieving points from lat/lon bounds
+ lat_bounds = np.array(
+ [
+ 44.731495,
+ 44.602596,
+ ],
+ dtype=np.float32,
+ )
+ lon_bounds = np.array(
+ [
+ -93.01798,
+ -92.85648,
+ ],
+ dtype=np.float32,
+ )
+ meso_nav_bounded = navigation.ABINavigation(
+ meso_data, lat_bounds=lat_bounds, lon_bounds=lon_bounds, degrees=True
)
- assert (abi_nav2.lat_deg == abi_nav.lat_deg).all()
- assert (abi_nav2.lon_deg == abi_nav.lon_deg).all()
- # test on a navigation dataset that contains NaNs (G16 CONUS)
- abi_data = load(abi_cc07_nc)
- abi_nav = navigation.ABINavigation(
- abi_data, lat_deg=44.72609499, lon_deg=-93.02279070
+ # check index
+ assert meso_nav_bounded.index == np.s_[92:97, 42:47]
+
+ # check bounds
+ assert meso_nav_bounded.lat_deg[0, 0], (
+ meso_nav_bounded.lat_deg[-1, -1] == lat_bounds
+ )
+ assert meso_nav_bounded.lon_deg[0, 0], (
+ meso_nav_bounded.lon_deg[-1, -1] == lon_bounds
)
- assert abi_nav.index == (192, 1152)
- # test retrieving indices fron multiple points
- abi_nav = navigation.ABINavigation(
- abi_data,
- lat_deg=np.array(
- [
- [44.73149871826172, 44.73030090332031],
- [44.700382232666016, 44.699188232421875],
- ],
- dtype=np.float32,
- ),
- lon_deg=np.array(
- [
- [-93.01798248291016, -92.98883819580078],
- [-93.006591796875, -92.97746276855469],
- ],
- dtype=np.float32,
- ),
+ # another nav object using the same derived index should have the same gridded data
+ meso_nav_indexed = navigation.ABINavigation(
+ meso_data, index=meso_nav_bounded.index, degrees=True
)
- assert abi_nav.index == np.s_[192:194, 1152:1154]
+ for i in ["lat_deg", "lon_deg", "sat_za", "sat_az", "sun_za", "sun_az", "area_m"]:
+ assert (getattr(meso_nav_indexed, i) == getattr(meso_nav_bounded, i)).all()
- # test that the indices work
- abi_nav2 = navigation.ABINavigation(
- abi_data,
- index=abi_nav.index,
+ # test conus
+ conus_data = load(abi_cc07_nc)
+
+ # test retrieving points fron multiple lat/lon
+ lat_bounds = np.array(
+ [
+ [44.7315, 44.7303, 44.729107, 44.727917, 44.726727],
+ [44.700382, 44.69919, 44.697994, 44.696804, 44.695618],
+ [44.669365, 44.668175, 44.666985, 44.6658, 44.66461],
+ [44.638294, 44.637104, 44.635918, 44.63473, 44.63355],
+ [44.60733, 44.606144, 44.60496, 44.60378, 44.602596],
+ ],
+ dtype=np.float32,
+ )
+ lon_bounds = np.array(
+ [
+ [-93.01798, -92.98884, -92.9597, -92.93057, -92.90145],
+ [-93.00659, -92.97746, -92.94835, -92.91924, -92.890144],
+ [-92.99527, -92.96617, -92.93707, -92.90799, -92.878914],
+ [-92.98393, -92.95485, -92.92577, -92.89671, -92.86766],
+ [-92.97267, -92.94361, -92.91456, -92.88552, -92.85648],
+ ],
+ dtype=np.float32,
+ )
+ conus_nav_bounded = navigation.ABINavigation(
+ conus_data, lat_bounds=lat_bounds, lon_bounds=lon_bounds, degrees=True
+ )
+
+ # check index
+ assert conus_nav_bounded.index == np.s_[192:197, 1152:1157]
+
+ # check bounds
+ assert (conus_nav_bounded.lat_deg == lat_bounds).all()
+ assert (conus_nav_bounded.lon_deg == lon_bounds).all()
+
+ # another nav object using the same derived index should have the same gridded data
+ conus_nav_indexed = navigation.ABINavigation(
+ conus_data, index=conus_nav_bounded.index, degrees=True
)
- assert (abi_nav2.lat_deg == abi_nav.lat_deg).all()
- assert (abi_nav2.lon_deg == abi_nav.lon_deg).all()
+ for i in ["lat_deg", "lon_deg", "sat_za", "sat_az", "sun_za", "sun_az", "area_m"]:
+ assert (getattr(conus_nav_indexed, i) == getattr(conus_nav_bounded, i)).all()
+
+ # meso and conus nav reference the same area, the time-invariant gridded data should match within tolerance
+ assert (
+ conus_nav_bounded.lat_deg - meso_nav_bounded.lat_deg < epsilon_degrees
+ ).all()
+ assert (
+ conus_nav_bounded.lon_deg - meso_nav_bounded.lon_deg < epsilon_degrees
+ ).all()
+ assert (conus_nav_bounded.sat_za - meso_nav_bounded.sat_za < epsilon_degrees).all()
+ assert (conus_nav_bounded.sat_az - meso_nav_bounded.sat_az < epsilon_degrees).all()
+ assert (conus_nav_bounded.area_m - meso_nav_bounded.area_m < epsilon_meters).all()
def test_ancillary():
slc = np.s_[250:1250, 500:2000]
-
abi_data = load(abi_cc07_nc)
- water = ancillary.WaterMask(abi_data, index=slc, rivers=True)
- assert water.data["water_mask"].dtype == bool
- assert water.data["water_mask"].shape == (
- 1000,
- 1500,
+ # test water mask
+ water_indexed = ancillary.WaterMask(abi_data, index=slc, rivers=True)
+ water_bounded = ancillary.WaterMask(
+ abi_data,
+ lat_bounds=np.array([44.224377, 19.44], dtype=np.float32),
+ lon_bounds=np.array([-113.55036, -71.337036], dtype=np.float32),
+ rivers=True,
)
- water.save(save_dir=output_dir)
- cv2.imwrite(str(output_dir.joinpath("water.png")), water.data["water_mask"] * 255)
+ assert (
+ water_indexed.data["water_mask"].dtype
+ == water_bounded.data["water_mask"].dtype
+ == bool
+ )
+ assert (
+ water_indexed.data["water_mask"].shape
+ == water_bounded.data["water_mask"].shape
+ == (
+ 1000,
+ 1500,
+ )
+ )
+ assert (water_indexed.data["water_mask"] == water_bounded.data["water_mask"]).all()
- iremis = ancillary.IREMIS(abi_data, index=slc)
+ water_bounded.save(save_dir=output_dir)
+ cv2.imwrite(
+ str(output_dir.joinpath("water.png")), water_bounded.data["water_mask"] * 255
+ )
- assert iremis.data["c07_land_emissivity"].dtype == np.float32
- assert iremis.data["c07_land_emissivity"].shape == (
- 1000,
- 1500,
+ # test IREMIS
+ iremis_indexed = ancillary.IREMIS(abi_data, index=slc)
+ iremis_bounded = ancillary.IREMIS(
+ abi_data,
+ lat_bounds=np.array([44.224377, 19.44], dtype=np.float32),
+ lon_bounds=np.array([-113.55036, -71.337036], dtype=np.float32),
)
- iremis.save(save_dir=output_dir)
+ assert (
+ iremis_indexed.data["c07_land_emissivity"].dtype
+ == iremis_bounded.data["c07_land_emissivity"].dtype
+ == np.float32
+ )
+ assert (
+ iremis_indexed.data["c07_land_emissivity"].shape
+ == iremis_bounded.data["c07_land_emissivity"].shape
+ == (
+ 1000,
+ 1500,
+ )
+ )
+ assert (
+ iremis_indexed.data["c07_land_emissivity"]
+ == iremis_bounded.data["c07_land_emissivity"]
+ ).all()
+ assert (
+ iremis_indexed.data["c14_land_emissivity"]
+ == iremis_bounded.data["c14_land_emissivity"]
+ ).all()
+
+ iremis_bounded.save(save_dir=output_dir)
cv2.imwrite(
str(output_dir.joinpath("iremis.png")),
- minmax(iremis.data["c07_land_emissivity"]) * 255,
+ minmax(iremis_bounded.data["c07_land_emissivity"]) * 255,
)
diff --git a/test/test_resources.py b/test/test_resources.py
new file mode 100644
index 0000000..6dd8648
--- /dev/null
+++ b/test/test_resources.py
@@ -0,0 +1,121 @@
+from pathlib import Path
+
+SCRIPT_PATH = Path(__file__).parent.resolve()
+
+input_dir = SCRIPT_PATH.joinpath("input")
+input_dir.mkdir(parents=True, exist_ok=True)
+
+# abi
+abi_mc01_nc = input_dir.joinpath(
+ "abi/OR_ABI-L1b-RadM1-M6C01_G16_s20211691942252_e20211691942310_c20211691942342.nc"
+)
+abi_mc02_nc = input_dir.joinpath(
+ "abi/OR_ABI-L1b-RadM1-M6C02_G16_s20211691942252_e20211691942310_c20211691942334.nc"
+)
+abi_mc03_nc = input_dir.joinpath(
+ "abi/OR_ABI-L1b-RadM1-M6C03_G16_s20211691942252_e20211691942310_c20211691942351.nc"
+)
+abi_mc04_nc = input_dir.joinpath(
+ "abi/OR_ABI-L1b-RadM1-M6C04_G16_s20211691942252_e20211691942310_c20211691942340.nc"
+)
+abi_mc05_nc = input_dir.joinpath(
+ "abi/OR_ABI-L1b-RadM1-M6C05_G16_s20211691942252_e20211691942310_c20211691942347.nc"
+)
+abi_mc06_nc = input_dir.joinpath(
+ "abi/OR_ABI-L1b-RadM1-M6C06_G16_s20211691942252_e20211691942315_c20211691942345.nc"
+)
+abi_mc07_nc = input_dir.joinpath(
+ "abi/OR_ABI-L1b-RadM1-M6C07_G16_s20211691942252_e20211691942321_c20211691942355.nc"
+)
+abi_mc08_nc = input_dir.joinpath(
+ "abi/OR_ABI-L1b-RadM1-M6C08_G16_s20211691942252_e20211691942310_c20211691942357.nc"
+)
+abi_mc09_nc = input_dir.joinpath(
+ "abi/OR_ABI-L1b-RadM1-M6C09_G16_s20211691942252_e20211691942315_c20211691942368.nc"
+)
+abi_mc10_nc = input_dir.joinpath(
+ "abi/OR_ABI-L1b-RadM1-M6C10_G16_s20211691942252_e20211691942322_c20211691942353.nc"
+)
+abi_mc11_nc = input_dir.joinpath(
+ "abi/OR_ABI-L1b-RadM1-M6C11_G16_s20211691942252_e20211691942310_c20211691942348.nc"
+)
+abi_mc12_nc = input_dir.joinpath(
+ "abi/OR_ABI-L1b-RadM1-M6C12_G16_s20211691942252_e20211691942316_c20211691942356.nc"
+)
+abi_mc13_nc = input_dir.joinpath(
+ "abi/OR_ABI-L1b-RadM1-M6C13_G16_s20211691942252_e20211691942322_c20211691942361.nc"
+)
+abi_mc14_nc = input_dir.joinpath(
+ "abi/OR_ABI-L1b-RadM1-M6C14_G16_s20211691942252_e20211691942310_c20211691942364.nc"
+)
+abi_mc15_nc = input_dir.joinpath(
+ "abi/OR_ABI-L1b-RadM1-M6C15_G16_s20211691942252_e20211691942316_c20211691942358.nc"
+)
+abi_mc16_nc = input_dir.joinpath(
+ "abi/OR_ABI-L1b-RadM1-M6C16_G16_s20211691942252_e20211691942322_c20211691942366.nc"
+)
+
+# add a few g16 conuses to test off-earth pixels
+abi_cc01_nc = input_dir.joinpath(
+ "abi/OR_ABI-L1b-RadC-M6C01_G16_s20211691941174_e20211691943547_c20211691943589.nc"
+)
+abi_cc02_nc = input_dir.joinpath(
+ "abi/OR_ABI-L1b-RadC-M6C02_G16_s20211691941174_e20211691943547_c20211691943571.nc"
+)
+abi_cc03_nc = input_dir.joinpath(
+ "abi/OR_ABI-L1b-RadC-M6C03_G16_s20211691941174_e20211691943547_c20211691943587.nc"
+)
+abi_cc07_nc = input_dir.joinpath(
+ "abi/OR_ABI-L1b-RadC-M6C07_G16_s20211691941174_e20211691943558_c20211691944002.nc"
+)
+
+abi_ncs = [
+ abi_cc01_nc,
+ abi_cc02_nc,
+ abi_cc03_nc,
+ abi_cc07_nc,
+ abi_mc01_nc,
+ abi_mc02_nc,
+ abi_mc03_nc,
+ abi_mc04_nc,
+ abi_mc05_nc,
+ abi_mc06_nc,
+ abi_mc07_nc,
+ abi_mc08_nc,
+ abi_mc09_nc,
+ abi_mc10_nc,
+ abi_mc11_nc,
+ abi_mc12_nc,
+ abi_mc13_nc,
+ abi_mc14_nc,
+ abi_mc15_nc,
+ abi_mc16_nc,
+]
+
+# suvi
+suvi_094_nc = input_dir.joinpath(
+ "suvi/OR_SUVI-L1b-Fe093_G16_s20203160623501_e20203160623511_c20203160624091.nc"
+)
+suvi_131_nc = input_dir.joinpath(
+ "suvi/OR_SUVI-L1b-Fe131_G16_s20203160623001_e20203160623011_c20203160623196.nc"
+)
+suvi_171_nc = input_dir.joinpath(
+ "suvi/OR_SUVI-L1b-Fe171_G16_s20203160624201_e20203160624211_c20203160624396.nc"
+)
+suvi_195_nc = input_dir.joinpath(
+ "suvi/OR_SUVI-L1b-Fe195_G16_s20203160623301_e20203160623311_c20203160623491.nc"
+)
+suvi_284_nc = input_dir.joinpath(
+ "suvi/OR_SUVI-L1b-Fe284_G16_s20203160624501_e20203160624511_c20203160625090.nc"
+)
+suvi_304_nc = input_dir.joinpath(
+ "suvi/OR_SUVI-L1b-He303_G16_s20203160622501_e20203160622511_c20203160623090.nc"
+)
+suvi_ncs = [
+ suvi_094_nc,
+ suvi_131_nc,
+ suvi_171_nc,
+ suvi_195_nc,
+ suvi_284_nc,
+ suvi_304_nc,
+]