mk_matchup_sst.py

Go to the documentation of this file.

#!/usr/bin/env python3
# coding: utf-8
 
"""
A Perl script to create and output satellite matchups from a SeaBASS file given an 
OB.DAAC L2 satellite SST file and a valid SeaBASS file containing
lat, lon, date, and time as /field entries.
written by J.Scott on 2018/07/20 (joel.scott@nasa.gov)
"""
 
def main():
 
    import argparse
    import os
    import re
    import subprocess
    from datetime import datetime, timedelta
    from copy import copy
    from math import isnan
    from collections import OrderedDict
    from seabass.SB_support import readSB
 
    parser = argparse.ArgumentParser(formatter_class=argparse.RawTextHelpFormatter,description='''\
      This program create and output satellite matchups from a given SeaBASS file.
 
      REQUIRED inputs:
          1) --sat_file        an OB.DAAC L2 satellite SST file(s)
          2) --seabass_file    a valid SeaBASS file(s) with latitude, longitude, and date-time information as field entries.
 
      Outputs:
          1) the original SeaBASS data
          AND
          2) collocated satellite products as additional columns appended to --seabass_file
 
      Example usage call:
         mk_matchup.py --sat_file=[file name].nc --seabass_file=[file name].sb
 
      Caveats:
        * This script is designed to work with files that have been properly
          formatted according to SeaBASS guidelines (i.e. Files that passed FCHECK).
          Some error checking is performed, but improperly formatted input files
          could cause this script to error or behave unexpectedly. Files
          downloaded from the SeaBASS database should already be properly formatted, 
          however, please email seabass@seabass.gsfc.nasa.gov and/or the contact listed
          in the metadata header if you identify problems with specific files.
 
        * It is always HIGHLY recommended that you check for and read any metadata
          header comments and/or documentation accompanying data files. Information 
          from those sources could impact your analysis.
 
        * Compatibility: This script was developed for Python 3.6.
 
      License:
        /*=====================================================================*/
                         NASA Goddard Space Flight Center (GSFC) 
                 Software distribution policy for Public Domain Software
 
         The fd_matchup.py code is in the public domain, available without fee for 
         educational, research, non-commercial and commercial purposes. Users may 
         distribute this code to third parties provided that this statement appears
         on all copies and that no charge is made for such copies.
 
         NASA GSFC MAKES NO REPRESENTATION ABOUT THE SUITABILITY OF THE SOFTWARE
         FOR ANY PURPOSE. IT IS PROVIDED "AS IS" WITHOUT EXPRESS OR IMPLIED
         WARRANTY. NEITHER NASA GSFC NOR THE U.S. GOVERNMENT SHALL BE LIABLE FOR
         ANY DAMAGE SUFFERED BY THE USER OF THIS SOFTWARE.
        /*=====================================================================*/
      ''',add_help=True)
 
    parser.add_argument('--sat_file', nargs='+', type=argparse.FileType('r'), required=True, help='''\
      REQUIRED: input OB.DAAC Level-2 satellite netCDF file
      ''')
 
    parser.add_argument('--seabass_file', nargs='+', type=argparse.FileType('r'), required=True, help='''\
      REQUIRED: input SeaBASS file
      Must be a valid SeaBASS file, passing FHCHECK with no errors.
      Matched-up satellite variables will be appended as additional fields to the data matrix and relevant headers.
      File must contain latitude and longitude and date-time expressed as FIELD entries.
      ''')
 
    parser.add_argument('--box_size', nargs=1, default=([5]), type=int, help=('''\
      OPTIONAL: box size of the satellite data extract made around the in situ point
      Valid values are odd numbers between 3 and 11, default = 5
      '''))
 
    parser.add_argument('--min_valid_sat_pix', nargs=1, default=([50.0]), type=float, help=('''\
      OPTIONAL: percent minimum valid satellite pixels required to create an extract
      Valid value: (0.0 - 100.0), default = 50.0
      '''))
 
    parser.add_argument('--max_time_diff', nargs=1, default=([0.5]), type=float, help=('''\
      OPTIONAL: maximum time difference between satellite and in situ point
      Valid value: decimal number of hours (0 - 36 hours), default = 3
      '''))
 
    parser.add_argument('--verbose', default=False, action='store_true', help=('''\
      OPTIONAL: Displays reason for failed matchup for each in situ target called.
      '''))
 
    parser.add_argument('--no_header_comment', default=False, action='store_true', help=('''\
      OPTIONAL: Flag to NOT append exclusion criteria to the OFILE header. Useful when running script repeatedly. 
      '''))
 
    args=parser.parse_args()
    dict_args=vars(args)
 
    # input verification
    if ((dict_args["box_size"][0] % 2) == 0) or (dict_args["box_size"][0] > 11) or (dict_args["box_size"][0] < 3):
        parser.error("invalid --box_size specified, must be an ODD integer between 3 and 11")
 
    if (dict_args["min_valid_sat_pix"][0] > 100.0) or (dict_args["min_valid_sat_pix"][0] < 0.0):
        parser.error("invalid --min_valid_sat_pix specified, must be a percentage expressed as a floating point number between 0.0 and 100.0")
 
    if (dict_args["max_time_diff"][0] > 36) or (dict_args["max_time_diff"][0] < 0):
        parser.error("invalid --max_time_diff specified, must be a decimal number between 0 and 36")
    else:
        twin_Hmin = -1 * int(dict_args['max_time_diff'][0])
        twin_Mmin = -60 * (dict_args['max_time_diff'][0] - int(dict_args['max_time_diff'][0]))
        twin_Hmax = 1 * int(dict_args['max_time_diff'][0])
        twin_Mmax = 60 * (dict_args['max_time_diff'][0] - int(dict_args['max_time_diff'][0]))
 
    for filein_sb in dict_args['seabass_file']:
 
        # read and verify SeaBASS file and required fields
        if os.path.isfile(filein_sb.name):
            ds = readSB(filename=filein_sb.name, 
                        mask_missing=False, 
                        mask_above_detection_limit=False, 
                        mask_below_detection_limit=False, 
                        no_warn=True)
        else:
            parser.error('ERROR: invalid --seabass_file specified; does ' + filein_sb.name + ' exist?')
 
        ds.datetime = ds.fd_datetime()
        if not ds.datetime:
            parser.error('missing fields in SeaBASS file -- file must contain a valid FIELDS combination of date/year/month/day/sdy and time/hour/minute/second')
 
        print('Looking for satellite/in situ match-ups for:',filein_sb.name)
 
        for filein_sat in dict_args['sat_file']:
 
            if not re.search('\.nc', filein_sat.name.lower()) or \
               not re.search('l2', filein_sat.name.lower()) or \
               not re.search('sst', filein_sat.name.lower()):
                parser.error("invalid --sat_file specified, must be a Level-2 (L2) OB.DAAC SST netCDF (nc) file")
            else:
                #set l2_flags to check for OC/IOP versus SST/SST4 product suites
                flag_arg = ' ignore_flags=LAND\ NAVFAIL\ NAVWARN' + \
                           ' count_flags=LAND\ NAVFAIL\ NAVWARN'
 
            print('Checking:',filein_sat.name)
            write_flag = 0
 
            # loop through input SeaBASS file data rows
            for row,dt in enumerate(ds.datetime):
 
                # create time range of satellite obs to extract
                tim_min = dt + timedelta(hours=twin_Hmin,minutes=twin_Mmin)
                tim_max = dt + timedelta(hours=twin_Hmax,minutes=twin_Mmax)
 
                # verify lat/lon inputs from file
                try:
                    ds.lon = [float(i) for i in ds.data['lon']]
                    ds.lat = [float(i) for i in ds.data['lat']]
                except Exception as E:
                    print(E)
                    parser.error('Missing fields in SeaBASS file. File must contain lat and lon as fields, or specify --slat and --slon.')
 
                if isnan(ds.lat[row]) or isnan(ds.lon[row]):
                    continue
 
                if abs(ds.lon[row]) > 180.0:
                    parser.error('invalid longitude input: all longitude values in ' + filein_sb.name + ' MUST be between -180/180E deg.')
                if abs(ds.lat[row]) > 90.0:
                    parser.error('invalid latitude input: all latitude values in ' + filein_sb.name + ' MUST be between -90/90N deg.')
 
                # construct sys call to sstval_extract
                sys_call_str = 'sstval_extract' + \
                               ' ifile=' + filein_sat.name + \
                               ' slon=' + str(ds.lon[row]) + \
                               ' slat=' + str(ds.lat[row]) + \
                               ' qual_check=qual_sst' + \
                               ' qual_check_distance=10' + \
                               ' global_att=1' + \
                               ' variable_att=1' + \
                               ' boxsize=' + str(dict_args['box_size'][0]) + flag_arg
                # variable_att flag needed to extract units
                # global_att flag needed to extract sensor/instrument names
 
                pid = subprocess.run(sys_call_str, shell=True, encoding='ascii', stdout=subprocess.PIPE, stderr=subprocess.PIPE)
                if pid.returncode == 99:
                    #if dict_args['verbose']:
                        #print('No matchup: in situ target not in granule.')
                    continue #no valid matchup
                elif pid.returncode == 101 or pid.returncode == 102:
                    parser.error('sstval_extract failed -- only accepts Level-2 (L2) satellite files. ' + \
                                    filein_sat.name + ' is not a valid L2 file')
                elif pid.returncode != 0:
                    print(pid.returncode)
                    print(pid.stdout)
                    print(pid.stderr)
                    parser.error('sstval_extract failed -- verify that the sstval_extract binary is compiled and on your PATH and that ' + \
                                    filein_sat.name + ' exists')
 
                # define structures to keep track of sstval_extract's output files
                file_ls = OrderedDict()
                file_del = []
                var_ls = []
 
                upix_ct = 0
                fpix_ct = dict_args['box_size'][0]^2
                pix_ct  = 0
 
                tims = 0
                tim_sat = 0
 
                cvs = []
 
                # parse the extract information
                file_del.append(filein_sat.name + '.qc');
                try:
                    fileobj = open(filein_sat.name + '.qc','r')
                    lines = fileobj.readlines()
                    for line in lines:
                        newline = re.sub("[\r\n]+",'',line)
                        if 'unflagged_pixel_count' in newline:
                            upix_ct = int(newline.split('=')[1])
                        elif 'flagged_pixel_count' in newline:
                            fpix_ct = int(newline.split('=')[1])
                        elif 'pixel_count' in newline:
                            pix_ct = int(newline.split('=')[1])
                        elif 'time' in newline:
                            try:
                                tims = re.search("(\d+)-(\d+)-(\d+)\s+(\d+):(\d+):(\d+)", newline.split('=')[1]);
                                tim_sat = datetime(year=int(tims.group(1)), \
                                          month=int(tims.group(2)), \
                                          day=int(tims.group(3)), \
                                          hour=int(tims.group(4)), \
                                          minute=int(tims.group(5)), \
                                          second=int(tims.group(6)))
                            except:
                                continue
                        elif 'variables' in newline:
                            var_ls = newline.split('=')[1].split(',')
                            for var in var_ls:
                                file_del.append(filein_sat.name + '.qc.' + var);
                                if 'l2_flags'   in var or \
                                   'cntl_pt_rows'   in var or \
                                   'clat'   in var or \
                                   'clon'   in var or \
                                   'day'   in var or \
                                   'elat'   in var or \
                                   'elon'   in var or \
                                   'msec'   in var or \
                                   'slat'   in var or \
                                   'slon'   in var or \
                                   'year'   in var:
                                    continue
                                file_ls[var] = filein_sat.name + '.qc.' + var
                    fileobj.close()
                except Exception as E:
                    print(E)
                    parser.error(' unable to open and read file ' + filein_sat.name + '.qc')
 
                # parse the satellite nc file information
                file_del.append(filein_sat.name + '.qc.global_attrs');
                [inst, plat] = readValEglobatt(filein_sat.name + '.qc.global_attrs', parser)
                if not inst:
                    inst = 'na'
                if not plat:
                    plat = 'na'
 
                # apply exclusion criteria
                # compute and evaluate the max time diff test
                if tim_sat > tim_max or tim_sat < tim_min:
                    clean_file_lis(file_del)
                    if dict_args['verbose']:
                        print('No matchup: failed MAX_TIME_DIFF, required =',dict_args["max_time_diff"][0],'Exclusion level = 1, Matrix row =',row)
                    continue #no valid matchup
 
                # compute and evaluate the min valid sat pix test
                if (pix_ct - fpix_ct) != 0:
                    if upix_ct >= dict_args['box_size'][0]:
                        pix_thresh = 100.0 * (upix_ct / (pix_ct - fpix_ct))
                        if pix_thresh < dict_args['min_valid_sat_pix'][0]:
                            clean_file_lis(file_del)
                            if dict_args['verbose']:
                                print('No matchup: failed MIN_VALID_SAT_PIX, required =',dict_args['min_valid_sat_pix'][0],'found =',pix_thresh,'Exclusion level = 4, Matrix row =',row)
                            continue #no valid matchup
                    else:
                        clean_file_lis(file_del)
                        if dict_args['verbose']:
                            print('No matchup: failed MIN_VALID_SAT_PIX, extracted satellite pixels less than box size, required =',dict_args['box_size'][0],'found =',upix_ct,'Exclusion level = 3, Matrix row =',row)
                        continue #no valid matchup
                else:
                    clean_file_lis(file_del)
                    if dict_args['verbose']:
                        print('No matchup: failed MIN_VALID_SAT_PIX, division by zero when deriving pix_thresh due to required L2FLAG criteria, Exclusion level = 2, Data row =',row)
                    continue #no valid matchup
 
                write_flag = 1 #only write out (write_flag == true), if matchups found
 
                #save L2_fname
                L2file_varname = inst + '_' + plat + '_l2fname'
                ds.addDataToOutput(row, L2file_varname.lower(), 'none', os.path.basename(filein_sat.name), True)
 
                #save tdiff
                tdiff_varname = inst + '_' + plat + '_tdiff'
                tdiff = tim_sat - dt
                ds.addDataToOutput(row, tdiff_varname.lower(), 'seconds', tdiff.total_seconds(), True)
 
                # save extract-variables
                for var in file_ls:
 
                    #save extracted values for each var in file_lis
                    [fmax, fmin, fmean, fmedian, fstdev, centerval, units] = readValEfile(file_ls[var], parser)
 
                    var_name = inst + '_' + plat + '_' + var.lower() + '_max'
                    ds.addDataToOutput(row, var_name.lower(), units, fmax, True)
 
                    var_name = inst + '_' + plat + '_' + var.lower() + '_min'
                    ds.addDataToOutput(row, var_name.lower(), units, fmin, True)
 
                    var_name = inst + '_' + plat + '_' + var.lower() + '_mean'
                    ds.addDataToOutput(row, var_name.lower(), units, fmean, True)
 
                    var_name = inst + '_' + plat + '_' + var.lower() + '_median'
                    ds.addDataToOutput(row, var_name.lower(), units, fmedian, True)
 
                    var_name = inst + '_' + plat + '_' + var.lower() + '_stddev'
                    ds.addDataToOutput(row, var_name.lower(), units, fstdev, True)
 
                    var_name = inst + '_' + plat + '_' + var.lower() + '_center_pixel_value'
                    ds.addDataToOutput(row, var_name.lower(), units, centerval, True)
 
                clean_file_lis(file_del)
 
            if write_flag == 1:
        
                for line in ds.comments:
                    if 'File ammended by OCSSW match-up maker script' in line:
                        comment_flag = True
                        break
            
                    else:
                        comment_flag = False
        
                if not dict_args['no_header_comment']:
                    ds.comments.append(' ')
                    ds.comments.append(' File ammended by OCSSW match-up maker script: mk_matchup_sst.py on ' + datetime.now().strftime('%Y-%m-%d %H:%M:%S') + '.')
                    ds.comments.append(' WARNINGS: This script does NOT adjust in situ data to water-leaving values')
                    ds.comments.append('           This script does NOT account for potential oversampling by the in situ data in time or space.')
                    ds.comments.append('           If successive calls to this script are made for a single in situ file AND multiple-valid-overpasses exist,')
                    ds.comments.append('               only the data from the last successive call will be saved to the output file. This may NOT be the best')
                    ds.comments.append('               quality satellite data in space and time.')
                    ds.comments.append(' EXCLUSION CRITERIA applied to this satellite file:')
                    ds.comments.append('     Box size of satellite extract = ' + str(dict_args['box_size'][0]) + ' pixels by ' + str(dict_args['box_size'][0]) + ' pixels')
                    ds.comments.append('     Minimum percent valid satellite pixels = ' + str(dict_args['min_valid_sat_pix'][0]))
                    ds.comments.append('     Maximum time difference between satellite and in situ = ' + str(dict_args['max_time_diff'][0]) + ' hours')
                    ds.comments.append('     The qual_sst value varies between 0 (best) and 4 (worst).')
                    ds.comments.append('     The qual_sst_mean (qual_sst_max) is the mean (max) of the ' + \
                                       str(dict_args['box_size'][0]) + ' by ' + str(dict_args['box_size'][0]) + ' pixel satellite extract.')
                    ds.comments.append(' ')
 
                print('Satellite/in situ match-up(s) found.')
                
                ds.writeSBfile(filein_sb.name)
        
            else:
                print('No valid satellite match-ups found.')
            
        print(' ')
 
    return
 
 
