# # # Distribution Statement A. Approved for public release. Distribution unlimited.
# # #
# # # Author:
# # # Naval Research Laboratory, Marine Meteorology Division
# # #
# # # This program is free software: you can redistribute it and/or modify it under
# # # the terms of the NRLMMD License included with this program. This program is
# # # distributed WITHOUT ANY WARRANTY; without even the implied warranty of
# # # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the included license
# # # for more details. If you did not receive the license, for more information see:
# # # https://github.com/U-S-NRL-Marine-Meteorology-Division/
"""Reader to read a grannual NOAA ATMS SDR TBs in h5 format.
Output variables in xarray object for geoips processing system
V0: August 25, 2021
The date is generated by the NOAA community satellite processing package
(CSPP), developed at CIMSS
Example of ATMS file names::
'SATMS_j01_d20210809_t0959306_e1000023_b19296_fnmoc_ops.h5'
SDR TBs variables
'GATMO_j01_d20210809_t0959306_e1000023_b19296_fnmoc_ops.h5'
SDR Geolocation variables
Dataset info::
TB[12,96,22]: for each granuel
CHAN# Center-Freq(GHz) POL
1 23.8 V
2 31.4 V
3 50.3 H
4 51.76 H
5 52.8 H
6 53.596+-0.115 H
7 54.4 H
8 54.94 H
9 55.5 H
10 57.290(f0) H
11 f0 +-0.217 H
12 f0 +-0.322+-0.048 H
13 f0 +-0.322+-0.022 H
14 f0 +-0.322+-0.010 H
15 f0 +-0.322+-0.0045 H
16 88.2 V
17 165.5 H
18 183.1+-7 H
19 183.1+-4.5 H (FNMOC used this chan for 183 GHz image)
20 183.1+-3.0 H
21 183.1+-1.8 H
22 183.1+-1.0 H
BeamTime[12,96]: microsecond, i.e., 1*10^-6. IET "IDPS Epoch Time" is used
It is a signed 64-bit integer giving microseconds since
00:00:00.000000 Jan 1 1958.
BrightnessTemperatureFactors[2]: 1: scale (unitless); 2: offset (K)
BrightnessTemperature[12,96,22]: [scan,pix,chans]
SDR geolocation Info
Latitude/Longitude[12,96]: for Chan 17 only (for initial product,
it is used for all channels)
BeamLatitude[12,96,5]: for chan 1,2,3,16,17. They will be used for
associated TBs at a later date.
BeamLongitude[12,96,5]: for chan 1,2,3,16,17
SatelliteAzimuthAngle, SatelliteZenithAngle,
SolarAzimuthAngle, SolarZenithAngle[12,96]
Notes
-----
Unix epoch time is defined as the number of seconds that have elapsed
since January 1, 1970 (midnight UTC/GMT). Thus, there is a 12 years
difference for the JPSS data when
datetime.datetime.utcfromtimestamp(epoch) is used
to convert the JPSS IDPS Epoch time to the humman-readable date.
This reader is developed to read one granual a time from ATMS npp and
jpss-1(n20) data files.
The example files are:
* ``SATMS_j01_d20210809_t0959306_e1000023_b19296_fnmoc_ops.h5``: for TBs. 'b': orbit#
* ``GATMO_j01_d20210809_t0959306_e1000023_b19296_fnmoc_ops.h5``: for geolocations
"""
# Python Standard Libraries
from os.path import basename
import h5py
import numpy as np
import datetime
import xarray as xr
from astropy.time import Time
from dateutil.relativedelta import relativedelta
import logging
LOG = logging.getLogger(__name__)
interface = "readers"
family = "standard"
name = "atms_hdf5"
# from IPython import embed as shell
# list of variables selected from input files
atms_vars = [
"BrightnessTemperature",
"BrightnessTemperatureFactors",
"BeamTime",
"Latitude",
"Longitude",
"SatelliteAzimuthAngle",
"SatelliteZenithAngle",
"SolarAzimuthAngle",
"SolarZenithAngle",
"StartTime",
]
# unify common names for sat/sun-zenith and azimuth angles
xvarnames = {
"SolarZenithAngle": "solar_zenith_angle",
"SolarAzimuthAngle": "solar_azimuth_angle",
"SatelliteZenithAngle": "satellite_zenith_angle",
"SatelliteAzimuthAngle": "satellite_azimuth_angle",
}
final_xarray = xr.Dataset() # define a xarray to hold all selected variables
[docs]def convert_epoch_to_datetime64(time_array, use_shape=None):
"""Convert time to datetime object.
Parameters
----------
time_array : array
Array of start time integers (multiplied by 1e-6 in function)
use_shape : tuple, optional
desired output shape of time array, by default None
Returns
-------
array
array of converted datetime objects
"""
# Convert TAI to UTC using astropy
utc_time = Time(time_array / 1e6, format="unix_tai").utc.datetime
# ATMS Epoch starts at 1958-01-01, not 1970-01-01
utc_time -= relativedelta(years=12)
# Either convert 1D array to 2D, or return original shape
if time_array.ndim == 1 and use_shape:
nscan, npix = use_shape
converted_time = np.zeros((nscan, npix)).astype(datetime.datetime)
for i in range(nscan):
converted_time[i, :] = utc_time[i]
else:
converted_time = utc_time
return converted_time.astype(np.datetime64)
[docs]def read_atms_file(fname, xarray_atms):
"""Read ATCF data from file fname."""
fileobj = h5py.File(fname, mode="r")
# check for available variables from nput file
if "ATMS-SDR_All" in fileobj["All_Data"].keys(): # for TB-data
data_select = fileobj["All_Data"]["ATMS-SDR_All"]
tb = data_select["BrightnessTemperature"][()]
tb_time = data_select["BeamTime"][()]
tb_factor = data_select["BrightnessTemperatureFactors"][()]
# convert tb to actual values
tbs = tb * tb_factor[0] + tb_factor[1]
# TBs for selected channels
V23 = tbs[:, :, 0]
V31 = tbs[:, :, 1]
H50 = tbs[:, :, 2]
V89 = tbs[:, :, 15]
H165 = tbs[:, :, 16]
H183 = tbs[:, :, 18] # to match the 183+-4.5 GHz channel used by FNMOC
# get UTC time in datetime64 format required by geoips for each pixel
time_scan = convert_epoch_to_datetime64(tb_time)
# make dict of numpy arrays
ingested = {
"V23": V23,
"V31": V31,
"H50": H50,
"V89": V89,
"H165": H165,
"H183": H183,
"sdr_time": time_scan,
}
if "ATMS-SDR-GEO_All" in fileobj["All_Data"].keys(): # for geo-data
data_select = fileobj["All_Data"]["ATMS-SDR-GEO_All"]
lat = data_select["Latitude"][()]
lon = data_select["Longitude"][()]
solar_zenith_angle = data_select["SolarZenithAngle"][()]
solar_azimuth_angle = data_select["SolarAzimuthAngle"][()]
satellite_zenith_angle = data_select["SatelliteZenithAngle"][()]
satellite_azimuth_angle = data_select["SatelliteAzimuthAngle"][()]
StartTime = data_select["StartTime"][()]
# get UTC time in datetime64 format required by geoips for each pixel
time_scan = convert_epoch_to_datetime64(StartTime, use_shape=lon.shape)
# make dict of numpy arrays
ingested = {
"latitude": lat,
"longitude": lon,
"solar_zenith_angle": solar_zenith_angle,
"solar_azimuth_angle": solar_azimuth_angle,
"satellite_zenith_angle": satellite_zenith_angle,
"satellite_azimuth_angle": satellite_azimuth_angle,
"geo_time": time_scan,
}
# close the h5 object
fileobj.close()
# ------ setup xarray variables ------
for key, data in ingested.items():
if key not in xarray_atms.variables.keys():
final_xarray[key] = xr.DataArray(data)
else:
merged_array = np.vstack([xarray_atms[key], data])
final_xarray[key] = xr.DataArray(
merged_array, dims=["dim_" + str(merged_array.shape[0]), "dim_1"]
)
# Set official time to either sdr_time or geo_time
final_keys = final_xarray.variables.keys()
if "sdr_time" in final_keys:
final_xarray["time"] = final_xarray["sdr_time"]
else:
final_xarray["time"] = final_xarray["geo_time"]
return final_xarray
[docs]def call(fnames, metadata_only=False, chans=None, area_def=None, self_register=False):
"""Read ATMS hdf5 data products.
Parameters
----------
fnames : list
* List of strings, full paths to files
metadata_only : bool, default=False
* NOT YET IMPLEMENTED
* Return before actually reading data if True
chans : list of str, default=None
* NOT YET IMPLEMENTED
* List of desired channels (skip unneeded variables as needed).
* Include all channels if None.
area_def : pyresample.AreaDefinition, default=None
* NOT YET IMPLEMENTED
* Specify region to read
* Read all data if None.
self_register : str or bool, default=False
* NOT YET IMPLEMENTED
* register all data to the specified dataset id (as specified in the
return dictionary keys).
* Read multiple resolutions of data if False.
Returns
-------
dict of xarray.Datasets
* dictionary of xarray.Dataset objects with required Variables and
Attributes.
* Dictionary keys can be any descriptive dataset ids.
See Also
--------
:ref:`xarray_standards`
Additional information regarding required attributes and variables
for GeoIPS-formatted xarray Datasets.
"""
LOG.info("Reading files %s", fnames)
if metadata_only is True: # read-in datafiles first time if 'metadata_only= True'
xarray_atms = xr.Dataset()
source_file_names = []
for fname in fnames:
xarray_atms = read_atms_file(fname, xarray_atms)
source_file_names += [basename(fname)] # name of last file from input files
xarray_atms.attrs["source_file_names"] = source_file_names
# setup attributors
fileobj = h5py.File(fname, mode="r")
from geoips.xarray_utils.time import (
get_max_from_xarray_time,
get_min_from_xarray_time,
)
xarray_atms.attrs["start_datetime"] = get_min_from_xarray_time(
xarray_atms, "time"
)
xarray_atms.attrs["end_datetime"] = get_max_from_xarray_time(
xarray_atms, "time"
)
xarray_atms.attrs["source_name"] = "atms"
xarray_atms.attrs["platform_name"] = fileobj.attrs["Platform_Short_Name"][
0, 0
].decode(
"utf-8"
) # could be changed if needed
xarray_atms.attrs["data_provider"] = "NOAA"
# MTIFs need to be "prettier" for PMW products, so 2km resolution for
# final image
xarray_atms.attrs["sample_distance_km"] = 2
xarray_atms.attrs[
"interpolation_radius_of_influence"
] = 30000 # could be tuned if needed
fileobj.close()
else: # if 'metadata_only= False', it is for the second time to read-in datafiles
xarray_atms = final_xarray # avoid read the data the second time
return {"METADATA": xarray_atms[[]], "ATMS": xarray_atms}