prismbil2oci.py

Go to the documentation of this file.

#!/usr/bin/env python3
# -*- coding: utf-8 -*-
'''
written by J.Scott on 2017/10/03 (joel.scott@nasa.gov)
'''
import sys
import argparse
import os
import logging
from collections import OrderedDict
from PRISM_module import rd_prism_hdr, rd_L1B_bil
import numpy as np
import operator
import re
import subprocess
from netCDF4 import Dataset
from datetime import datetime, timedelta
 
__version__ = '1.1'
 
#==========================================================================================================================================
 
def ParseCommandLine(args):
    '''
    Defines and parses command line arguments
    '''
    ocdataroot = os.getenv('OCDATAROOT')
    parser = argparse.ArgumentParser(formatter_class=argparse.RawTextHelpFormatter, 
        description=' Converts PRISM data from the CORAL project to netCDF proxy data in the format of PACE-OCI ', 
        add_help=True)
    parser.add_argument('-v', '--version', action='version', version='%(prog)s ' + __version__)
    parser.add_argument('-i', '--ifile', nargs=1, required=True, type=str, 
        help=' PRISM L1 prm*img.hdr file (required) ')
    parser.add_argument('-o', '--ofile', nargs=1, required=True, type=str, 
        help=' Output netCDF filename (required) ')
    parser.add_argument('--cdlfile', nargs=1, type=str, 
        default=os.path.join(ocdataroot, 'prism', 'OCIP_Level-1C_Data_Structure.cdl'), 
        help=' L1C format spec CDL file name ')
    parser.add_argument('-d', '--debug', action='store_true', default=False)
    parsedArgs=parser.parse_args(args)
    return parsedArgs
 
#==========================================================================================================================================
 
def SetLogger(pargs):
    '''
    Initiates a logger instance
    '''
    lgrName = 'prismbil2oci_%s_T_%s' % (datetime.date(datetime.now()), datetime.time(datetime.now()))
    lgr = logging.getLogger(lgrName)
    fmt = '%(message)s'
    if pargs.debug:
        level = logging.DEBUG
        fmt = '%(asctime)s - %(name)s - %(levelname)s -\
               [%(module)s..%(funcName)s..%(lineno)d] -\
               %(message)s'
        formatter = logging.Formatter(fmt)
        fh = logging.FileHandler('%s.log' % lgrName)
        fh.setLevel(logging.DEBUG)
        fh.setFormatter(formatter)
        lgr.addHandler(fh)
 
    level = logging.INFO
    formatter = logging.Formatter(fmt)
    ch = logging.StreamHandler()
    ch.setLevel(logging.INFO)
    ch.setFormatter(formatter)
    lgr.addHandler(ch)
    lgr.setLevel(level)
    lgr.debug('Logger initialized')
    return lgr
 
#==========================================================================================================================================
 
def read_prism_L1(logger, fname_in):
    '''
    A subroutine to read PRISM L1 data
    '''
    L1_hdr   = OrderedDict()
    L1_data  = OrderedDict()
 
    #create list of hdr fnames
    ls_hdr = []
    ls_hdr.append(fname_in) #img.hdr
 
    if os.path.exists(fname_in.rsplit('_',1)[0] + '_loc_ort.hdr'):
        ls_hdr.append(fname_in.rsplit('_',1)[0] + '_loc_ort.hdr') #loc_ort.hdr
    else:
        ls_hdr.append(fname_in.rsplit('_',1)[0] + '_loc.hdr') #loc.hdr
 
    if os.path.exists(fname_in.rsplit('_',1)[0] + '_obs_ort.hdr'):
        ls_hdr.append(fname_in.rsplit('_',1)[0] + '_obs_ort.hdr') #obs_ort.hdr
    else:
        ls_hdr.append(fname_in.rsplit('_',1)[0] + '_obs.hdr') #obs.hdr
 
    #create list of bil fnames
    ls_bil = [fname.rsplit('.')[0] for fname in ls_hdr]
 
    for fname_hdr,fname_bil in zip(ls_hdr,ls_bil):
        if not os.path.exists(fname_hdr):
            logger.warning('Missing PRISM header file: {:}'.format(fname_hdr))
 
        if not os.path.exists(fname_bil):
            logger.warning('Missing PRISM BIL file: {:}'.format(fname_bil))
 
        logger.info('Reading HDR file ' + fname_hdr)
        [hdr_type, dict_hdr] = rd_prism_hdr(fname_hdr)
 
        #save all hdr content
        L1_hdr[hdr_type] = OrderedDict()
        for key in dict_hdr:
            L1_hdr[hdr_type][key] = dict_hdr[key]
 
        logger.info('Reading BIL file ' + fname_bil)
        dict_data = rd_L1B_bil(fname_bil, dict_hdr)
        L1_data.update(dict_data)
 
    return [L1_hdr, L1_data]
 
