Skip to content

Sprint: Updating two GOES and three marine converters for naming conventions #1121

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
srherbener merged 31 commits into from
Jan 3, 2023
Merged

Sprint: Updating two GOES and three marine converters for naming conventions #1121

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
31 commits
Select commit Hold shift + click to select a range
bf75378
update GOES converter to data conventions
gthompsnWRF Nov 14, 2022
c76dfad
updated some marine converters for dateTime and naming conventions
gthompsnWRF Nov 16, 2022
d22f8a1
updated MetaData group vars to new naming convention
gthompsnWRF Nov 20, 2022
5465b36
switch preqc from float to int
gthompsnWRF Dec 6, 2022
872e4e2
quick change of Location to 32-bit integer, not 64-bit
gthompsnWRF Dec 7, 2022
6651f4d
changing wording to albedo from reflectance factor
gthompsnWRF Dec 9, 2022
8225561
fix latlon Location dimension to int32 type, not int64
gthompsnWRF Dec 19, 2022
2470c66
fix latlon Location dimension to int32 type, not int64
gthompsnWRF Dec 19, 2022
f1783b7
fix indent for coding norm
gthompsnWRF Dec 20, 2022
69de9ed
another coding norm fix
gthompsnWRF Dec 20, 2022
c9e6041
import fix that was added to develop is copied here also
gthompsnWRF Dec 20, 2022
f5b4f9b
updating few marine converters to use dateTime properly
gthompsnWRF Dec 20, 2022
db36535
updating a couple test reference files
gthompsnWRF Dec 20, 2022
17d7f37
updating smap land converters and test reference files for naming con…
gthompsnWRF Dec 20, 2022
399d736
amsr2 converter and test reference file updated for naming conventions
gthompsnWRF Dec 20, 2022
7196ca8
more marine and land converters updated naming convention
gthompsnWRF Dec 21, 2022
fa370ee
updated omi O3 converter for naming conventions
gthompsnWRF Dec 22, 2022
ed0c1d4
update NSIDC ice fraction converter to naming conventions
gthompsnWRF Dec 22, 2022
6a86620
quick fixes for coding norms
gthompsnWRF Dec 22, 2022
fa0596b
updating mopitt CO converter to naming conventions
gthompsnWRF Dec 23, 2022
6a34890
fix MODIS AOD converter and adopt naming conventions
gthompsnWRF Dec 23, 2022
1db00df
update MLS ozone converter to naming conventions
gthompsnWRF Dec 23, 2022
6839b75
updating AVHRR converter to naming conventions
gthompsnWRF Dec 23, 2022
57af4ac
update 1 of 2 sig-wave-height converters to naming convention
gthompsnWRF Dec 23, 2022
e73a61b
update copernicus absoluteDynamicTopography convert to naming convent…
gthompsnWRF Dec 23, 2022
3771c01
fix a coding norm problem
gthompsnWRF Dec 23, 2022
f4c8a82
fix overindent coding norm
gthompsnWRF Dec 24, 2022
91613c8
change ncdiff tolerance check on some (to match develop branch)
gthompsnWRF Dec 24, 2022
afbca6c
update GNSSRO python converter to naming convention
gthompsnWRF Dec 24, 2022
b8947cf
updating testoutput file for change to fillvalue
gthompsnWRF Dec 24, 2022
778b1ba
empty commit
BenjaminRuston Dec 30, 2022
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
42 changes: 19 additions & 23 deletions src/compo/mls_o3_nc2ioda.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -29,36 +29,37 @@
locationKeyList = [
("latitude", "float"),
("longitude", "float"),
("air_pressure", "float"),
("dateTime", "integer"),
("pressure", "float"),
("dateTime", "long"),
]

varname_ozone = 'ozoneProfile'

ioda2nc = {}
ioda2nc['latitude'] = 'HDFEOS/SWATHS/O3/Geolocation Fields/Latitude'
ioda2nc['longitude'] = 'HDFEOS/SWATHS/O3/Geolocation Fields/Longitude'
ioda2nc['dateTime'] = 'HDFEOS/SWATHS/O3/Geolocation Fields/Time'
ioda2nc['air_pressure'] = 'HDFEOS/SWATHS/O3/Geolocation Fields/Pressure'
ioda2nc['pressure'] = 'HDFEOS/SWATHS/O3/Geolocation Fields/Pressure'
ioda2nc['valKey'] = 'HDFEOS/SWATHS/O3/Data Fields/O3'
ioda2nc['precision'] = 'HDFEOS/SWATHS/O3/Data Fields/O3Precision'
ioda2nc['convergence'] = 'HDFEOS/SWATHS/O3/Data Fields/Convergence'
ioda2nc['status'] = 'HDFEOS/SWATHS/O3/Data Fields/Status'
ioda2nc['quality'] = 'HDFEOS/SWATHS/O3/Data Fields/Quality'
ioda2nc['solar_zenith_angle'] = 'HDFEOS/SWATHS/O3/Geolocation Fields/SolarZenithAngle'
ioda2nc['solarZenithAngle'] = 'HDFEOS/SWATHS/O3/Geolocation Fields/SolarZenithAngle'

