argoClim2ioda.py

Go to the documentation of this file.

#!/usr/bin/env python3
 
#
# (C) Copyright 2019-2019 UCAR
#
# 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.
#
 
from __future__ import print_function
import sys
import re
import bisect
from datetime import datetime
import netCDF4 as nc4
import numpy as np
from argparse import ArgumentParser, ArgumentDefaultsHelpFormatter
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 orddicts import DefaultOrderedDict
 
 
class argoClim(object):
 
    def __init__(self, filename, begindate=None, enddate=None):
 
        self.filenamefilename = filename
        self.begindatebegindate = begindate
        self.enddateenddate = enddate
 
        # Read data from file
        self._readData_readData()
 
        return
 
    def _readData(self):
        '''
        Read and store the data.
        This file can be downloaded from:
        http://www.argo.ucsd.edu/Gridded_fields.html
        Global gridded NetCDF Argo only dataset produced by optimal interpolation
        doi:10.1002/2017GL073426
        '''
 
        try:
            nc = nc4.Dataset(self.filenamefilename, 'r')
        except IOError:
            raise IOError('%s file not found!' % self.filenamefilename)
        except Exception:
            raise Exception('Unknown error opening %s' % self.filenamefilename)
 
        # First find the variable in this file
        res = [x for x in list(nc.variables.keys()) if re.search("^ARGO_.*MEAN$", x)]
 
        if res:
            self.varnamevarname = res[0].split('_')[1]
            print("%s contains the variable: %s" % (self.filenamefilename, self.varnamevarname))
        else:
            raise NameError("No useable variable was found in the file %s" % self.filenamefilename)
 
        assert self.varnamevarname in ['TEMPERATURE', 'SALINITY'],\
            "%s is not a valid variable name" % self.varnamevarname
 
        lon = nc.variables['LONGITUDE'][:]
        lat = nc.variables['LATITUDE'][:]
        pres = nc.variables['PRESSURE'][:]
 
        # Get absolute time instead of time since epoch
        dtime = nc.variables['TIME']
        time360day = nc4.num2date(dtime[:], dtime.units, calendar='360_day')
 
        # Convert Datetime360Day calendar to a regular datetime object
        timeArray = np.array([datetime.strptime(x.strftime('%Y%m%d%H'), '%Y%m%d%H') for x in time360day])
 
        bI = 0
        if self.begindatebegindate is not None and timeArray[0] <= self.begindatebegindate <= timeArray[-1]:
            bI = bisect.bisect_left(timeArray, self.begindatebegindate)
 
        eI = -1
        if self.enddateenddate is not None and timeArray[0] <= self.enddateenddate <= timeArray[-1]:
            eI = bisect.bisect_left(timeArray, self.enddateenddate) + 1
 
        if bI == 0 and eI == -1:  # all times
            anomaly = nc.variables['ARGO_%s_ANOMALY' % self.varnamevarname][:]
            time = timeArray[:]
        elif bI == 0 and eI > -1:  # upto a given time
            anomaly = nc.variables['ARGO_%s_ANOMALY' % self.varnamevarname][:eI, :]
            time = timeArray[:eI]
        elif bI > 0 and eI == -1:  # from a given time
            anomaly = nc.variables['ARGO_%s_ANOMALY' % self.varnamevarname][bI:, :]
            time = timeArray[bI:]
        elif bI > 0 and eI > -1:  # from and till given times
            anomaly = nc.variables['ARGO_%s_ANOMALY' % self.varnamevarname][bI:eI, :]
            time = timeArray[bI:eI]
 
        mean = nc.variables['ARGO_%s_MEAN' % self.varnamevarname][:]
 
        # create a full field from mean and anomaly
        fullField = anomaly + np.tile(mean, (anomaly.shape[0], 1, 1, 1))
 
        try:
            nc.close()
        except IOError:
            raise IOError('%s file could not be closed!' % self.filenamefilename)
        except Exception:
            raise Exception('Unknown error closing %s' % self.filenamefilename)
 
        # self.data is the data structure
        self.datadata = {}
        self.datadata['lat'] = lat
        self.datadata['lon'] = lon
        self.datadata['pres'] = pres
        self.datadata['time'] = time
        self.datadata['field'] = fullField.data
 
        return
 
 
class IODA(object):
 
    def __init__(self, filename, date, argo):
        '''
        Initialize IODA writer class,
        transform to IODA data structure and,
        write out to IODA file.
        '''
 
        self.filenamefilename = filename
        self.datedate = date
 
        self.locKeyListlocKeyList = [
            ("latitude", "float"),
            ("longitude", "float"),
            ("pressure", "float"),
            ("datetime", "string")
        ]
 
        self.AttrDataAttrData = {
            'odb_version': 1,
            'date_time_string': self.datedate.strftime("%Y-%m-%dT%H:%M:%SZ")
        }
 
        self.writerwriter = iconv.NcWriter(self.filenamefilename, self.locKeyListlocKeyList)
 
        # data is the dictionary containing IODA friendly data structure
        self.datadata = DefaultOrderedDict(lambda: DefaultOrderedDict(dict))
 
        recKey = 0
        valKey = argo.varname, self.writerwriter.OvalName()
        errKey = argo.varname, self.writerwriter.OerrName()
        qcKey = argo.varname, self.writerwriter.OqcName()
 
        # There has to be a better way than explicitly looping over 4! dimensions
        for t, time in enumerate(argo.data['time']):
            for y, lat in enumerate(argo.data['lat']):
                for x, lon in enumerate(argo.data['lon']):
                    for z, pres in enumerate(argo.data['pres']):
 
                        locKey = lat, lon, pres, time.strftime('%Y-%m-%dT%H:%M:%SZ')
 
                        val = argo.data['field'][t, z, y, x]
                        err = 0.
                        qc = 1
 
                        self.datadata[recKey][locKey][valKey] = val
                        self.datadata[recKey][locKey][errKey] = err
                        self.datadata[recKey][locKey][qcKey] = qc
 
                    recKey += 1
 
        (ObsVars, LocMdata, VarMdata) = self.writerwriter.ExtractObsData(self.datadata)
        self.writerwriter.BuildNetcdf(ObsVars, LocMdata, VarMdata, self.AttrDataAttrData)
 
        return
 
 
def main():
 
    desc = 'Convert ARGO gridded global data to IODA netCDF4 format'
    parser = ArgumentParser(
        description=desc,
        formatter_class=ArgumentDefaultsHelpFormatter)
    parser.add_argument(
        '-i', '--input', help='name of the Global ARGO file',
        type=str, required=True)
    parser.add_argument(
        '-o', '--output', help='name of the output IODA netCDF ARGO file',
        type=str, required=True, default=None)
    parser.add_argument(
        '-d', '--date', help='file date', metavar='YYYYMMDDHH',
        type=str, required=True)
    parser.add_argument(
        '-b', '--begindate', help='end date for data time window', metavar='YYYYMMDDHH',
        type=str, required=False, default=None)
    parser.add_argument(
        '-e', '--enddate', help='end date for data time window', metavar='YYYYMMDDHH',
        type=str, required=False, default=None)
 
    args = parser.parse_args()
 
    filename = args.input
    foutput = args.output
    fdate = datetime.strptime(args.date, '%Y%m%d%H')
    bdate = None if args.begindate is None else datetime.strptime(args.begindate, '%Y%m%d%H')
    edate = None if args.enddate is None else datetime.strptime(args.enddate, '%Y%m%d%H')
 
    argo = argoClim(filename, begindate=bdate, enddate=edate)
 
    IODA(foutput, fdate, argo)
 
 
if __name__ == '__main__':
    main()