Skip to content

Commit

Permalink
ABI imagery and NCInterface improvements
Browse files Browse the repository at this point in the history 
- Reduce memory usage by lazy-loading slices of ABI images when desired
- Reveal the original mask of an ABI image with the `.mask` attribute
- By user request: set the space background of ABI imagery to be black with the `black_space` boolean argument in ABIImage and ABINaturalRGB
- Improve consistency of NCInterface variable behavior when when indices are changed
- Syntax updates
  • Loading branch information
@wx-star
wx-star committed Nov 5, 2023
1 parent ecdc3ad commit 2d7fdd4
Show file tree
Hide file tree
Showing 8 changed files with 395 additions and 256 deletions.
2 changes: 1 addition & 1 deletion README.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -79,7 +79,7 @@ With an L1b SUVI netCDF file `suvi_nc`, view the exposure time in seconds:
```python
from heregoes import load
suvi_data = load(suvi_nc)
print(suvi_data["CMD_EXP"][:])
print(suvi_data["CMD_EXP"][...])
```

Which returns:
Expand Down
12 changes: 6 additions & 6 deletions heregoes/ancillary.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -109,7 +109,7 @@ class IREMIS(AncillaryDataset):
"""

def __init__(self, abi_data, iremis_dir=IREMIS_DIR):
super(IREMIS, self).__init__()
super().__init__()

self.abi_data = abi_data
month = self.abi_data.time_coverage_start.month
Expand Down Expand Up @@ -152,15 +152,15 @@ def __init__(self, abi_data, iremis_dir=IREMIS_DIR):
self.data["c07_land_emissivity"] = linear_interp(
3.7,
4.3,
iremis["emis1"][:],
iremis["emis2"][:],
iremis["emis1"][...],
iremis["emis2"][...],
3.9,
).astype(np.float32)
self.data["c14_land_emissivity"] = linear_interp(
10.8,
12.1,
iremis["emis8"][:],
iremis["emis9"][:],
iremis["emis8"][...],
iremis["emis9"][...],
11.2,
).astype(np.float32)

Expand Down Expand Up @@ -210,7 +210,7 @@ class WaterMask(AncillaryDataset):
"""

def __init__(self, abi_data, gshhs_scale="intermediate", rivers=False):
super(WaterMask, self).__init__()
super().__init__()

self.abi_data = abi_data
self.dataset_name = "gshhs_" + gshhs_scale
Expand Down
12 changes: 6 additions & 6 deletions heregoes/goesr/__init__.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -94,12 +94,12 @@ def __init__(self, *args, **kwargs):
self.resolution_ifov = 56.0e-6
self.resolution_km = 2.0

self.band_id_safe = "C" + str(self["band_id"][:].item()).zfill(2)
self.band_id_safe = "C" + str(self["band_id"][...].item()).zfill(2)

self.midpoint_time = self.epoch2timestamp(seconds=float(self["t"][:].item()))
self.midpoint_time = self.epoch2timestamp(seconds=float(self["t"][...].item()))

self.instrument_coefficients = coefficients.ABICoeff(
self.platform_ID, self["band_id"][:].item()
self.platform_ID, self["band_id"][...].item()
)


Expand All @@ -111,11 +111,11 @@ def __init__(self, *args, **kwargs):

# Wavelength for SUVI 304 is masked in netCDF
self["WAVELNTH"].set_fill_value(0)
self.wavelength_safe = str(int(self["WAVELNTH"][:].item())).zfill(3)
self.wavelength_safe = str(int(self["WAVELNTH"][...].item())).zfill(3)
if self.wavelength_safe == "000":
self.wavelength_safe = "304"
self["WAVELNTH"][:] = 304
self["WAVELNTH"][...] = 304

self.instrument_coefficients = coefficients.SUVICoeff(
self.platform_ID, self["WAVELNTH"][:].item()
self.platform_ID, self["WAVELNTH"][...].item()
)
100 changes: 70 additions & 30 deletions heregoes/image.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -85,6 +85,7 @@ def __init__(
abi_nc,
index=slice(None, None),
gamma=1.0,
black_space=False,
):
"""
Creates Cloud Moisture Imagery (CMI) following the CMIP ATBD: https://www.star.nesdis.noaa.gov/goesr/docs/ATBD/Imagery.pdf
Expand All @@ -102,16 +103,21 @@ def __init__(
- `abi_nc`: String or Path object pointing to a GOES-R ABI L1b Radiance netCDF file
- `index`: Optionally process an ABI image for a single array index or slice
- `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`
"""

super(ABIImage, self).__init__()
super().__init__()
self.gamma = gamma
self.black_space = black_space

self._cmi = None

self.abi_data = load(abi_nc)

self.rad = self.abi_data["Rad"][index]
self.dqf = self.abi_data["DQF"][index]
self.quality = self.abi_data["Rad"].quality
self.mask = self.abi_data["Rad"].mask
self.quality = self.abi_data["Rad"].pct_unmasked