obsvars = {
'mole_fraction_of_ozone_in_air': 'mole_fraction_of_ozone_in_air',
'mole_fraction_of_ozone_in_air': varname_ozone,
}

AttrData = {
'converter': os.path.basename(__file__),
'nvars': np.int32(len(obsvars)),
'converter': os.path.basename(__file__)
}

DimDict = {
}

VarDims = {
'mole_fraction_of_ozone_in_air': ['nlocs'],
varname_ozone: ['Location'],
}


Expand All @@ -79,10 +80,10 @@ def __init__(self, filenames, lvmin, lvmax, sTAI, eTAI, nrt, qcOn, errorOn):
if(v != 'valKey' and v != 'errKey'):
self.outdata[(v, 'MetaData')] = []
self.outdata[('level', 'MetaData')] = []
self._setVarDict('mole_fraction_of_ozone_in_air')
self.outdata[self.varDict['mole_fraction_of_ozone_in_air']['valKey']] = []
self._setVarDict(varname_ozone)
self.outdata[self.varDict[varname_ozone]['valKey']] = []
if(self.qcOn):
self.outdata[self.varDict['mole_fraction_of_ozone_in_air']['errKey']] = []
self.outdata[self.varDict[varname_ozone]['errKey']] = []

self._read()

Expand Down Expand Up @@ -112,12 +113,8 @@ def _setVarAttr(self, iodavar):
self.varAttrs[vkey]['units'] = 'seconds since 1993-01-01T00:00:00Z'
elif('angle' in v.lower() or 'latitude' in v.lower() or 'longitude' in v.lower()):
self.varAttrs[vkey]['units'] = 'degrees'
elif('flag' in v.lower()):
self.varAttrs[vkey]['units'] = 'unitless'
elif('prior' in v.lower()):
self.varAttrs[vkey]['units'] = 'ppmv'
else:
self.varAttrs[vkey]['units'] = 'unitless'

# Read data needed from raw MLS file.
def _read_nc(self, filename, ifile, maxfile):
Expand All @@ -139,7 +136,7 @@ def _read_nc(self, filename, ifile, maxfile):

d = {}
for k in list(ioda2nc.keys()):
if (k == 'air_pressure'):
if (k == 'pressure'):
d[k] = ncd[ioda2nc[k]][...]*100. # convert to Pa
d[k].mask = False
else:
Expand Down Expand Up @@ -179,13 +176,13 @@ def _just_flatten(self, d):
dd['precision'] = d['precision'][idx, self.lmin:self.lmax+1]
lvec = np.arange(self.lmin+1, self.lmax+2)
dd['level'], dd['status'] = np.meshgrid(np.arange(self.lmin+1, self.lmax+2), d['status'][idx])
dd['air_pressure'], dd['dateTime'] = np.meshgrid(d['air_pressure'][self.lmin:self.lmax+1], d['dateTime'][idx])
dd['pressure'], dd['dateTime'] = np.meshgrid(d['pressure'][self.lmin:self.lmax+1], d['dateTime'][idx])
dd['quality'] = np.tile(d['quality'][idx], (lvec.shape[0], 1)).T
dd['convergence'] = np.tile(d['convergence'][idx], (lvec.shape[0], 1)).T
dd['status'] = np.tile(d['status'][idx], (lvec.shape[0], 1)).T
dd['latitude'] = np.tile(d['latitude'][idx], (lvec.shape[0], 1)).T
dd['longitude'] = np.tile(d['longitude'][idx], (lvec.shape[0], 1)).T
dd['solar_zenith_angle'] = np.tile(d['solar_zenith_angle'][idx], (lvec.shape[0], 1)).T
dd['solarZenithAngle'] = np.tile(d['solarZenithAngle'][idx], (lvec.shape[0], 1)).T
for k in list(dd.keys()):
dd[k] = np.asarray(dd[k])
dd[k] = dd[k].flatten().tolist()
Expand All @@ -209,9 +206,9 @@ def _do_qc(self, d):
if(d['dateTime'][irec] < self.startTAI or d['dateTime'][irec] > self.endTAI):
continue
for k in list(d.keys()):
if (len(d[k].shape) == 1 and k != 'air_pressure'):
if (len(d[k].shape) == 1 and k != 'pressure'):
dd[k].append(d[k][irec])
elif (k == 'air_pressure'):
elif (k == 'pressure'):
dd[k].append(d[k][ilev])
elif k != 'errKey':
dd[k].append(d[k][irec, ilev])
Expand All @@ -220,7 +217,7 @@ def _do_qc(self, d):

