smap_ssm2ioda.py

Go to the documentation of this file.

#!/usr/bin/env python3
#
# (C) Copyright 2021 EMC/NCEP/NWS/NOAA
#
# This software is licensed under the terms of the Apache Licence Version 2.0
# which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
#
import time, os, sys
import argparse
import netCDF4 as nc
import numpy as np
import re
from datetime import datetime, timedelta
from pathlib import Path
 
IODA_CONV_PATH = Path(__file__).parent/"@SCRIPT_LIB_PATH@"
if not IODA_CONV_PATH.is_dir():
    IODA_CONV_PATH = Path(__file__).parent/'..'/'lib-python'
sys.path.append(str(IODA_CONV_PATH.resolve()))
 
import ioda_conv_ncio as iconv
from collections import defaultdict, OrderedDict
from orddicts import DefaultOrderedDict
 
locationKeyList = [
    ("latitude", "float"),
    ("longitude", "float"),
    ("datetime", "string")
]
 
obsvars = {
    'soil_moisture': 'soilMoistureVolumetric',
}
 
AttrData = {
    'converter': os.path.basename(__file__),
}
 
 
class smap(object):
    def __init__(self, filename, mask, writer):
        self.filenamefilename = filename
        self.maskmask = mask
        self.writerwriter = writer
        self.varDictvarDict = defaultdict(lambda: defaultdict(dict))
        self.outdataoutdata = defaultdict(lambda: DefaultOrderedDict(OrderedDict))
        self.loc_mdataloc_mdata = defaultdict(lambda: DefaultOrderedDict(OrderedDict))
        self.var_mdatavar_mdata = defaultdict(lambda: DefaultOrderedDict(OrderedDict))
        self.unitsunits = {}
        self._read_read()
 
    # Open input file and read relevant info
    def _read(self):
        # set up variable names for IODA
        for iodavar in ['soilMoistureVolumetric']:
            self.varDictvarDict[iodavar]['valKey'] = iodavar, self.writerwriter.OvalName()
            self.varDictvarDict[iodavar]['errKey'] = iodavar, self.writerwriter.OerrName()
            self.varDictvarDict[iodavar]['qcKey'] = iodavar, self.writerwriter.OqcName()
            self.unitsunits[iodavar] = 'm3m-3'
        # open input file name
        ncd = nc.Dataset(self.filenamefilename, 'r')
        # set and get global attributes
        self.satellitesatellite = "SMAP"
        self.sensorsensor = "radar and radiometer"
        AttrData["satellite"] = self.satellitesatellite
        AttrData["sensor"] = self.sensorsensor
 
        data = ncd.groups['Soil_Moisture_Retrieval_Data'].variables['soil_moisture'][:]
        vals = data[:].ravel()
        _FillValue = ncd.groups['Soil_Moisture_Retrieval_Data'].variables['soil_moisture'].getncattr('_FillValue')
        valid_max = ncd.groups['Soil_Moisture_Retrieval_Data'].variables['soil_moisture'].getncattr('valid_max')
        valid_min = ncd.groups['Soil_Moisture_Retrieval_Data'].variables['soil_moisture'].getncattr('valid_min')
 
        lats = ncd.groups['Soil_Moisture_Retrieval_Data'].variables['latitude'][:].ravel()
        lons = ncd.groups['Soil_Moisture_Retrieval_Data'].variables['longitude'][:].ravel()
        errs = ncd.groups['Soil_Moisture_Retrieval_Data'].variables['soil_moisture_error'][:].ravel()
        qflg = ncd.groups['Soil_Moisture_Retrieval_Data'].variables['retrieval_qual_flag'][:].ravel()
        times = np.empty_like(vals, dtype=object)
 
        if self.maskmask == "maskout":
            with np.errstate(invalid='ignore'):
                mask = (vals > valid_min) & (vals < valid_max)
            vals = vals[mask]
            lats = lats[mask]
            lons = lons[mask]
            errs = errs[mask]
            qflg = qflg[mask]
            times = times[mask]
 
        # get datetime from filename
        str_split = self.filenamefilename.split("_")
        str_datetime = str_split[7]
        my_datetime = datetime.strptime(str_datetime, "%Y%m%dT%H%M%S")
        base_datetime = my_datetime.strftime('%Y-%m-%dT%H:%M:%SZ')
 
        qflg = qflg.astype('int32')
        AttrData['date_time_string'] = base_datetime
 
        for i in range(len(lons)):
 
            if vals[i] > 0.0:
                # assumed 4% SM rathern than -999.0
                errs[i] = 0.04*vals[i]
                if qflg[i] > 5:
                    qflg[i] = 0
                else:
                    qflg[i] = 1
            else:
                qflg[i] = 1
 
            times[i] = base_datetime
        self.loc_mdataloc_mdata['datetime'] = self.writerwriter.FillNcVector(times, "datetime")
        self.loc_mdataloc_mdata['latitude'] = lats
        self.loc_mdataloc_mdata['longitude'] = lons
        for iodavar in ['soilMoistureVolumetric']:
            self.outdataoutdata[self.varDictvarDict[iodavar]['valKey']] = vals
            self.outdataoutdata[self.varDictvarDict[iodavar]['errKey']] = errs
            self.outdataoutdata[self.varDictvarDict[iodavar]['qcKey']] = qflg
        self.writerwriter._nvars = len(obsvars)
        self.writerwriter._nlocs = len(self.loc_mdataloc_mdata['datetime'])
 
 
def main():
 
    parser = argparse.ArgumentParser(
        description=('Read SMAP surface soil moisture file(s) and Converter'
                     ' of native h5 format for observations of surface'
                     ' soil moisture to IODA netCDF format.')
    )
    parser.add_argument('-i', '--input',
                        help="name of smap surface soil moisture input file(s)",
                        type=str, required=True)
    parser.add_argument('-o', '--output',
                        help="name of ioda output file",
                        type=str, required=True)
    optional = parser.add_argument_group(title='optional arguments')
    optional.add_argument(
        '-m', '--mask',
        help="maskout missing values: maskout/default, default=none",
        type=str, required=True)
 
    args = parser.parse_args()
 
    writer = iconv.NcWriter(args.output, locationKeyList)
 
    # Read in the profiles
    ssm = smap(args.input, args.mask, writer)
 
    writer.BuildNetcdf(ssm.outdata, ssm.loc_mdata, ssm.var_mdata, AttrData, ssm.units)
 
 
if __name__ == '__main__':
    main()