#==========================================================================================================================================
 
def subsample_spectral(band_data, band_weights):
    '''
    a subroutine to subsample via weighted averaging
    spectral data about a desired band center
 
    Inputs:    band_data    -- 3-dim numpy array of shape-order (number_of_bands(_to_average), number_of_scans, pixels_per_scan)
               band_weights -- list or 1-dim numpy array of length (number_of_bands(_to_average))
 
    Outputs:   data         -- 2-dim numpy array of shape-order (number_of_scans, pixels_per_scan)
               width        -- floating point value of sub-sampled FWHM value
    '''
    width = np.sum(band_weights)
    
    data_numerator = np.empty_like(band_data)
    for i,weight in zip(range(data_numerator.shape[0]), band_weights):
        data_numerator[i,:,:] = np.multiply(band_data[i,:,:], weight)
    
    data = np.divide(np.sum(data_numerator, axis=0), width)
    
    return [data, width]
 
#==========================================================================================================================================
 
def resample_prism2oci(bands_prism, widths_prism, lt_prism):
    '''
    A subroutine to resample PRISM data to OCI targets to create OCIP data
    '''
 
    #map bands and FWHMs from PRISM to OCI for 350:890nm
    index_ocip  = list(range(0,191))
    bands_ocip  = list(bands_prism[0:191])
    widths_ocip = list(widths_prism[0:191])
 
    #map bands and FWHMs from PRISM to OCI for 940nm
    lis_940bands  = list(range(202,215))
    band_array = np.ndarray((len(lis_940bands),1,1), dtype='f8')
    band_array[:,0,0] = [bands_prism[i] for i in lis_940bands]
    [data,width] = subsample_spectral(band_array, [widths_prism[i] for i in lis_940bands])
    bands_ocip.append(data[0,0])
    widths_ocip.append(width)
 
    #map bands and FWHMs from PRISM to OCI for 1038nm
    lis_1038bands = list(range(239,246))
    band_array = np.ndarray((len(lis_1038bands),1,1), dtype='f8')
    band_array[:,0,0] = [bands_prism[i] for i in lis_1038bands]
    [data,width] = subsample_spectral(band_array, [widths_prism[i] for i in lis_1038bands])
    bands_ocip.append(data[0,0])
    widths_ocip.append(width)
 
    #resample Lt
    lt_ocip = np.ndarray((len(bands_ocip),lt_prism.shape[1],lt_prism.shape[2]), dtype="f8")
    lt_ocip[0:len(index_ocip),:,:] = lt_prism[index_ocip,:,:]
 
    [data,width] = subsample_spectral(lt_prism[lis_940bands,:,:], [widths_prism[i] for i in lis_940bands])
    lt_ocip[-2,:,:] = data[:,:]
 
    [data,width] = subsample_spectral(lt_prism[lis_1038bands,:,:], [widths_prism[i] for i in lis_1038bands])
    lt_ocip[-1,:,:] = data[:,:]
 
    return [bands_ocip, widths_ocip, lt_ocip]
 
#==========================================================================================================================================
 
def read_OCIP_L1C_CDL(fname_in, num_scans, num_pix, num_bands):
    '''
    A subroutine to read a OCIP L1C CDL file and fill the proper dimensions
    '''
    with open(fname_in,'r') as fid:
        lines = fid.readlines()
 
    #define dims and chunksizes in CDL file
    lines = [re.sub('NUMSSCANS',  str(num_scans), line) for line in lines]
    lines = [re.sub('NUMSPIXELS', str(num_pix),   line) for line in lines]
    lines = [re.sub('NUMSBANDS',  str(num_bands), line) for line in lines]
 
    return lines
 
