l1bextract_safe_nc.py

Go to the documentation of this file.

#! /usr/bin/env python3
 
# Extractor for OLCI Sentinel 3 and MERIS SAFE formatted L1B files
 
import argparse
from datetime import datetime, timedelta
import netCDF4
import pathlib
import sys
import time
from shutil import copyfile as cp
import xml.etree.ElementTree as ET
 
from seadasutils.netcdf_utils import ncsubset_vars
from seadasutils.pixlin_utils import pixlin
from seadasutils.setupenv import env
 
def parseManifest(manifestfile):
    tree = ET.parse(manifestfile)
    root = tree.getroot()
    radfiles = []
    tiefiles = []
    engfiles = []
    for child in root.find('dataObjectSection'):
        filename = pathlib.PurePath(child.find('byteStream').find('fileLocation').attrib['href']).name
        if 'radiance' in filename:
            radfiles.append(filename)
        elif 'tie_' in filename:
            tiefiles.append(filename)
        elif 'time_coordinates' in filename:
            continue
        else:
            engfiles.append(filename)
 
    return (radfiles, tiefiles, engfiles)
 
def minmax(arr):
    return arr.min(), arr.max()
 
def epoch2000(usec):
    # format Epoch 2000 time (microseconds since 2000-01-01)
    base = datetime(2000, 1, 1, 0, 0, 0)
    t = base + timedelta(microseconds=int(usec))
    return t.strftime('%Y-%m-%dT%H:%M:%S.%fZ')
 
