godae_ship2ioda.py

Go to the documentation of this file.

#!/usr/bin/env python3
 
#
# (C) Copyright 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 numpy as np
from datetime import datetime
from scipy.io import FortranFile
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 ship(object):
 
    def __init__(self, filename, date):
 
        self.filenamefilename = filename
        self.datedate = date
 
        # Read ship data
        self._rd_ship_rd_ship()
 
        return
 
    def _rd_ship(self):
        '''
        Read the surface obs data
        Based on subroutine rd_ship in ocn_obs.f
        '''
 
        try:
            fh = FortranFile(self.filenamefilename, mode='r', header_dtype='>u4')
        except IOError:
            raise IOError('%s file not found!' % self.filenamefilename)
        except Exception:
            raise Exception('Unknown error opening %s' % self.filenamefilename)
 
        # data is the dictionary with data structure as in ocn_obs.f
        data = {}
 
        data['n_obs'], data['n_lvl'], data['n_vrsn'] = fh.read_ints('>i4')
 
        print('    number ship obs: %d' % data['n_obs'])
        print('  max number levels: %d' % data['n_lvl'])
        print('file version number: %d' % data['n_vrsn'])
 
        if data['n_obs'] <= 0:
            print('No ship observations to process from %s' % self.filenamefilename)
            return
 
        data['ob_wm'] = fh.read_reals('>i4')
        data['ob_glb'] = fh.read_reals('>f4')
        data['ob_lat'] = fh.read_reals('>f4')
        data['ob_lon'] = fh.read_reals('>f4')
        data['ob_age'] = fh.read_reals('>f4')
        data['ob_clm'] = fh.read_reals('>f4')
        data['ob_qc'] = fh.read_reals('>f4')
        data['ob_rgn'] = fh.read_reals('>f4')
        data['ob_sst'] = fh.read_reals('>f4')
        data['ob_typ'] = fh.read_reals('>i4')
        data['ob_dtg'] = fh.read_record('>S12').astype('U12')
        data['ob_rcpt'] = fh.read_record('>S12').astype('U12')
        data['ob_scr'] = fh.read_record('>S1').astype('U1')
 
        if data['n_vrsn'] <= 2:
            print('verify ob_sign for version = %d' % data['n_vrsn'])
            data['ob_sign'] = fh.read_record('S6').astype('U6')
        else:
            data['ob_sign'] = fh.read_record('>S7').astype('U7')
 
        if data['n_vrsn'] > 1:
            data['ob_csgm'] = fh.read_reals('>f4')
            data['ob_gsgm'] = fh.read_reals('>f4')
            data['ob_rsgm'] = fh.read_reals('>f4')
        else:
            minus999 = np.ones(data['n_obs'], dtype=np.float32) * -999.
            data['ob_csgm'] = minus999
            data['ob_gsgm'] = minus999
            data['ob_rsgm'] = minus999
 
        fh.close()
 
        self.datadata = data
 
        return
 
 
class IODA(object):
 
    def __init__(self, filename, date, varDict, obsList):
        '''
        Initialize IODA writer class,
        transform to IODA data structure and,
        write out to IODA file.
        '''
 
        self.filenamefilename = filename
        self.datedate = date
        self.varDictvarDict = varDict
 
        self.locKeyListlocKeyList = [
            ("latitude", "float"),
            ("longitude", "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)
 
        self.keyDictkeyDict = DefaultOrderedDict(lambda: DefaultOrderedDict(dict))
        for key in self.varDictvarDict.keys():
            value = self.varDictvarDict[key]
            self.keyDictkeyDict[key]['valKey'] = value, self.writerwriter.OvalName()
            self.keyDictkeyDict[key]['errKey'] = value, self.writerwriter.OerrName()
            self.keyDictkeyDict[key]['qcKey'] = value, self.writerwriter.OqcName()
 
        # data is the dictionary containing IODA friendly data structure
        self.datadata = DefaultOrderedDict(lambda: DefaultOrderedDict(dict))
 
        recKey = 0
 
        for obs in obsList:
 
            if obs.data['n_obs'] <= 0:
                print('No ship observations for IODA!')
                continue
 
            for n in range(obs.data['n_obs']):
 
                lat = obs.data['ob_lat'][n]
                lon = obs.data['ob_lon'][n]
                dtg = datetime.strptime(obs.data['ob_dtg'][n], '%Y%m%d%H%M')
 
                locKey = lat, lon, dtg.strftime("%Y-%m-%dT%H:%M:%SZ")
 
                for key in self.varDictvarDict.keys():
 
                    val = obs.data[key][n]
                    err = 0.5
                    qc = (100*obs.data['ob_qc'][n]).astype('i4')
 
                    valKey = self.keyDictkeyDict[key]['valKey']
                    errKey = self.keyDictkeyDict[key]['errKey']
                    qcKey = self.keyDictkeyDict[key]['qcKey']
 
                    self.datadata[recKey][locKey][valKey] = val
                    self.datadata[recKey][locKey][errKey] = err
                    self.datadata[recKey][locKey][qcKey] = qc
 
        (ObsVars, LocMdata, VarMdata) = self.writerwriter.ExtractObsData(self.datadata)
        self.writerwriter.BuildNetcdf(ObsVars, LocMdata, VarMdata, self.AttrDataAttrData)
 
        return
 
 
def main():
 
    desc = 'Convert GODAE binary ship data to IODA netCDF4 format'
    parser = ArgumentParser(description=desc,
                            formatter_class=ArgumentDefaultsHelpFormatter)
    parser.add_argument(
        '-i', '--input', help='name of the binary GODAE ship file',
        type=str, nargs='+', required=True)
    parser.add_argument(
        '-o', '--output', help='name of the output netCDF GODAE ship file',
        type=str, required=True, default=None)
    parser.add_argument(
        '-d', '--date', help='file date',
        type=str, metavar='YYYYMMDDHH', required=True)
 
    args = parser.parse_args()
 
    fList = args.input
    foutput = args.output
    fdate = datetime.strptime(args.date, '%Y%m%d%H')
 
    obsList = []
    for fname in fList:
        obsList.append(ship(fname, fdate))
 
    varDict = {
        'ob_sst': 'sea_surface_temperature',
    }
 
    IODA(foutput, fdate, varDict, obsList)
 
 
if __name__ == '__main__':
    main()