self.rad_range = np.array(
self.abi_data["Rad"].valid_range * self.abi_data["Rad"].scale_factor
Expand All @@ -132,20 +138,20 @@ def __init__(
@property
def cmi(self):
if self._cmi is None:
if 1 <= self.abi_data["band_id"][:] <= 6:
if 1 <= self.abi_data["band_id"][...] <= 6:
self._cmi = abi.rad2rf(
self.rad,
self.abi_data["earth_sun_distance_anomaly_in_AU"][:].item(),
self.abi_data["esun"][:].item(),
self.abi_data["earth_sun_distance_anomaly_in_AU"][...].item(),
self.abi_data["esun"][...].item(),
)

elif 7 <= self.abi_data["band_id"][:] <= 16:
elif 7 <= self.abi_data["band_id"][...] <= 16:
self._cmi = abi.rad2bt(
self.rad,
self.abi_data["planck_fk1"][:].item(),
self.abi_data["planck_fk2"][:].item(),
self.abi_data["planck_bc1"][:].item(),
self.abi_data["planck_bc2"][:].item(),
self.abi_data["planck_fk1"][...].item(),
self.abi_data["planck_fk2"][...].item(),
self.abi_data["planck_bc1"][...].item(),
self.abi_data["planck_bc2"][...].item(),
)

return self._cmi
Expand All @@ -157,20 +163,23 @@ def cmi(self, value):
@property
def bv(self):
if self._bv is None:
if 1 <= self.abi_data["band_id"][:] <= 6:
if 1 <= self.abi_data["band_id"][...] <= 6:
# calculate the range of possible reflectance factors from the provided valid range of radiance, and use it to normalize before the gamma correction
self.rf_min, self.rf_max = (
self.rad_range
* np.pi
* np.square(self.abi_data["earth_sun_distance_anomaly_in_AU"][:])
) / self.abi_data["esun"][:]
* np.square(self.abi_data["earth_sun_distance_anomaly_in_AU"][...])
) / self.abi_data["esun"][...]
self._bv = abi.rf2bv(
self.cmi, min=self.rf_min, max=self.rf_max, gamma=self.gamma
)

elif 7 <= self.abi_data["band_id"][:] <= 16:
elif 7 <= self.abi_data["band_id"][...] <= 16:
self._bv = abi.bt2bv(self.cmi)

if self.black_space:
self._bv[self.mask] = 0

return self._bv

@bv.setter
Expand All @@ -190,6 +199,7 @@ def __init__(
upscale=False,
upscale_algo=cv2.INTER_CUBIC,
gamma=1.0,
black_space=False,
):
"""
Creates the "natural" color RGB for ABI following https://doi.org/10.1029/2018EA000379 in BGR order
Expand All @@ -200,13 +210,14 @@ def __init__(
- `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`
- `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`
"""

super(ABINaturalRGB, self).__init__()
super().__init__()

red_image = ABIImage(red_nc, gamma=gamma)
green_image = ABIImage(green_nc, gamma=gamma)
blue_image = ABIImage(blue_nc, gamma=gamma)
red_image = ABIImage(red_nc, gamma=gamma, black_space=black_space)
green_image = ABIImage(green_nc, gamma=gamma, black_space=black_space)
blue_image = ABIImage(blue_nc, gamma=gamma, black_space=black_space)

if upscale:
# upscale green and blue to the size of red
Expand All @@ -216,20 +227,37 @@ def __init__(
interpolation=upscale_algo,
)
green_image.dqf = cv2.resize(
green_image.bv,
green_image.dqf,
(red_image.dqf.shape[1], red_image.dqf.shape[0]),
interpolation=cv2.INTER_NEAREST,
)
green_image.mask = (
cv2.resize(
green_image.mask.astype(np.uint8),
(red_image.mask.shape[1], red_image.mask.shape[0]),
interpolation=cv2.INTER_NEAREST,
)
== 1
)

blue_image.bv = cv2.resize(
blue_image.bv,
(red_image.bv.shape[1], red_image.bv.shape[0]),
interpolation=upscale_algo,
)
blue_image.dqf = cv2.resize(
blue_image.bv,
blue_image.dqf,
(red_image.dqf.shape[1], red_image.dqf.shape[0]),
interpolation=cv2.INTER_NEAREST,
)
blue_image.mask = (
cv2.resize(
blue_image.mask.astype(np.uint8),
(red_image.mask.shape[1], red_image.mask.shape[0]),
interpolation=cv2.INTER_NEAREST,
)
== 1
)

else:
# downscale red to the size of green and blue
Expand All @@ -239,10 +267,18 @@ def __init__(
interpolation=cv2.INTER_AREA,
)
red_image.dqf = cv2.resize(
red_image.bv,
red_image.dqf,
(green_image.dqf.shape[1], green_image.dqf.shape[0]),
interpolation=cv2.INTER_NEAREST,
)
red_image.mask = (
cv2.resize(
red_image.mask.astype(np.uint8),
(green_image.mask.shape[1], green_image.mask.shape[0]),
interpolation=cv2.INTER_NEAREST,
)
== 1
)

green_image.bv = (
(red_image.bv * r_coeff)
Expand All @@ -257,14 +293,17 @@ def __init__(
sum([red_image.quality, green_image.quality, blue_image.quality]) / 3
)
self.dqf = np.stack([blue_image.dqf, green_image.dqf, red_image.dqf], axis=2)
self.mask = np.stack(
[blue_image.mask, green_image.mask, red_image.mask], axis=2
)

if upscale:
self.abi_data = red_image.abi_data

else:
self.abi_data = green_image.abi_data

self.abi_data["band_id"][:] = np.atleast_1d(0)
self.abi_data["band_id"][...] = np.atleast_1d(0)
self.abi_data.band_id_safe = "Color"

self.abi_data.dataset_name = "RGB from " + ", ".join(
Expand Down Expand Up @@ -303,20 +342,21 @@ def __init__(
- `dqf_correction`: Whether to interpolate over bad pixels marked by DQF. Default `True`
"""

super(SUVIImage, self).__init__()
super().__init__()

self.suvi_data = load(suvi_nc)

if self.suvi_data["CMD_EXP"][:] != 1.0:
if self.suvi_data["CMD_EXP"][...] != 1.0:
logger.warning(
"Short SUVI exposure detected: SUVI exposures shorter than 1 second are not officially supported."
"Short SUVI exposure detected: SUVI exposures shorter than 1 second are not officially supported.",
extra={"caller": f"{__name__}.{self.__class__.__name__}"},
)

self.rad = self.suvi_data["RAD"][:]
self.dqf = self.suvi_data["DQF"][:]
self.quality = self.suvi_data["RAD"].quality
x_offset = 640 - self.suvi_data["CRPIX1"][:]
y_offset = 640 - self.suvi_data["CRPIX2"][:]
self.rad = self.suvi_data["RAD"][...]
self.dqf = self.suvi_data["DQF"][...]
self.quality = self.suvi_data["RAD"].pct_unmasked
x_offset = 640 - self.suvi_data["CRPIX1"][...]
y_offset = 640 - self.suvi_data["CRPIX2"][...]

self.default_filename = "_".join(
(
Expand Down
10 changes: 5 additions & 5 deletions heregoes/navigation.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -76,8 +76,8 @@ def __init__(

if self.index == slice(None, None):
self.x_rad, self.y_rad = np.meshgrid(
self.abi_data["x"][:],
self.abi_data["y"][:],
self.abi_data["x"][...],
self.abi_data["y"][...],
)

else:
Expand Down Expand Up @@ -218,9 +218,9 @@ def area_m(self, value):

def _calc_sat(self):
self._sat_az, self._sat_za = orbital.get_observer_look(
sat_lon=np.atleast_1d(self.abi_data["nominal_satellite_subpoint_lon"][:]),
sat_lat=np.atleast_1d(self.abi_data["nominal_satellite_subpoint_lat"][:]),
sat_alt=np.atleast_1d(self.abi_data["nominal_satellite_height"][:]),
sat_lon=np.atleast_1d(self.abi_data["nominal_satellite_subpoint_lon"][...]),
sat_lat=np.atleast_1d(self.abi_data["nominal_satellite_subpoint_lat"][...]),
sat_alt=np.atleast_1d(self.abi_data["nominal_satellite_height"][...]),
jdays2000=orbital.jdays2000(self.time),
lon=self.lon_deg,
lat=self.lat_deg,
Expand Down
8 changes: 4 additions & 4 deletions heregoes/projection.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -56,10 +56,10 @@ def __init__(self, abi_data):
lon_0 = self.abi_data["goes_imager_projection"].longitude_of_projection_origin
sweep = self.abi_data["goes_imager_projection"].sweep_angle_axis

ul_x = self.abi_data["x_image_bounds"][0] * h
ul_y = self.abi_data["y_image_bounds"][0] * h
lr_x = self.abi_data["x_image_bounds"][1] * h
lr_y = self.abi_data["y_image_bounds"][1] * h
ul_x = (self.abi_data["x_image_bounds"][0] * h).item()
ul_y = (self.abi_data["y_image_bounds"][0] * h).item()
lr_x = (self.abi_data["x_image_bounds"][1] * h).item()
lr_y = (self.abi_data["y_image_bounds"][1] * h).item()
self.abi_bounds = [ul_x, ul_y, lr_x, lr_y]

self._intermediate_format = "GTiff"
Expand Down
Loading

0 comments on commit 2d7fdd4

Please sign in to comment.