l1mapgen.py

Go to the documentation of this file.

#!/usr/bin/env python3
import argparse
import sys
import re
import subprocess
import os
import seadasutils.MetaUtils as MetaUtils
from seadasutils.ParamUtils import ParamProcessing
 
def getBinRes(resolution):
    resvalue = 2000
    if "km" in resolution:
        resvalue = float(resolution.strip('km')) * 1000.
    elif "m" in resolution:
        resvalue = float(resolution.strip('m'))
    elif "deg" in resolution:
        resvalue = float(resolution.strip('deg')) * 111312.
 
    if resvalue <= 99.0 :
        return 'HH'
    elif resvalue < 250.0:
        return 'HQ'
    elif resvalue < 500.0:
        return 'Q'
    elif resvalue < 1150.0:
        return 'H'
    elif resvalue < 2300.0:
        return '1'
    elif resvalue < 4600.0:
        return '2'
    elif resvalue < 9200.  :
        return '4'
    else:
        return '9'
 
def get_rhots(binfile):
    rhot = {'MODIS':'rhot_645,rhot_555,rhot_469',
            'SeaWiFS':'rhot_670,rhot_555,rhot_412',
            'OCTS':'rhot_670,rhot_565,rhot_412',
            'VIIRS':'rhot_671,rhot_551,rhot_486',
            'MERIS':'rhot_665,rhot_560,rhot_412',
            'OLCI':'rhot_665,rhot_560,rhot_412',
            'CZCS':'rhot_670,rhot_550,rhot_443',
            'OLI':'rhot_655,rhot_561,rhot_442',
            'HAWKEYE':'rhot_670,rhot_556,rhot_447',
            }
    cmd = ' '.join(['ncdump', '-h', binfile,'|','grep','instrument'])
    result = subprocess.Popen(cmd, shell=True,
                                   stdout=subprocess.PIPE).communicate()[0]
    #sensor = ((((result.splitlines()[0]).strip()).rstrip(' ;').split('= "'))[1]).strip('"')
    sensor = (result.splitlines()[0]).decode().split('"')[1]
    return rhot[sensor]
 