def clean_file_lis(file_ls):
    import os
    for d in file_ls:
        try:
            os.remove(d)
        except Exception as E:
            print(E)
            print('WARNING: Cleanup of ',d,' failed. Verify that you have read/write priviledges in the current working directory.')
    return
 
 
def readValEglobatt(fname, parser):
    import re
    inst = ''
    plat = ''
    try:
        fileobj = open(fname,'r')
        lines = fileobj.readlines()
        for line in lines:
            newline = re.sub("[\r\n]+",'',line)
            if 'instrument=' in newline:
                inst = newline.lower().split('=')[1]
            elif 'platform=' in newline:
                plat = newline.lower().split('=')[1]
        fileobj.close()
    except Exception as E:
        print(E)
        parser.error(' unable to open and read file ' + fname)
    return(inst, plat)
 
 
def readValEfile(fname, parser):
 
    import re
 
    missing   = ''
    fmax      = ''
    fmin      = ''
    fmean     = ''
    fmedian   = ''
    fstdev    = ''
    centerval = ''
    units     = ''
 
    try:
 
        fileobj = open(fname,'r')
        lines = fileobj.readlines()
        fileobj.close()
 
        for line in lines:
 
            newline = re.sub("[\r\n]+",'',line)
 
            if '_FillValue' in newline:
                missing = newline.split('=')[1]
 
        for line in lines:
 
            newline = re.sub("[\r\n]+",'',line)
 
            if 'filtered_max' in newline:
                fmax = newline.split('=')[1]
                if fmax == missing:
                    fmax = ''
 
            elif 'filtered_min' in newline:
                fmin = newline.split('=')[1]
                if fmin == missing:
                    fmin = ''
 
            elif 'filtered_mean' in newline:
                fmean = newline.split('=')[1]
                if fmean == missing:
                    fmean = ''
 
            elif 'filtered_median' in newline:
                fmedian = newline.split('=')[1]
                if fmedian == missing:
                    fmedian = ''
 
            elif 'filtered_stddev' in newline:
                fstdev = newline.split('=')[1]
                if fstdev == missing:
                    fstdev = ''
 
            elif 'center_value' in newline:
                centerval = newline.split('=')[1]
                if centerval == missing:
                    centerval = ''
 
            elif 'units' in newline:
                units = re.sub('\s', '_', newline.split('=')[1])
 
    except Exception as E:
 
        print(E)
        parser.error(' unable to open and read file ' + fname)
 
    return(fmax, fmin, fmean, fmedian, fstdev, centerval, units)
 
 
if __name__ == "__main__": main()