#==========================================================================================================================================
 
def PRISMmain(args):
    '''
    Coordinates calls to PRISM L1 BIL to L1C netCDF conversion process
    '''
    #parse input arguments and instantiate logger
    pArgs = ParseCommandLine(args)
    mainLogger = SetLogger(pArgs)
    mainLogger.info("prismbil2oci.py %s" % __version__)
    ifile = ''.join(pArgs.ifile)
    ofile = ''.join(pArgs.ofile)
    cfile = ''.join(pArgs.cdlfile)
 
    #load L1 data
    [dsL1_hdr, dsL1_data] = read_prism_L1(mainLogger, ifile)
 
    #resample PRISM band centers and widths to more closely match OCI
    [bands_ocip, widths_ocip, lt_ocip] = resample_prism2oci([float(d) for d in dsL1_hdr['rdn_img']['wavelength']], 
                                                            [float(d) for d in dsL1_hdr['rdn_img']['fwhm'][:len(dsL1_hdr['rdn_img']['wavelength'])]], 
                                                            dsL1_data['lt'])
 
    #read and parse CDL file
    cdl_lines = read_OCIP_L1C_CDL(cfile, dsL1_data['utc time'].shape[0], dsL1_data['utc time'].shape[1], len(bands_ocip))
 
    #mk netCDF4 file via NCGEN subprocess call
    if '/' in ofile:
        nc_fname = ofile.rsplit('/',1)[-1]
    else:
        nc_fname = ofile
 
    mainLogger.info('Creating netCDF4 file {:}'.format(ofile))
    pid = subprocess.run('ncgen -b -v3 -o ' + ofile, input=''.join(cdl_lines), encoding='ascii', 
                         shell=True, stdout=subprocess.PIPE, stderr=subprocess.PIPE)
    if pid.stderr:
        logger.warning('Error in ncgen call: {:}'.format(pid.stderr))
 
    nc_fid = Dataset(ofile, 'a')
 
 
    #fill global attributes
    mainLogger.info('Filling global attributes')
    nc_fid.product_name = nc_fname
    nc_fid.prism_version = dsL1_hdr['rdn_img']['prism_version']
    dt_now = datetime.now()
    nc_fid.date_created = dt_now.strftime("%Y-%m-%dT%H:%M:%SZ")
    dt_start = dsL1_hdr['rdn_img']['start datetime']
    nc_fid.flight_line = dt_start.strftime("prm%Y%m%dt%H%M%S")
    nc_fid.time_coverage_start = dt_start.strftime("%Y-%m-%dT%H:%M:%SZ")
    dt_end = dt_start + timedelta(microseconds=np.nanmax(dsL1_data['utc time']))
    nc_fid.time_coverage_end = dt_end.strftime("%Y-%m-%dT%H:%M:%SZ")
    nc_fid.northernmost_latitude = float(np.nanmax(dsL1_data['latitude (wgs-84)']))
    nc_fid.southernmost_latitude = float(np.nanmin(dsL1_data['latitude (wgs-84)']))
    nc_fid.easternmost_longitude = float(np.nanmax(dsL1_data['longitude (wgs-84)']))
    nc_fid.westernmost_longitude = float(np.nanmin(dsL1_data['longitude (wgs-84)']))
    nc_fid.geospatial_lat_max = float(np.nanmax(dsL1_data['latitude (wgs-84)']))
    nc_fid.geospatial_lat_min = float(np.nanmin(dsL1_data['latitude (wgs-84)']))
    nc_fid.geospatial_lon_max = float(np.nanmax(dsL1_data['longitude (wgs-84)']))
    nc_fid.geospatial_lon_min = float(np.nanmin(dsL1_data['longitude (wgs-84)']))
    nc_fid.day_night_flag = "Day"
    nc_fid.earth_sun_distance_correction = float(np.nanmean(dsL1_data['earth-sun distance (au)']))
 
    #fill sbp-group vars
    mainLogger.info('Filling and writing sensor_band_parameters')
    nc_fid.groups['sensor_band_parameters'].variables['wavelength'][:] = bands_ocip
    nc_fid.groups['sensor_band_parameters'].variables['fwhm'][:] = widths_ocip
 
    #fill sla-group vars
    mainLogger.info('Filling and writing scan_line_attributes')
    dsL1_data['scan_start_time'][np.isnan(dsL1_data['scan_start_time'])] = nc_fid.groups['scan_line_attributes'].variables['scan_start_time']._FillValue
    nc_fid.groups['scan_line_attributes'].variables['scan_start_time'][:] = dsL1_data['scan_start_time']
 
    dsL1_data['scan_end_time'][np.isnan(dsL1_data['scan_end_time'])] = nc_fid.groups['scan_line_attributes'].variables['scan_end_time']._FillValue
    nc_fid.groups['scan_line_attributes'].variables['scan_end_time'][:] = dsL1_data['scan_end_time']
 
    #fill nd-group vars
    mainLogger.info('Filling and writing navigation_data')
    dsL1_data['longitude (wgs-84)'][np.isnan(dsL1_data['longitude (wgs-84)'])] = nc_fid.groups['navigation_data'].variables['longitude']._FillValue
    nc_fid.groups['navigation_data'].variables['longitude'][:,:] = dsL1_data['longitude (wgs-84)']
 
    dsL1_data['latitude (wgs-84)'][np.isnan(dsL1_data['latitude (wgs-84)'])] = nc_fid.groups['navigation_data'].variables['latitude']._FillValue
    nc_fid.groups['navigation_data'].variables['latitude'][:,:] = dsL1_data['latitude (wgs-84)']
 
    dsL1_data['elevation (m)'][np.isnan(dsL1_data['elevation (m)'])] = nc_fid.groups['navigation_data'].variables['height']._FillValue
    nc_fid.groups['navigation_data'].variables['height'][:,:] = dsL1_data['elevation (m)']
 
    dsL1_data['path length (m)'][np.isnan(dsL1_data['path length (m)'])] = nc_fid.groups['navigation_data'].variables['range']._FillValue
    nc_fid.groups['navigation_data'].variables['range'][:,:] = dsL1_data['path length (m)']
 
    dsL1_data['to-sensor zenith (0 to 90 degrees from zenith)'][np.isnan(dsL1_data['to-sensor zenith (0 to 90 degrees from zenith)'])] = nc_fid.groups['navigation_data'].variables['sensor_zenith']._FillValue
    nc_fid.groups['navigation_data'].variables['sensor_zenith'][:,:] = dsL1_data['to-sensor zenith (0 to 90 degrees from zenith)']
 
    dsL1_data['to-sensor azimuth (0 to 360 degrees cw from n)'][np.isnan(dsL1_data['to-sensor azimuth (0 to 360 degrees cw from n)'])] = nc_fid.groups['navigation_data'].variables['sensor_azimuth']._FillValue
    nc_fid.groups['navigation_data'].variables['sensor_azimuth'][:,:] = dsL1_data['to-sensor azimuth (0 to 360 degrees cw from n)']
 
    dsL1_data['to-sun zenith (0 to 90 degrees from zenith)'][np.isnan(dsL1_data['to-sun zenith (0 to 90 degrees from zenith)'])] = nc_fid.groups['navigation_data'].variables['solar_zenith']._FillValue
    nc_fid.groups['navigation_data'].variables['solar_zenith'][:,:] = dsL1_data['to-sun zenith (0 to 90 degrees from zenith)']
 
    dsL1_data['to-sun azimuth (0 to 360 degrees cw from n)'][np.isnan(dsL1_data['to-sun azimuth (0 to 360 degrees cw from n)'])] = nc_fid.groups['navigation_data'].variables['solar_azimuth']._FillValue
    nc_fid.groups['navigation_data'].variables['solar_azimuth'][:,:] = dsL1_data['to-sun azimuth (0 to 360 degrees cw from n)']
 
    #fill od-group vars
    mainLogger.info('Filling and writing observation_data')
    lt_ocip[np.isnan(lt_ocip)] = nc_fid.groups['observation_data'].variables['Lt']._FillValue
    nc_fid.groups['observation_data'].variables['Lt'][:,:,:] = lt_ocip
 
    #close nc file
    nc_fid.close()
    mainLogger.info('Successfully created and filled: {:}'.format(ofile))
 
    return
 
if __name__ == "__main__":
    PRISMmain(sys.argv[1:])