ocssw V2020
l1mapgen.py
Go to the documentation of this file.
1 #!/usr/bin/env python3
2 import argparse
3 import sys
4 import re
5 import subprocess
6 import os
7 import seadasutils.MetaUtils as MetaUtils
8 from seadasutils.ParamUtils import ParamProcessing
9 
10 def getBinRes(resolution):
11  resvalue = 2000
12  if "km" in resolution:
13  resvalue = float(resolution.strip('km')) * 1000.
14  elif "m" in resolution:
15  resvalue = float(resolution.strip('m'))
16  elif "deg" in resolution:
17  resvalue = float(resolution.strip('deg')) * 111312.
18 
19  if resvalue <= 99.0 :
20  return 'HH'
21  elif resvalue < 250.0:
22  return 'HQ'
23  elif resvalue < 500.0:
24  return 'Q'
25  elif resvalue < 1150.0:
26  return 'H'
27  elif resvalue < 2300.0:
28  return '1'
29  elif resvalue < 4600.0:
30  return '2'
31  elif resvalue < 9200. :
32  return '4'
33  else:
34  return '9'
35 
36 def get_rhots(binfile):
37  rhot = {'MODIS':'rhot_645,rhot_555,rhot_469',
38  'SeaWiFS':'rhot_670,rhot_555,rhot_412',
39  'OCTS':'rhot_670,rhot_565,rhot_412',
40  'VIIRS':'rhot_671,rhot_551,rhot_486',
41  'MERIS':'rhot_665,rhot_560,rhot_412',
42  'OLCI':'rhot_665,rhot_560,rhot_412',
43  'CZCS':'rhot_670,rhot_550,rhot_443',
44  'OLI':'rhot_655,rhot_561,rhot_442',
45  'HAWKEYE':'rhot_670,rhot_556,rhot_447',
46  }
47  cmd = ' '.join(['ncdump', '-h', binfile,'|','grep','instrument'])
48  result = subprocess.Popen(cmd, shell=True,
49  stdout=subprocess.PIPE).communicate()[0]
50  #sensor = ((((result.splitlines()[0]).strip()).rstrip(' ;').split('= "'))[1]).strip('"')
51  sensor = (result.splitlines()[0]).decode().split('"')[1]
52  return rhot[sensor]
53 
54 if __name__ == "__main__":
55 
56  parser = argparse.ArgumentParser(formatter_class=argparse.RawTextHelpFormatter,description='''\
57  This program takes a product_rgb (rhos by default) from a L1 file (or list of files), bins them
58  then maps the L2 binned file using a Plate Carree cylindrical projection,
59  and produces a faux-color png file.
60 
61  The arguments can be specified on the commandline, or put into a parameter file,
62  or the two methods can be used together, with commandline over-riding.''',add_help=True)
63  parser.add_argument('--parfile', type=str, help='input parameter file')
64  parser.add_argument('ifile', nargs='?', type=str, help='input file (L1A/B file or file with list)')
65  parser.add_argument('--geofile', nargs='?', type=str,help='input file geolocation file name')
66  parser.add_argument('--ofile', type=str, help='output file name; default = <ifile>.TC_MAP.<oformat extension>')
67  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")
68  parser.add_argument('--resolution', type=str, default='2.0km', help='''\
69 
70  resolution
71  -----------------
72  #.#: width of a pixel in meters
73  #.#km: width of a pixel in kilometers
74  integer value (e.g. 9km) will result in SMI nominal dimensions
75  (9km == 9.2km pixels, 4km == 4.6km pixels, etc.)
76  #.#deg: width of a pixel in degrees
77  ''')
78  parser.add_argument('--oformat', choices=['netcdf4','hdf4','png','ppm','tiff'], type=str, default="png", help=('''\
79  output file format
80  --------------------------------------------------------
81  netcdf4: netCDF4 file, can contain more than one product
82  hdf4: HDF4 file (old SMI format)
83  png: PNG image file
84  ppm: PPM image file
85  tiff: TIFF file with georeference tags
86  '''))
87  parser.add_argument('--north', type=float, help=('Northern most Latitude (-999=file north)'))
88  parser.add_argument('--south', type=float, help=('Southern most Latitude (-999=file south)'))
89  parser.add_argument('--east', type=float, help=('Eastern most Latitude (-999=file east)'))
90  parser.add_argument('--west', type=float, help=('Western most Latitude (-999=file west)'))
91  parser.add_argument('--projection', type=str,default='platecarree', help='''\
92  enter a proj.4 projection string or choose one of the following
93  predefined projections:
94  --------------------------------------------------------------
95  smi: Standard Mapped image, cylindrical projection, uses
96  central_meridian. n,s,e,w default to whole globe.
97  projection="+proj=eqc +lat_0=<central_meridian>"
98  platecarree: Plate Carree image, cylindrical projection, uses
99  central_meridian
100  projection="+proj=eqc +lat_0=<central_meridian>"
101  mollweide: Mollweide projection
102  projection="+proj=moll +lat_0=<central_meridian>"
103  lambert: Lambert conformal conic projection
104  projection="+proj=lcc +lat_0=<central_meridian>"
105  albersconic: Albers equal-area conic projection
106  projection="+proj=aea +lat_0=<central_meridian>"
107  mercator: Mercator cylindrical map projection
108  projection="+proj=merc +lat_0=<central_meridian>"
109  ease2: Ease Grid 2 projection
110  projection="+proj=cea +lon_0=0 +lat_ts=30 +ellps=WGS84
111  +datum=WGS84 +units=m +lat_0=<central_meridian>"
112  conus: USA Contiguous Albers Equal Area-Conic projection, USGS"
113  projection="+proj=aea +lat_1=29.5 +lat_2=45.5
114  +lat_0=23.0 +lon_0=-96 +x_0=0 +y_0=0
115  +ellps=GRS80 +datum=NAD83 +units=m +no_defs"
116  ''')
117  parser.add_argument('--central_meridian', type=float, help='central meridian for projection in deg east. Only used for smi, mollweide and raw projection')
118  parser.add_argument('--atmocor', action="store_true", help='apply Rayleigh correction')
119  parser.add_argument('--palfile', type=str, help='palette filename. Default uses file for product in product.xml')
120  parser.add_argument('--fudge', type=float,default=1.0, help='fudge factor used to modify size of L3 pixels')
121  parser.add_argument('--threshold', type=float,default=0, help='minimum percentage of filled pixels before an image is generated')
122  parser.add_argument('--keep-intermediates', action='store_true',default=False, help='Do not delete the intermediate L2/L3B files produced')
123  parser.add_argument('-v','--verbose',action='count',default=0, help="Let's get chatty. Two -v's a better than one.")
124 
125  args=parser.parse_args()
126 
127  if args.parfile:
128  param_proc = ParamProcessing(parfile=args.parfile)
129  param_proc.parseParFile()
130  for param in (param_proc.params['main'].keys()):
131  value = param_proc.params['main'][param]
132  if hasattr(args,param):
133  setattr(args,param,value)
134  elif param == 'file':
135  args.ifile = value
136 
137  if not args.ifile:
138  parser.error("You must specify an input file either on the command line or in a parameter file")
139 
140  if not args.ofile:
141  extension = args.oformat
142  if extension == 'netcdf4':
143  extension = 'nc'
144  setattr(args,'ofile',args.ifile + '.TC_MAP.' + extension)
145 
146  default_opts=["product_rgb", "palfile", "north", "south", "east", "west", "central_meridian"]
147  geo_opts = ["north", "south", "east", "west"]
148  l2bin_opts = ["resolution"]
149  l3map_opts = ["resolution", "product_rgb", "oformat", "north", "south", "west", "east", "projection", "central_meridian", "palfile", "fudge", "threshold"]
150  tmpfile_l2gen = "tmp.l2gen"
151  tmpfile_l2bin = "tmp.l2bin"
152 
153  if args.verbose > 1:
154  print(args)
155 
156  ifiles = []
157  geofiles = []
158  if 'text' in MetaUtils.get_mime_data(args.ifile):
159  with open(args.ifile,'r') as lstfile:
160  ifiles = lstfile.readlines()
161  else:
162  ifiles.append(args.ifile)
163 
164  if args.geofile:
165  if 'text' in MetaUtils.get_mime_data(args.geofile):
166  with open(args.geofile,'r') as geolistfile:
167  geofiles = geolistfile.readlines()
168  else:
169  geofiles.append(args.geofile)
170  if len(geofiles) != len(ifiles):
171  parser.error("Number of provided geofiles does not match number of provided ifiles...")
172 
173  l2filelst = open(tmpfile_l2gen,'w')
174  for i, l1file in enumerate(ifiles):
175  print(i,l1file)
176  # Build the l2gen command
177  l1file = l1file.strip()
178  clo = ["l2gen"]
179  clo.append('ifile='+l1file)
180  if args.geofile:
181  clo.append("geofile=" + geofiles[i])
182 
183  ofile = tmpfile_l2gen + '_'+str(i)
184  clo.append("ofile"+"="+ofile)
185 
186  if args.atmocor:
187  clo.append('l2prod=rhos_nnn')
188  else:
189  clo.append('l2prod=rhot_nnn')
190 
191  clo.append("atmocor=0")
192  for geo_opt in geo_opts:
193  if getattr(args,geo_opt):
194  clo.append(geo_opt + "=" + str(getattr(args,geo_opt)))
195 
196  if args.verbose:
197  print("l2gen parameters for %s:" % l1file)
198  print(clo)
199 
200  try:
201  if args.verbose > 1:
202  subprocess.run(clo,check=True)
203  else:
204  subprocess.run(clo,stdout=subprocess.PIPE,stderr=subprocess.STDOUT,check=True)
205  if args.verbose:
206  print("Generated L2 file for input: %s" % l1file)
207  except subprocess.CalledProcessError as e:
208  print("Process error ({0}): message:{1}".format(str(e.returncode), e.output))
209  sys.exit()
210 
211  l2filelst.write(ofile+'\n')
212  l2filelst.close()
213 
214  # Build the l2bin command line
215  clo = ["l2bin"]
216  clo.append('flaguse=ATMFAIL,BOWTIEDEL')
217  clo.append('ifile=' + tmpfile_l2gen)
218  clo.append("ofile=" + tmpfile_l2bin)
219  clo.append("area_weighting=1")
220  for l2bin_opt in l2bin_opts:
221  if getattr(args,l2bin_opt):
222  if l2bin_opt == 'resolution':
223  clo.append(l2bin_opt + "=" + getBinRes(str(getattr(args,l2bin_opt))))
224  else:
225  clo.append(l2bin_opt + "=" + str(getattr(args,l2bin_opt)))
226 
227  if args.verbose:
228  print("l2bin parameters:")
229  print(clo)
230 
231  try:
232  if args.verbose > 1:
233  subprocess.run(clo,check=True)
234  else:
235  subprocess.run(clo,stdout=subprocess.PIPE,stderr=subprocess.STDOUT,check=True)
236  if args.verbose:
237  print("Generated bin file %s" % tmpfile_l2bin)
238  except subprocess.CalledProcessError as e:
239  print("Process error ({0}): message:{1}".format(str(e.returncode), e.output))
240  sys.exit()
241 
242  # Build the l3mapgen command line
243  clo = ["l3mapgen"]
244  clo.append('ifile=' + tmpfile_l2bin)
245  clo.append('ofile=' + args.ofile)
246  clo.append('interp=nearest')
247  clo.append('use_rgb=1')
248  for opt in l3map_opts:
249  if opt == 'product_rgb':
250  if not getattr(args,'atmocor'):
251  clo.append('product_rgb=' + get_rhots(tmpfile_l2bin))
252  elif getattr(args,opt):
253  if type(getattr(args,opt)) is bool: # handle boolean options
254  clo.append(opt + "=1" )
255  else: # handle non-boolean options
256  clo.append(opt + "=" + str(getattr(args,opt)))
257 
258  if args.verbose:
259  print("l3mapgen parameters:")
260  print(clo)
261  try:
262  if args.verbose > 1:
263  subprocess.run(clo,check=True)
264  else:
265  subprocess.run(clo,stdout=subprocess.PIPE,stderr=subprocess.STDOUT,check=True)
266  if args.verbose:
267  print("Generated map file: %s" % args.ofile)
268  except subprocess.CalledProcessError as e:
269  print("Process error ({0}): message:{1}".format(str(e.returncode), e.output))
270  sys.exit()
271 
272  if not args.keep_intermediates:
273  for f in [tmpfile_l2bin,tmpfile_l2gen]:
274  if f == tmpfile_l2gen:
275  l2filelst = open(tmpfile_l2gen,'r')
276  for l2f in l2filelst.readlines():
277  l2f = l2f.strip()
278  if args.verbose:
279  print("removing: " + l2f)
280  os.remove(l2f)
281  if args.verbose:
282  print("removing: " + f)
283  os.remove(f)
def get_rhots(binfile)
Definition: l1mapgen.py:36
def getBinRes(resolution)
Definition: l1mapgen.py:10