def _read(self):
# set up variable names for IODA
self._setVarAttr('mole_fraction_of_ozone_in_air')
self._setVarAttr(varname_ozone)

# loop through input filenames
for i, f in enumerate(self.filenames):
Expand All @@ -246,8 +243,7 @@ def _read(self):
if(self.errorOn):
self.outdata[self.varDict[iodavar]['errKey']].extend(d['errKey'])

DimDict['nlocs'] = np.float32(len(self.outdata[('dateTime', 'MetaData')]))
AttrData['nlocs'] = np.int32(DimDict['nlocs'])
DimDict['Location'] = np.float32(len(self.outdata[('dateTime', 'MetaData')]))

for k in self.outdata.keys():
self.outdata[k] = np.asarray(self.outdata[k])
Expand Down
69 changes: 39 additions & 30 deletions src/compo/modis_aod2ioda.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -14,6 +14,7 @@
from datetime import datetime, date, timedelta
import os
from pathlib import Path
from pyhdf.SD import SD, SDC
Copy link
Copy Markdown
Collaborator

@srherbener srherbener Jan 3, 2023

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Why switch to pyhdf API? What problem is this solving?

Copy link
Copy Markdown
Contributor Author

@gthompsnJCSDA gthompsnJCSDA Jan 3, 2023

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Because I had to. The input file is hdf4, not hdf5/netcdf, so I had to switch. This is basically coming from a develop branch change that someone made - so I just followed it thru into the sprint branch so it functions again.

Copy link
Copy Markdown
Collaborator

@BenjaminRuston BenjaminRuston Jan 3, 2023

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

also found this strange was this something from the aerosol group? but if pyhdf is supported in the stack would potentially get this through and then have the aerosol group review and hopefully consolidate and minimize the number of libraries used

Copy link
Copy Markdown
Collaborator

@srherbener srherbener Jan 3, 2023

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@gthompsnJCSDA thanks for the explanation. pyhdf is supported/included in spack-stack. It would be nice in the future to get providers updated from hdf4 to netcdf/hdf5, but that is certainly outside the scope of this PR.


IODA_CONV_PATH = Path(__file__).parent/"@SCRIPT_LIB_PATH@"
if not IODA_CONV_PATH.is_dir():
Expand All @@ -33,8 +34,10 @@
obsvars = ["aerosolOpticalDepth"]

# A dictionary of global attributes. More filled in further down.
AttrData = {}
AttrData['ioda_object_type'] = 'AOD at 550nm'
AttrData = {
'converter': os.path.basename(__file__),
'description': 'AOD at 550nm'
}

# A dictionary of variable dimensions.
DimDict = {}
Expand All @@ -52,12 +55,19 @@


class AOD(object):
def __init__(self, filenames, obs_time_arg):
def __init__(self, filenames, obs_time, pltfrm):
self.filenames = filenames
self.obs_time = obs_time_arg
self.obs_time = obs_time
self.pltfrm = pltfrm
self.varDict = defaultdict(lambda: defaultdict(dict))
self.outdata = defaultdict(lambda: DefaultOrderedDict(OrderedDict))
self.varAttrs = DefaultOrderedDict(lambda: DefaultOrderedDict(dict))
# there's absolutely no difference in the hdf4 files attributes
# between Terra and Aqua files. So it is user specified
AttrData['platform'] = pltfrm
# sensor would be always MODIS for this converter
AttrData['sensor'] = 'MODIS'
AttrData['datetimeReference'] = obs_time
self._read()

def _read(self):
Expand All @@ -79,50 +89,45 @@ def _read(self):
self.varAttrs[('dateTime', metaDataName)]['units'] = 'seconds since 1993-01-01T00:00:00Z'

# Make empty lists for the output vars
self.outdata[('latitude', metaDataName)] = np.array([], dtype=np.float64)
self.outdata[('longitude', metaDataName)] = np.array([], dtype=np.float64)
self.outdata[('latitude', metaDataName)] = np.array([], dtype=np.float32)
self.outdata[('longitude', metaDataName)] = np.array([], dtype=np.float32)
self.outdata[('dateTime', metaDataName)] = np.array([], dtype=np.int64)
for iodavar in obsvars:
self.outdata[self.varDict[iodavar]['valKey']] = np.array([], dtype=np.float64)
self.outdata[self.varDict[iodavar]['errKey']] = np.array([], dtype=np.float64)
self.outdata[self.varDict[iodavar]['valKey']] = np.array([], dtype=np.float32)
self.outdata[self.varDict[iodavar]['errKey']] = np.array([], dtype=np.float32)
self.outdata[self.varDict[iodavar]['qcKey']] = np.array([], dtype=np.int32)

