OB.DAAC Logo
NASA Logo
Ocean Color Science Software

ocssw V2022
mapgen.py
Go to the documentation of this file.
1 #!/usr/bin/env python3
2 import argparse
3 from datetime import datetime
4 import logging
5 from math import ceil, floor
6 import os
7 import re
8 from shlex import quote
9 import subprocess
10 import sys
11 
12 import seadasutils.MetaUtils as MetaUtils
13 from seadasutils.ParamUtils import ParamProcessing
14 from seadasutils.setupenv import build_executable_path
15 
16 def getBinRes(resolution):
17  resvalue = 2000
18  if "km" in resolution:
19  resvalue = float(resolution.strip('km')) * 1000.
20  elif "m" in resolution:
21  resvalue = float(resolution.strip('m'))
22  elif "deg" in resolution:
23  resvalue = float(resolution.strip('deg')) * 111312.
24 
25  if resvalue <= 99.0 :
26  return 'HH'
27  elif resvalue < 250.0:
28  return 'HQ'
29  elif resvalue < 500.0:
30  return 'Q'
31  elif resvalue < 1150.0:
32  return 'H'
33  elif resvalue < 2300.0:
34  return '1'
35  elif resvalue < 4600.0:
36  return '2'
37  elif resvalue < 9200. :
38  return '4'
39  else:
40  return '9'
41 
42 def get_ftype(filepath):
43  level = 1
44  cmd = ['ncdump', '-h', filepath]
45  result = subprocess.run(cmd, shell=False, capture_output=True, text=True)
46  levelstr = None
47  if not result.returncode:
48  for line in result.stdout.splitlines():
49  if 'processing_level' in line:
50  levelstr = line
51  break
52  if levelstr:
53  if re.search(r'"L2',levelstr):
54  level = 2
55  elif re.search(r'"L3',levelstr):
56  level = 3
57 
58  return level
59 
60 def whoops(errormsg, status=1):
61  logging.error(errormsg)
62  sys.exit(status)
63 
64 def execute_command(command):
65  """
66  Execute a process on the system
67  """
68  logging.info("Running: " + ' '.join(command))
69 
70  result = subprocess.run(command, shell=False, capture_output=True, text=True)
71  std_out = result.stdout
72  err_out = result.stderr
73 
74  if result.returncode:
75  logging.debug(std_out)
76  if err_out:
77  logging.error(err_out)
78  whoops("Exiting: {command} returned a status {status}".format(command=' '.join(command),status=result.returncode), result.returncode)
79  else:
80  if std_out:
81  logging.debug(std_out)
82  if err_out:
83  logging.debug(err_out)
84 
85 def run_clo_program(program, inputs):
86  """
87  Set up and run a program that accepts keyword=value command line arguemnts ala CLO
88  """
89  logging.info('\n###### {program} ######\n'.format(program=program))
90  cmd = [str(build_executable_path(program))]
91 
92  params = ["%s=%s" % (k,v) for k, v in inputs.items()]
93 
94  cmd.extend(params)
95  execute_command(cmd)
96 
97 if __name__ == "__main__":
98 
99  parser = argparse.ArgumentParser(formatter_class=argparse.RawTextHelpFormatter,description='''\
100 generate mapped output from a SeaDAS supported satellite data files
101 arguments can be specified on the commandline or in a parameter file
102 the two methods can be used together, with commandline over-riding the parfile''',add_help=True)
103  parser.add_argument('--parfile', '-p', type=str, help='input parameter file')
104  parser.add_argument('--ifile', '-i', type=str, help='input file or text file list of files')
105  parser.add_argument('--geofile', '-g', type=str,help='geolocation file or text file list of files')
106  parser.add_argument('--ofile', '-o', type=str, help='output file name; default: <ifile>.MAP.<oformat ext>')
107  parser.add_argument('--logfile','-l', type=str, default='mapgen_<timestamp>.log', help='''\
108 log file
109 default: mapgen_<timestamp>.log
110 <timestamp> is in seconds since Jan 1, 1970 00:00:00
111 this file is deleted if verbose is not set and no errors
112 occur during processing''')
113 
114  parser.add_argument('--use_rgb', action='store_true', default=False, help='''\
115 generate an RGB image output
116 default: a pseudo-true color image with bands to use
117  controlled by --product_rgb option''')
118  prodGroup = parser.add_mutually_exclusive_group()
119  prodGroup.add_argument('--product', type=str, help='product(s) to map; comma separated')
120  prodGroup.add_argument('--product_rgb', type=str,help='''\
121 comma separated string of RGB products
122 e.g., product_rgb=rhos_645,rhos_555,rhos_469
123 default: sensor specific, see
124 $OCDATAROOT/<sensor>/l1mapgen_defaults.par''')
125 
126  parser.add_argument('--resolution', '-r', type=str, default='2.0km', help='''\
127  #.#: width of a pixel in meters
128  #.#km: width of a pixel in kilometers
129  #.#deg: width of a pixel in degrees''')
130  parser.add_argument('--oformat', choices=['netcdf4','png','ppm','tiff'], type=str, default="png", help='''\
131 netcdf4: Network Common Data Form v4 file
132  can contain more than one product
133 png: Portable Network Graphics format image
134 ppm: Portable PixMap format image
135 tiff: Tagged Image File Format with georeference tags''')
136  parser.add_argument('--use_transparency','-t',action='store_true', default=False, help='make missing data transparent\nonly valid for color PNG and TIFF output')
137  parser.add_argument('--north', '-n', type=float, help=('northern-most latitude; default: input file max lßatitude'))
138  parser.add_argument('--south', '-s', type=float, help=('southern-most latitude; default: input file min latitude'))
139  parser.add_argument('--east', '-e', type=float, help=('eastern-most latitude; default: input file max longitude'))
140  parser.add_argument('--west', '-w', type=float, help=('western-most latitude; default: input file min longitude'))
141  parser.add_argument('--projection', type=str,default='platecarree', help='''\
142  "proj" projection string or one of the following:
143  platecarree: Plate Carree (cylindrical) projection
144  projection="+proj=eqc +lat_0=<central_meridian>"
145  mollweide: Mollweide projection
146  projection="+proj=moll +lat_0=<central_meridian>"
147  lambert: Lambert conformal conic projection
148  projection="+proj=lcc +lat_0=<central_meridian>"
149  albersconic: Albers equal-area conic projection
150  projection="+proj=aea +lat_0=<central_meridian>"
151  mercator: Mercator cylindrical map projection
152  projection="+proj=merc +lat_0=<central_meridian>"
153  ease2: Ease Grid 2 projection
154  projection="+proj=cea +lon_0=0 +lat_ts=30 +ellps=WGS84
155  +datum=WGS84 +units=m +lat_0=<central_meridian>"''')
156  parser.add_argument('--central_meridian', type=float, help='central meridian to use for projection in degrees east')
157  parser.add_argument('--palfile', type=str, help='palette filename\ndefault: see $OCDATAROOT/common/product.xml')
158  parser.add_argument('--fudge', type=float,default=1.0, help='factor used to modify pixel search radius for mapping')
159  parser.add_argument('--datamin', type=float, help='minimum value for scaling (default from product.xml)')
160  parser.add_argument('--datamax', type=float, help='maximum value for scaling (default from product.xml)')
161  parser.add_argument('--scale_type', choices=['linear','log','arctan'],type=str, help='data scaling method (default from product.xml)')
162  parser.add_argument('--threshold', type=float, help='minimum percentage of filled pixels for image generation\ndefault: 0')
163  parser.add_argument('--trimNSEW', action='store_false',default=True, help='do not trim output to match input NSEW range')
164  parser.add_argument('--write_projtext', action='store_true',default=False, help='write projection information to a text file (for mapgen_overlay script)')
165  parser.add_argument('--keep-intermediates', action='store_true',default=False, help='do not delete the intermediate L2/L3B files produced')
166  parser.add_argument('--verbose', '-v',action='count',default=0, help="let's get chatty; each occurrence increases verbosity\ndefault: error\n-v info -vv debug'")
167 
168  args=parser.parse_args()
169 
170  sanitize = ['ifile','ofile','logfile','geofile','parfile','product','resolution']
171  for arg in vars(args):
172  if arg in sanitize:
173  value = getattr(args,arg)
174  if value:
175  setattr(args,arg,quote(value))
176 
177  if 'timestamp' in args.logfile:
178  setattr(args,'logfile',''.join(['mapgen_',datetime.now().strftime('%s'),'.log']))
179  # Set up the logging levels
180  levels = [logging.ERROR, logging.INFO, logging.DEBUG]
181  # limit verbosity to the number of levels, in case someone thinks if 2 v's are good 20 is better
182  verbosity = max(0, min(args.verbose, len(levels)-1))
183 
184  logging.basicConfig(filename=args.logfile,
185  format='-%(levelname)s- %(asctime)s - %(message)s',
186  level=levels[verbosity])
187  logging.getLogger().addHandler(logging.StreamHandler(sys.stdout))
188 
189  if not os.getenv("OCSSWROOT"):
190  whoops("You must define the OCSSW enviroment for me to help you...")
191 
192  if args.parfile:
193  param_proc = ParamProcessing(parfile=args.parfile)
194  param_proc.parseParFile()
195  for param in (param_proc.params['main'].keys()):
196  value = param_proc.params['main'][param]
197  if param in sanitize:
198  value = quote(value)
199  if hasattr(args,param):
200  setattr(args,param,value)
201  elif param == 'file':
202  args.ifile = value
203 
204  if not args.ifile:
205  parser.print_help()
206  whoops("\nYou must specify an input file either on the command line or in a parameter file")
207 
208  if not os.path.isfile(args.ifile):
209  whoops("\nYou must specify a valid input file either on the command line or in a parameter file")
210 
211  if not args.ofile:
212  extension = args.oformat
213  if extension == 'netcdf4':
214  extension = 'nc'
215  setattr(args,'ofile',args.ifile + '.MAP.' + extension)
216 
217  geo_opts = ["north", "south", "east", "west"]
218  l2bin_opts = ["resolution"]
219  l3map_opts = ["product", "use_rgb","product_rgb", "oformat", "resolution",
220  "projection", "central_meridian", "palfile", "fudge", "threshold",
221  "datamin", "datamax", "scale_type", "palfile", "use_transparency",
222  "trimNSEW","write_projtext"]
223  tmpfile_l2gen = "mapgen_l2gen_list.txt"
224  tmpfile_l2bin = "mapgen_tmp.l2bin.nc"
225 
226  if args.verbose > 1:
227  msg = "Parameters passed to or assumed by mapgen:\n"
228  for key,value in (vars(args)).items():
229  msg += "\t{k}={v}\n".format(k=key,v=value)
230  logging.debug(msg)
231 
232  ftype = None
233  ifiles = []
234  geofiles = []
235  if 'text' in MetaUtils.get_mime_data(args.ifile):
236  with open(args.ifile,'r') as lstfile:
237  ifile_list = lstfile.readlines()
238  for ifile in ifile_list:
239  if not os.path.isfile(ifile.strip()):
240  whoops("Hmmm...seems I don't recognize this input: {ifile}".format(ifile=ifile))
241 
242  ifiles.append(ifile.strip())
243  if ftype:
244  if get_ftype(ifile.strip()) != ftype:
245  whoops("Seems you are mixing input types. While I might be able to make sense of your request, it'd be nicer to pass me just on type of file at a time...")
246  else:
247  ftype = get_ftype(ifile.strip())
248  else:
249  ifiles.append(args.ifile)
250  ftype = get_ftype(args.ifile)
251 
252  if not args.use_rgb and not args.product:
253  whoops("My cousin has the artificial intelligence in the family...I'll need you to tell me what product(s) you wish to map...")
254  if ftype == 3 and len(ifiles) > 1:
255  whoops("I can only map one Level 3 file at a time...")
256 
257  # initialize the command-line option dictionary
258  clo = {}
259  # For L1 inputs we need to process to L2 for the rhos/t
260  if ftype == 1:
261  if not args.use_rgb:
262  whoops("You gave me Level 1 data, but didn't tell me you wanted RGB output, and I don't want to assume... If you want another product, I'll need Level 2 inputs")
263 
264  if args.oformat == 'netcdf4':
265  whoops("netCDF output is not supported for RGB images")
266 
267  # we only care about geolocation files if we need to run l2gen...
268  if args.geofile:
269  if 'text' in MetaUtils.get_mime_data(args.geofile):
270  with open(args.geofile,'r') as geolistfile:
271  geofiles = geolistfile.readlines()
272  else:
273  geofiles.append(args.geofile)
274  if len(geofiles) != len(ifiles):
275  whoops("Number of provided geofiles does not match number of provided ifiles...")
276 
277  l2filelst = open(tmpfile_l2gen,'w')
278  for i, l1file in enumerate(ifiles):
279  if args.verbose > 1:
280  logging.info("{cnt}: {ifile}".format(cnt=i, ifile=l1file))
281 
282  # Build the l2gen command
283  l1file = l1file.strip()
284  clo.clear()
285  clo['ifile']=l1file
286  if args.geofile:
287  clo['geofile'] = geofiles[i].strip()
288 
289  ofile = "{filebase}_{cnt:02}.nc".format(filebase=tmpfile_l2gen.replace('_list.txt',''),cnt=i+1)
290  clo['ofile'] = ofile
291  clo['l2prod'] = 'rhos_nnn'
292  clo['atmocor'] = '0'
293  clo['proc_sst'] = '0'
294 
295  for geo_opt in geo_opts:
296  if getattr(args,geo_opt) is not None:
297  clo[geo_opt] = str(getattr(args,geo_opt))
298 
299  if args.verbose > 1:
300  msg = "l2gen parameters for {ofile}\n".format(ofile=ofile)
301  for key,value in (vars(args)).items():
302  msg += "\t{k}={v}\n".format(k=key,v=value)
303  logging.debug(msg)
304 
305  run_clo_program('l2gen', clo)
306  logging.info("Generated L2 file for input: {ifile}".format(ifile=l1file))
307 
308  l2filelst.write(ofile+'\n')
309  l2filelst.close()
310 
311  # Bin the inputs
312  if ftype != 3:
313  # Build the l2bin command line
314  clo.clear()
315  if args.use_rgb:
316  clo['flaguse'] = 'ATMFAIL,BOWTIEDEL'
317  else:
318  if args.product and re.search(r'sst.?',args.product):
319  clo['suite'] = 'SST'
320  clo['flaguse'] = 'LAND,BOWTIEDEL'
321  else:
322  clo['flaguse'] = 'ATMFAIL,CLDICE,BOWTIEDEL'
323  clo['l3bprod'] = args.product
324  if ftype == 1:
325  clo['ifile'] = tmpfile_l2gen
326  else:
327  clo['ifile'] = args.ifile
328 
329  clo['ofile'] = tmpfile_l2bin
330  clo['area_weighting'] = '1'
331  clo['prodtype'] = "regional"
332 
333  if args.north is not None:
334  clo['latnorth'] = ceil(args.north)
335  if args.south is not None:
336  clo['latsouth'] = floor(args.south)
337  if args.west is not None:
338  clo['lonwest'] = floor(args.west)
339  if args.east is not None:
340  clo['loneast'] = ceil(args.east)
341  for l2bin_opt in l2bin_opts:
342  if getattr(args,l2bin_opt) is not None:
343  if l2bin_opt == 'resolution':
344  clo[l2bin_opt] = getBinRes(str(getattr(args,l2bin_opt)))
345  else:
346  clo[l2bin_opt] = str(getattr(args,l2bin_opt))
347 
348  if args.verbose > 1:
349  msg = "l2bin parameters:\n"
350  for key,value in (vars(args)).items():
351  msg += "\t{k}={v}\n".format(k=key,v=value)
352  logging.debug(msg)
353 
354  run_clo_program('l2bin', clo)
355  logging.info("Generated bin file {binfile}".format(binfile=tmpfile_l2bin))
356 
357  # Build the l3mapgen command line
358  clo.clear()
359  if not args.use_rgb:
360  clo['mask_land'] = 1
361  if ftype != 3:
362  clo['ifile'] = tmpfile_l2bin
363  else:
364  clo['ifile'] = args.ifile
365 
366  clo['interp'] = 'nearest'
367 
368  for opt in l3map_opts + geo_opts:
369  if getattr(args,opt) is not None:
370  # set boolean options if defined
371  if type(getattr(args,opt)) is bool:
372  if getattr(args,opt) is True:
373  clo[opt] = '1'
374  else:
375  clo[opt] = str(getattr(args,opt))
376 
377  if args.product:
378  clo['product'] = args.product
379 
380  if args.oformat == 'netcdf4':
381  clo['ofile'] = args.ofile
382  else:
383  if args.product:
384  if len(args.product.split(',')) > 1:
385  ofile = args.ofile.replace(args.oformat,'') + 'PRODUCT' + '.' + args.oformat
386  clo['ofile'] = ofile
387  else:
388  clo['product'] = args.product
389  clo['ofile'] = args.ofile
390  else:
391  clo['ofile'] = args.ofile
392 
393  # Generate the maps!
394  if args.verbose > 1:
395  msg = "l3mapgen parameters:\n"
396  for key,value in (vars(args)).items():
397  msg += "\t{k}={v}\n".format(k=key,v=value)
398  logging.debug(msg)
399 
400  run_clo_program('l3mapgen', clo)
401  logging.info("Generated map file(s)")
402 
403  # Clean up our crumbs...
404  if not args.keep_intermediates:
405  if os.path.exists(tmpfile_l2gen):
406  l2filelst = open(tmpfile_l2gen,'r')
407  for l2f in l2filelst.readlines():
408  l2f = l2f.strip()
409  if os.path.exists(l2f):
410  if args.verbose:
411  logging.info("removing: {rmfile}".format(rmfile=l2f))
412  os.remove(l2f)
413  l2filelst.close()
414  os.remove(tmpfile_l2gen)
415  if os.path.exists(tmpfile_l2bin):
416  if args.verbose:
417  logging.info("removing: {rmfile}".format(rmfile=tmpfile_l2bin))
418  os.remove(tmpfile_l2bin)
419 
420  # Remove the logfile if empty (because verbosity was zero and no errors were thrown...)
421  if os.stat(args.logfile).st_size == 0:
422  os.remove(args.logfile)
type
Definition: mapgen.py:103
float
Definition: mapgen.py:137
def execute_command(command)
Definition: mapgen.py:64
def build_executable_path(prog_name)
def run_clo_program(program, inputs)
Definition: mapgen.py:85
def getBinRes(resolution)
Definition: mapgen.py:16
format
Definition: mapgen.py:185
def whoops(errormsg, status=1)
Definition: mapgen.py:60
def get_ftype(filepath)
Definition: mapgen.py:42