class extract:
 
    def __init__(self, idir, odir=None,
                 north=None, south=None, west=None, east=None,
                 spixl=None, epixl=None, sline=None, eline=None,
                 verbose=False):
        # inputs
        self.idir = pathlib.Path(idir)
        self.odir = pathlib.Path(odir)
        self.north = north
        self.south = south
        self.west = west
        self.east = east
        self.spixl = spixl
        self.epixl = epixl
        self.sline = sline
        self.eline = eline
        self.verbose = verbose
        self.geofile = self.idir / 'geo_coordinates.nc'
        self.timefile = self.idir / 'time_coordinates.nc'
        self.tiefile = self.idir / 'tie_geo_coordinates.nc'
        self.manifest = self.idir / 'xfdumanifest.xml'
 
        # defaults
        self.runtime = None
        self.attrs = None
 
        # unused, but needed by setupenv.py
        self.dirs = {}
        self.ancdir = None
        self.curdir = False
        self.sensor = None
        env(self)  # run setupenv
 
    def runextract(self, files, subset):
        # subset each file
        for filename in files:
            srcfile = self.idir / filename
            if srcfile.exists():
                dstfile = self.odir / filename
                if self.verbose:
                    print('Extracting', srcfile)
 
                ncsubset_vars(srcfile, dstfile, subset, timestamp=self.runtime)
 
                # update global attributes
                with netCDF4.Dataset(dstfile, mode='a') as dst:
                    dst.setncatts(self.attrs)
        return 0
 
    def getpixlin(self):
        if self.verbose:
            print("north={} south={} west={} east={}".
                  format(self.north, self.south, self.west, self.east))
 
        # run lonlat2pixline
        pl = pixlin(geofile=self.geofile,
                    north=self.north, south=self.south,
                    west=self.west, east=self.east,
                    verbose=self.verbose)
        pl.lonlat2pixline(zero=False)  # using 1-based indices
        self.spixl, self.epixl, self.sline, self.eline = \
        (pl.spixl, pl.epixl, pl.sline, pl.eline)
        return pl.status
 
    def run(self):
        # convert to zero-based index
        self.spixl, self.epixl, self.sline, self.eline = \
        (v-1 for v in (self.spixl, self.epixl, self.sline, self.eline))
 
        # check/create output directory
        if not self.odir:
            self.odir = '.'.join([self.idir, 'subset'])
        pathlib.Path(self.odir).mkdir(parents=True, exist_ok=True)
 
        # find tie file endpoints
        with netCDF4.Dataset(self.tiefile, 'r') as src:
            npixl = src.dimensions['tie_columns'].size
            nline = src.dimensions['tie_rows'].size
            dpixl = getattr(src, 'ac_subsampling_factor', 1)  # tie_col_pts
            dline = getattr(src, 'al_subsampling_factor', 1)  # tie_row_pts
        spixl, epixl = [self.spixl, self.epixl + dpixl - 1] // dpixl
        sline, eline = [self.sline, self.eline + dline - 1] // dline
 
        # Make sure tie files have enough points for interpolation.
        # l1_olci.c uses type gsl_interp_cspline, which requires
        # at least 3 pixels in each dimension.
        points_required = 3
 
        pixls_needed = points_required - (epixl - spixl + 1)
        if pixls_needed > 0:
            spixl = max(0, spixl - pixls_needed // 2)
            epixl = min(spixl + points_required, npixl) - 1
            if epixl == npixl - 1:
                spixl = npixl - points_required
 
        lines_needed = points_required - (eline - sline + 1)
        if lines_needed > 0:
            sline = max(0, sline - lines_needed // 2)
            eline = min(sline + points_required, nline) - 1
            if eline == nline - 1:
                sline = nline - points_required
 
        # scale to full-resolution coordinates
        # will be aligned with tie files
        spixl, epixl = [dpixl*v for v in (spixl, epixl)]
        sline, eline = [dline*v for v in (sline, eline)]
        if self.verbose:
            print("spixl={} epixl={} sline={} eline={}".
                  format(spixl+1, epixl+1, sline+1, eline+1))
 
        # find new start, stop times
        with netCDF4.Dataset(self.timefile, 'r') as src:
            ts = src['time_stamp'][[sline, eline]]
            start_time = epoch2000(ts[0])
            stop_time = epoch2000(ts[1])
 
        # find new lat/lon ranges
        with netCDF4.Dataset(self.geofile, 'r') as src:
            lat_min, lat_max = minmax(src['latitude']
                                    [sline:eline, spixl:epixl])
            lon_min, lon_max = minmax(src['longitude']
                                    [sline:eline, spixl:epixl])
 
        # define global attributes
        self.attrs = {'start_time': start_time,
                      'stop_time':  stop_time,
                      'geospatial_lat_min': lat_min,
                      'geospatial_lat_max': lat_max,
                      'geospatial_lon_min': lon_min,
                      'geospatial_lon_max': lon_max }
        self.runtime = time.gmtime()  # same for all files
 
        # extract time_coordinates
        subset = {'rows':   [sline, eline]}
        status = self.runextract([self.timefile.name], subset)
 
        # extract full-resolution files
        subset = {'columns':[spixl, epixl],
                  'rows':   [sline, eline]}
        (radfiles, tiefiles, engfiles) = parseManifest(self.manifest)
        status += self.runextract(radfiles + engfiles, subset)
 
        # extract lower-resolution (tie) files
        subset = {'tie_columns':[spixl, epixl] // dpixl,
                  'tie_rows':   [sline, eline] // dline}
        status += self.runextract(tiefiles, subset)
 
        return status
 
 
if __name__ == "__main__":
 
    # parse command line
    parser = argparse.ArgumentParser(
        description='Extract specified area from OLCI Level 1B files.',
        epilog='Specify either geographic limits or pixel/line ranges, not both.')
    parser.add_argument('-v', '--verbose', help='print status messages',
                        action='store_true')
    parser.add_argument('idir',
                        help='directory containing OLCI Level 1B files')
    parser.add_argument('odir', nargs='?',
                        help='output directory (defaults to "idir.subset")')
 
    group1 = parser.add_argument_group('geographic limits')
    group1.add_argument('-n', '--north', type=float, help='northernmost latitude')
    group1.add_argument('-s', '--south', type=float, help='southernmost latitude')
    group1.add_argument('-w', '--west', type=float, help='westernmost longitude')
    group1.add_argument('-e', '--east', type=float, help='easternmost longitude')
 
    group2 = parser.add_argument_group('pixel/line ranges (1-based)')
    group2.add_argument('--spixl', type=int, help='start pixel')
    group2.add_argument('--epixl', type=int, help='end pixel')
    group2.add_argument('--sline', type=int, help='start line')
    group2.add_argument('--eline', type=int, help='end line')
 
    if len(sys.argv) == 1:
        parser.print_help()
        sys.exit(1)
    args = parser.parse_args()
 
    # initialize
    this = extract(idir=args.idir,
                   odir=args.odir,
                   north=args.north,
                   south=args.south,
                   west=args.west,
                   east=args.east,
                   spixl=args.spixl,
                   epixl=args.epixl,
                   sline=args.sline,
                   eline=args.eline,
                   verbose=args.verbose)
 
    # file checks
    if not this.idir.exists():
        print("ERROR: Directory '{}' does not exist. Exiting.".format(this.idir))
        sys.exit(1)
    if not this.timefile.exists():
        print("ERROR: Timestamp file ({}) not found!".format(this.timefile))
        sys.exit(1)
    if not this.geofile.exists():
        print("ERROR: Geolocation file ({}) not found!".format(this.geofile))
        sys.exit(1)
    if not this.tiefile.exists():
        print("ERROR: Tie file ({}) not found!".format(this.tiefile))
        sys.exit(1)
 
    # input value checks
    goodlatlons = None not in (this.north, this.south, this.west, this.east)
    goodindices = None not in (this.spixl, this.epixl, this.sline, this.eline)
    if (goodlatlons and goodindices):
        print("ERROR: Specify either geographic limits or pixel/line ranges, not both.")
        sys.exit(1)
    elif goodlatlons:
        status = this.getpixlin()
        if status not in (0, 110):
            print("No extract; lonlat2pixline status =", status)
            exit(status)
    elif goodindices:
        pass
    else:
        print("ERROR: Specify all values for either geographic limits or pixel/line ranges.")
        sys.exit(1)
 
    # run
    status = this.run()
 
    # copy the manifest in case we ever need it
    cp(this.manifest, this.odir / 'xfdumanifest.xml')
 
    exit(status)