# loop through input filenames
for f in self.filenames:
ncd = nc.Dataset(f, 'r')
ncd_str = str(ncd)
if 'Terra' in ncd_str:
AttrData['platform'] = 'Terra'
elif 'Aqua' in ncd_str:
AttrData['platform'] = 'Aqua'
if 'MODIS' in ncd_str:
AttrData['sensor'] = 'v.modis_terra'
hdf = SD(f, SDC.READ)

# Get variables
modis_time = ncd.variables['Scan_Start_Time'][:].ravel()
modis_time = modis_time.astype('float64')
lats = ncd.variables['Latitude'][:].ravel()
modis_time = hdf.select('Scan_Start_Time')[:].ravel()
print(f"length of time var: {len(modis_time)}")
modis_time = modis_time.astype('float32')
lats = hdf.select('Latitude')[:].ravel()
lats = lats.astype('float32')
lons = ncd.variables['Longitude'][:].ravel()
lons = hdf.select('Longitude')[:].ravel()
lons = lons.astype('float32')
aod = ncd.variables['AOD_550_Dark_Target_Deep_Blue_Combined'][:].ravel()
land_sea_flag = ncd.variables['Land_sea_Flag'][:].ravel()
QC_flag = ncd.variables['Land_Ocean_Quality_Flag'][:].ravel()
aod = hdf.select('AOD_550_Dark_Target_Deep_Blue_Combined')[:].ravel()
aod = aod.astype('float64')
land_sea_flag = hdf.select('Land_sea_Flag')[:].ravel()
QC_flag = hdf.select('Land_Ocean_Quality_Flag')[:].ravel()
QC_flag = QC_flag.astype('int32')
sol_zen = ncd.variables['Solar_Zenith'][:].ravel()
sen_zen = ncd.variables['Sensor_Zenith'][:].ravel()
unc_land = ncd.variables['Deep_Blue_Aerosol_Optical_Depth_550_Land_Estimated_Uncertainty'][:].ravel()
sol_zen = hdf.select('Solar_Zenith')[:].ravel()
sen_zen = hdf.select('Sensor_Zenith')[:].ravel()
unc_land = hdf.select('Deep_Blue_Aerosol_Optical_Depth_550_Land_Estimated_Uncertainty')[:].ravel()

# Remove undefined values
pos_index = np.where(aod > 0)
lats = lats[pos_index]
lons = lons[pos_index]
aod = aod[pos_index]
aod = aod[pos_index] * 1E-3 # see scale factor
land_sea_flag = land_sea_flag[pos_index]
QC_flag = QC_flag[pos_index]
sol_zen = sol_zen[pos_index]
sen_zen = sen_zen[pos_index]
unc_land = unc_land[pos_index]
unc_land = unc_land[pos_index] * 1E-3 # see scale factor
modis_time = modis_time[pos_index]
obs_time = np.full(len(modis_time), long_missing_value, dtype=np.int64)
for n, t in enumerate(modis_time):
Expand Down Expand Up @@ -151,7 +156,7 @@ def _read(self):
def main():

# get command line arguments
# Usage: python blah.py -i /path/to/obs/2021060801.nc /path/to/obs/2021060802.nc ... -t Analysis_time /path/to/obs/2021060823.nc
# Usage: python blah.py -i /path/to/obs/2021060801.nc /path/to/obs/2021060802.nc ... -p <Terra or Aqua> -t Analysis_time /path/to/obs/2021060823.nc
# -o /path/to/ioda/20210608.nc
# where the input obs could be for any desired interval to concatenated together. Analysis time is generally the midpoint of
# analysis window.
Expand All @@ -171,6 +176,10 @@ def main():
'-t', '--time',
help="Observation time in global attributes (YYYYMMDDHH)",
type=int, required=True)
required.add_argument(
'-p', '--platform',
help="AQUA or TERRA satellite?",
type=str, required=True)
required.add_argument(
'-o', '--output',
help="path of IODA output file",
Expand All @@ -184,7 +193,7 @@ def main():
obs_time = args.time

# Read in the AOD data
aod_class = AOD(args.input, args.time)
aod_class = AOD(args.input, args.time, args.platform)
# write everything out
writer = iconv.IodaWriter(args.output, locationKeyList, DimDict)
writer.BuildIoda(aod_class.outdata, VarDims, aod_class.varAttrs, AttrData)
Expand Down
Loading