if __name__ == "__main__":
 
    parser = argparse.ArgumentParser(formatter_class=argparse.RawTextHelpFormatter,description='''\
      This program takes a product_rgb (rhos by default) from a L1 file (or list of files), bins them
      then maps the L2 binned file using a Plate Carree cylindrical projection,
      and produces a faux-color png file.
 
      The arguments can be specified on the commandline, or put into a parameter file,
      or the two methods can be used together, with commandline over-riding.''',add_help=True)
    parser.add_argument('--parfile', type=str, help='input parameter file')
    parser.add_argument('ifile', nargs='?', type=str, help='input file (L1A/B file or file with list)')
    parser.add_argument('--geofile', nargs='?', type=str,help='input file geolocation file name')
    parser.add_argument('--ofile', type=str, help='output file name; default = <ifile>.TC_MAP.<oformat extension>')
    parser.add_argument('--product_rgb', type=str,help="3 products (e.g.,product_rgb=rhos_645,rhos_555,rhos_469) to use for RGB.  Default is sensor specific")
    parser.add_argument('--resolution', type=str, default='2.0km', help='''\
 
    resolution
    -----------------
         #.#:  width of a pixel in meters
       #.#km:  width of a pixel in kilometers
               integer value (e.g. 9km) will result in SMI nominal dimensions
               (9km == 9.2km pixels, 4km == 4.6km pixels, etc.)
      #.#deg:  width of a pixel in degrees
    ''')
    parser.add_argument('--oformat', choices=['netcdf4','hdf4','png','ppm','tiff'], type=str, default="png", help=('''\
     output file format
     --------------------------------------------------------
     netcdf4: netCDF4 file, can contain more than one product
     hdf4:    HDF4 file (old SMI format)
     png:     PNG image file
     ppm:     PPM image file
     tiff:    TIFF file with georeference tags
     '''))
    parser.add_argument('--north', type=float, help=('Northern most Latitude (-999=file north)'))
    parser.add_argument('--south', type=float, help=('Southern most Latitude (-999=file south)'))
    parser.add_argument('--east', type=float, help=('Eastern most Latitude (-999=file east)'))
    parser.add_argument('--west', type=float, help=('Western most Latitude (-999=file west)'))
    parser.add_argument('--projection', type=str,default='platecarree',  help='''\
        enter a proj.4 projection string or choose one of the following
        predefined projections:
        --------------------------------------------------------------
        smi:       Standard Mapped image, cylindrical projection, uses
                   central_meridian.  n,s,e,w default to whole globe.
                   projection="+proj=eqc +lat_0=<central_meridian>"
        platecarree: Plate Carree image, cylindrical projection, uses
                   central_meridian
                   projection="+proj=eqc +lat_0=<central_meridian>"
        mollweide: Mollweide projection
                   projection="+proj=moll +lat_0=<central_meridian>"
        lambert:   Lambert conformal conic projection
                   projection="+proj=lcc +lat_0=<central_meridian>"
        albersconic: Albers equal-area conic projection
                   projection="+proj=aea +lat_0=<central_meridian>"
        mercator:  Mercator cylindrical map projection
                   projection="+proj=merc +lat_0=<central_meridian>"
        ease2:     Ease Grid 2 projection
                   projection="+proj=cea +lon_0=0 +lat_ts=30 +ellps=WGS84
                         +datum=WGS84 +units=m +lat_0=<central_meridian>"
        conus:     USA Contiguous Albers Equal Area-Conic projection, USGS"
                   projection="+proj=aea +lat_1=29.5 +lat_2=45.5
                         +lat_0=23.0 +lon_0=-96 +x_0=0 +y_0=0
                         +ellps=GRS80 +datum=NAD83 +units=m +no_defs"
    ''')
    parser.add_argument('--central_meridian', type=float, help='central meridian for projection in deg east.  Only used for smi, mollweide and raw projection')
    parser.add_argument('--atmocor', action="store_true", help='apply Rayleigh correction')
    parser.add_argument('--palfile', type=str, help='palette filename. Default uses file for product in product.xml')
    parser.add_argument('--fudge', type=float,default=1.0, help='fudge factor used to modify size of L3 pixels')
    parser.add_argument('--threshold', type=float,default=0, help='minimum percentage of filled pixels before an image is generated')
    parser.add_argument('--keep-intermediates', action='store_true',default=False, help='Do not delete the intermediate L2/L3B files produced')
    parser.add_argument('-v','--verbose',action='count',default=0, help="Let's get chatty.  Two -v's a better than one.")
 
    args=parser.parse_args()
 
    if args.parfile:
        param_proc = ParamProcessing(parfile=args.parfile)
        param_proc.parseParFile()
        for param in (param_proc.params['main'].keys()):
            value = param_proc.params['main'][param]
            if hasattr(args,param):
                    setattr(args,param,value)
            elif param == 'file':
                args.ifile = value
 
    if not args.ifile:
        parser.error("You must specify an input file either on the command line or in a parameter file")
 
    if not args.ofile:
        extension = args.oformat
        if extension == 'netcdf4':
            extension = 'nc'
        setattr(args,'ofile',args.ifile + '.TC_MAP.' + extension)
 
    default_opts=["product_rgb", "palfile", "north", "south", "east", "west", "central_meridian"]
    geo_opts = ["north", "south", "east", "west"]
    l2bin_opts = ["resolution"]
    l3map_opts = ["resolution", "product_rgb", "oformat", "north", "south", "west", "east", "projection", "central_meridian", "palfile", "fudge", "threshold"]
    tmpfile_l2gen = "tmp.l2gen"
    tmpfile_l2bin = "tmp.l2bin"
 
    if args.verbose > 1:
        print(args)
 
    ifiles = []
    geofiles = []
    if 'text' in MetaUtils.get_mime_data(args.ifile):
        with open(args.ifile,'r') as lstfile:
            ifiles = lstfile.readlines()
    else:
        ifiles.append(args.ifile)
 
    if args.geofile:
        if 'text' in MetaUtils.get_mime_data(args.geofile):
            with open(args.geofile,'r') as geolistfile:
                geofiles = geolistfile.readlines()
        else:
            geofiles.append(args.geofile)
        if len(geofiles) != len(ifiles):
            parser.error("Number of provided geofiles does not match number of provided ifiles...")
 
    l2filelst = open(tmpfile_l2gen,'w')
    for i, l1file in enumerate(ifiles):
        print(i,l1file)
        # Build the l2gen command
        l1file = l1file.strip()
        clo = ["l2gen"]
        clo.append('ifile='+l1file)
        if args.geofile:
            clo.append("geofile=" + geofiles[i])
 
        ofile = tmpfile_l2gen + '_'+str(i)
        clo.append("ofile"+"="+ofile)
 
        if args.atmocor:
            clo.append('l2prod=rhos_nnn')
        else:
            clo.append('l2prod=rhot_nnn')
 
        clo.append("atmocor=0")
        for geo_opt in geo_opts:
            if getattr(args,geo_opt):
                clo.append(geo_opt + "=" + str(getattr(args,geo_opt)))
 
        if args.verbose:
            print("l2gen parameters for %s:" % l1file)
            print(clo)
 
        try:
            if args.verbose > 1:
                subprocess.run(clo,check=True)
            else:
                subprocess.run(clo,stdout=subprocess.PIPE,stderr=subprocess.STDOUT,check=True)
            if args.verbose:
                print("Generated L2 file for input: %s" % l1file)
        except subprocess.CalledProcessError as e:
            print("Process error ({0}): message:{1}".format(str(e.returncode), e.output))
            sys.exit()
 
        l2filelst.write(ofile+'\n')
    l2filelst.close()
 
    # Build the l2bin command line
    clo = ["l2bin"]
    clo.append('flaguse=ATMFAIL,BOWTIEDEL')
    clo.append('ifile=' + tmpfile_l2gen)
    clo.append("ofile=" + tmpfile_l2bin)
    clo.append("area_weighting=1")
    for l2bin_opt in l2bin_opts:
        if getattr(args,l2bin_opt):
            if l2bin_opt == 'resolution':
                clo.append(l2bin_opt + "=" + getBinRes(str(getattr(args,l2bin_opt))))
            else:
                clo.append(l2bin_opt + "=" + str(getattr(args,l2bin_opt)))
 
    if args.verbose:
        print("l2bin parameters:")
        print(clo)
 
    try:
        if args.verbose > 1:
            subprocess.run(clo,check=True)
        else:
            subprocess.run(clo,stdout=subprocess.PIPE,stderr=subprocess.STDOUT,check=True)
        if args.verbose:
            print("Generated bin file %s" % tmpfile_l2bin)
    except subprocess.CalledProcessError as e:
        print("Process error ({0}): message:{1}".format(str(e.returncode), e.output))
        sys.exit()
 
    # Build the l3mapgen command line
    clo = ["l3mapgen"]
    clo.append('ifile=' + tmpfile_l2bin)
    clo.append('ofile=' + args.ofile)
    clo.append('interp=nearest')
    clo.append('use_rgb=1')
    for opt in l3map_opts:
        if opt == 'product_rgb':
            if not getattr(args,'atmocor'):
                clo.append('product_rgb=' + get_rhots(tmpfile_l2bin))
        elif getattr(args,opt):
            if type(getattr(args,opt)) is bool: # handle boolean options
                clo.append(opt + "=1" )
            else: # handle non-boolean options
                clo.append(opt + "=" + str(getattr(args,opt)))
 
    if args.verbose:
        print("l3mapgen parameters:")
        print(clo)
    try:
        if args.verbose > 1:
            subprocess.run(clo,check=True)
        else:
            subprocess.run(clo,stdout=subprocess.PIPE,stderr=subprocess.STDOUT,check=True)
        if args.verbose:
            print("Generated map file: %s" % args.ofile)
    except subprocess.CalledProcessError as e:
        print("Process error ({0}): message:{1}".format(str(e.returncode), e.output))
        sys.exit()
 
    if not args.keep_intermediates:
        for f in [tmpfile_l2bin,tmpfile_l2gen]:
            if f == tmpfile_l2gen:
                l2filelst = open(tmpfile_l2gen,'r')
                for l2f in l2filelst.readlines():
                    l2f = l2f.strip()
                    if args.verbose:
                        print("removing: " + l2f)
                    os.remove(l2f)
            if args.verbose:
                print("removing: " + f)
            os.remove(f)