3 from datetime
import datetime
5 from math
import ceil, floor
8 from shlex
import quote
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.
27 elif resvalue < 250.0:
29 elif resvalue < 500.0:
31 elif resvalue < 1150.0:
33 elif resvalue < 2300.0:
35 elif resvalue < 4600.0:
37 elif resvalue < 9200. :
44 cmd = [
'ncdump',
'-h', filepath]
45 result = subprocess.run(cmd, shell=
False, capture_output=
True, text=
True)
47 if not result.returncode:
48 for line
in result.stdout.splitlines():
49 if 'processing_level' in line:
53 if re.search(
r'"L2',levelstr):
55 elif re.search(
r'"L3',levelstr):
61 logging.error(errormsg)
66 Execute a process on the system
68 logging.info(
"Running: " +
' '.join(command))
70 result = subprocess.run(command, shell=
False, capture_output=
True, text=
True)
71 std_out = result.stdout
72 err_out = result.stderr
75 logging.debug(std_out)
77 logging.error(err_out)
78 whoops(
"Exiting: {command} returned a status {status}".
format(command=
' '.join(command),status=result.returncode), result.returncode)
81 logging.debug(std_out)
83 logging.debug(err_out)
87 Set up and run a program that accepts keyword=value command line arguemnts ala CLO
89 logging.info(
'\n###### {program} ######\n'.
format(program=program))
92 params = [
"%s=%s" % (k,v)
for k, v
in inputs.items()]
97 if __name__ ==
"__main__":
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=
'''\
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''')
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''')
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'")
168 args=parser.parse_args()
170 sanitize = [
'ifile',
'ofile',
'logfile',
'geofile',
'parfile',
'product',
'resolution']
171 for arg
in vars(args):
173 value = getattr(args,arg)
175 setattr(args,arg,quote(value))
177 if 'timestamp' in args.logfile:
178 setattr(args,
'logfile',
''.join([
'mapgen_',datetime.now().strftime(
'%s'),
'.log']))
180 levels = [logging.ERROR, logging.INFO, logging.DEBUG]
182 verbosity = max(0, min(args.verbose, len(levels)-1))
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))
189 if not os.getenv(
"OCSSWROOT"):
190 whoops(
"You must define the OCSSW enviroment for me to help you...")
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:
199 if hasattr(args,param):
200 setattr(args,param,value)
201 elif param ==
'file':
206 whoops(
"\nYou must specify an input file either on the command line or in a parameter file")
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")
212 extension = args.oformat
213 if extension ==
'netcdf4':
215 setattr(args,
'ofile',args.ifile +
'.MAP.' + extension)
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"
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)
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))
242 ifiles.append(ifile.strip())
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...")
249 ifiles.append(args.ifile)
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...")
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")
264 if args.oformat ==
'netcdf4':
265 whoops(
"netCDF output is not supported for RGB images")
269 if 'text' in MetaUtils.get_mime_data(args.geofile):
270 with open(args.geofile,
'r')
as geolistfile:
271 geofiles = geolistfile.readlines()
273 geofiles.append(args.geofile)
274 if len(geofiles) != len(ifiles):
275 whoops(
"Number of provided geofiles does not match number of provided ifiles...")
277 l2filelst = open(tmpfile_l2gen,
'w')
278 for i, l1file
in enumerate(ifiles):
280 logging.info(
"{cnt}: {ifile}".
format(cnt=i, ifile=l1file))
283 l1file = l1file.strip()
287 clo[
'geofile'] = geofiles[i].strip()
289 ofile =
"{filebase}_{cnt:02}.nc".
format(filebase=tmpfile_l2gen.replace(
'_list.txt',
''),cnt=i+1)
291 clo[
'l2prod'] =
'rhos_nnn'
293 clo[
'proc_sst'] =
'0'
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))
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)
306 logging.info(
"Generated L2 file for input: {ifile}".
format(ifile=l1file))
308 l2filelst.write(ofile+
'\n')
316 clo[
'flaguse'] =
'ATMFAIL,BOWTIEDEL'
318 if args.product
and re.search(
r'sst.?',args.product):
320 clo[
'flaguse'] =
'LAND,BOWTIEDEL'
322 clo[
'flaguse'] =
'ATMFAIL,CLDICE,BOWTIEDEL'
323 clo[
'l3bprod'] = args.product
325 clo[
'ifile'] = tmpfile_l2gen
327 clo[
'ifile'] = args.ifile
329 clo[
'ofile'] = tmpfile_l2bin
330 clo[
'area_weighting'] =
'1'
331 clo[
'prodtype'] =
"regional"
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)))
346 clo[l2bin_opt] =
str(getattr(args,l2bin_opt))
349 msg =
"l2bin parameters:\n"
350 for key,value
in (vars(args)).items():
351 msg +=
"\t{k}={v}\n".
format(k=key,v=value)
355 logging.info(
"Generated bin file {binfile}".
format(binfile=tmpfile_l2bin))
362 clo[
'ifile'] = tmpfile_l2bin
364 clo[
'ifile'] = args.ifile
366 clo[
'interp'] =
'nearest'
368 for opt
in l3map_opts + geo_opts:
369 if getattr(args,opt)
is not None:
371 if type(getattr(args,opt))
is bool:
372 if getattr(args,opt)
is True:
375 clo[opt] =
str(getattr(args,opt))
378 clo[
'product'] = args.product
380 if args.oformat ==
'netcdf4':
381 clo[
'ofile'] = args.ofile
384 if len(args.product.split(
',')) > 1:
385 ofile = args.ofile.replace(args.oformat,
'') +
'PRODUCT' +
'.' + args.oformat
388 clo[
'product'] = args.product
389 clo[
'ofile'] = args.ofile
391 clo[
'ofile'] = args.ofile
395 msg =
"l3mapgen parameters:\n"
396 for key,value
in (vars(args)).items():
397 msg +=
"\t{k}={v}\n".
format(k=key,v=value)
401 logging.info(
"Generated map file(s)")
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():
409 if os.path.exists(l2f):
411 logging.info(
"removing: {rmfile}".
format(rmfile=l2f))
414 os.remove(tmpfile_l2gen)
415 if os.path.exists(tmpfile_l2bin):
417 logging.info(
"removing: {rmfile}".
format(rmfile=tmpfile_l2bin))
418 os.remove(tmpfile_l2bin)
421 if os.stat(args.logfile).st_size == 0:
422 os.remove(args.logfile)