Updated running version

bin/ReadAndPlotIASI.py

+#!/usr/bin/env python
+
+"""
+Produce netCDF files of IASI ash detections for met.no
+
+IMPORTANT:
+setenv PYTHONPATH ~/lib/python2.7/site-packages
+in order to use eccodes
+
+"""
+
 import IASITools as IASITools
 from glob import glob
 import os
+import sys
+import argparse
+import atexit
+import numpy as np
+reload(IASITools)
+
+
+### Location of input satellite data
+# at fred:
+#IASI_PATH = "/sat_data/eumetcast/EPS-11/"
+
+# at betty with nsfv4:
+IASI_PATH = "/sat_data_from_fred/eumetcast/EPS-11/"
+
+### Paths to output
+# locally for testing
+#NETCDF_OUT = "/viper/nadir/home/satdata/ash-iasi/netCDF/"
+#PNG_OUT = "/homevip/espen/repos/ash-iasi/Plots/"
+# for met.no
+NETCDF_OUT = "/viper/nadir/nvap/IASI/"
+PNG_OUT = "/viper/nadir/nvap/IASI/"
 
-def ReadAndPlotSingleIASIFile(fn, Projection='lcc', PlotFile=None):
+def at_exit(file_name):
+    """ Delete lock file on exit.
+    Does not quite work as intended, as it doesn't catch signals properly.
+	
+    """
+
+    try:
+        os.unlink(file_name)
+    except Exception as any_err:
+        print "at_exit: error removing lock file: ", any_err
+
+
+def ReadAndPlotSingleIASIFile(fn, Projection='lcc', PlotFile=None, ncfFile=None, verbose=None):
 
         Scene = IASITools.Scene()
-        Scene.ReadBufr(fn, test=False, verbose=False)#True)#
+        ret = Scene.ReadBufr(fn, test=False, verbose=verbose)#True)#
+	if ret != 0:
+		return None
         Scene.IfCrossDateLineCorrectLongitude()
         Scene.CalculateFootprint()
         Scene.ConvertRadiance2BT()
@@ -13,48 +59,168 @@ def ReadAndPlotSingleIASIFile(fn, Projection='lcc', PlotFile=None):
         #    Scene.CalculateSEVIRIChannels()
 
 
+	if (args.write_png):
         # Type is plot type
         # dBT is brightness temperature difference between two IASI channels
-        Scene.PlotLatLonVals(Type='dBT', #'BT', # 'SEVIRIdBT',
-                             Projection=Projection,
+		Scene.PlotLatLonVals(Type='dBT', #'BT', # 'SEVIRIdBT',
+				     Projection=Projection,
                              # llcrnrlon=181,
                              # llcrnrlat=50,
                              # urcrnrlon=120,
                              # urcrnrlat=70,
-                             verbose=False,
-                             Title=False,
-                             test=True, PlotFile=PlotFile)#,
+				     verbose=verbose,
+				     Title=True,
+				     test=True,
+				     PlotFile=PlotFile,
+				     PlotAshContour=False)
+
 
-        return
+	Scene.WriteNetcdf(fnncf)
+
+        return Scene
 
 #############################################################################################
 
-if __name__ == "__main__":
+if __name__ == "__main__":	
+## Program options
+    parser = argparse.ArgumentParser(prog="ReadAndPlotIASI", 
+				     description=":TODO:",
+				     formatter_class=argparse.
+				     ArgumentDefaultsHelpFormatter)
+	
+    parser.add_argument('-v', '--verbose', 
+			default=False,
+			action='store_true',
+			dest='verbose',
+			help="verbose output")
+
+    parser.add_argument('-p', '--plot', 
+			default=False,
+			action='store_true',
+			dest='write_png',
+			help="If true, produce plots")
+
+    parser.add_argument('-n', '--ncf', 
+			default=False,
+			action='store_true',
+			dest='write_ncf',
+			help="If true, produce netcdf files")
+
+    parser.add_argument('-f', '--force', 
+			default=False,
+			action='store_true',
+			dest='force_run',
+			help="If true, force run (reprocessing old files)")
+	
+    args = parser.parse_args()
 
-    path = '/home/aky/NILU/xnilu_wrk/sat_data/events/2011-02_SaharanDust/eumetcast/EPS-11/'
+# Check input arguments
+    if (not args.write_png) and (not args.write_ncf):
+	    print "Nothing to do (give at least one of -p, -n)"
+#	    sys.exit(0)
 
-    fs  = path+"iasi_*_eps_o.l1_bufr"
+# Lock file handling (avoid more than 1 running instance)
+    pid = str(os.getpid())
+    lock_file = os.path.join('/var/lock/', 'ReadAndPlotIASI.lock')
+    if os.path.isfile(lock_file):
+            print "%s already exists, exiting" % lock_file
+#            sys.exit(0)
+    
+    # file(lock_file, 'w').write(pid)
+    # atexit.register(at_exit, lock_file)
+    
+	    
+
+# Start by clearing the netcdf output directory for files older that 14 days
+    IASITools.remove_old_files(NETCDF_OUT, 14, False, dry_run=False)
+
+
+
+
+#    path = "/xnilu_wrk/projects/sat_data/events/2014-02_Kelut/eumetcast/EPS-11/"
+    path = IASI_PATH
+
+    fs  = path+"W_XX-*_eps_o_l1.bin"
     # All file names
     fns =  sorted(glob(fs))
 
+# Just test with 1 file
+#    fns = [fns[0]]
+
+# Invalid files (know to cause core dumps with message like 
+#GEOS_ERROR: TopologyException: found non-noded intersection between LINESTRING (..."
+    invalid_files = ['iasi_20110225_225353_metopa_22592_eps_o.l1_bufr',	
+		     'iasi_20110226_112952_metopa_22600_eps_o.l1_bufr', 
+		     'W_XX-EUMETSAT-Darmstadt,SOUNDING+SATELLITE,METOPA+IASI_C_EUMC_20140214110255_38000_eps_o_l1.bin']
+
+
+    print "Reading BUFR files", fns
+
     ii=0
-    for fn in fns:
-        fnbits = fn.split('/')
-        fnpng  = fnbits[len(fnbits)-1]
-        fnpng  = 'Plots/'+fnpng.replace('.l1_bufr','.png')
+    processed_count = 0
+# Set upper limit of files to process
+    max_process = 150
+    max_process = np.minimum(len(fns), max_process)
+ #   for fn in fns:
+#    while processed_count < max_process:
+    while processed_count < max_process and ii < len(fns):
+	fn = fns[ii]    
+	# if fn.split('/')[-1] in invalid_files:
+	#     continue
+# Plot filename
+	fnbits = fn.split('/')
+	fnpng  = fnbits[len(fnbits)-1]
+#	fnpng  = '../Plots/Kelud/'+fnpng.replace('.l1_bufr','.png')
+#	fnpng  = '../Plots/Kelud/'+fnpng.replace('.bin','.png')
+	fnpng  = PNG_OUT + fnpng.replace(fn.split('.')[-1],'png')
+#	fnpng  = '../Plots/Calbuco/'+fnpng.replace(fn.split('.')[-1],'png')
+
+# ncf filename
+	fnbits = fn.split('/')
+	fnncf  = fnbits[len(fnbits)-1]
+	fnncf  = NETCDF_OUT + fnncf.replace(fn.split('.')[-1],'nc')
+
+
+
+# Skip if netcdf already processed
+	if args.write_ncf and os.path.exists(fnncf) and not args.force_run:
+		print "Netcdf file already exists, skipping"
+		print fnncf
+		ii = ii + 1
+		continue
+
 
-        if os.path.exists(fnpng):
-            continue
+# Skip if input file less limit (some files seems to cantain only 1 scan line
+# and gives problems plotting?)
+	if os.path.getsize(fn) < 2.9*1.e7:
+		print "Skipping file (likely corrupted due to file size):", fn
+		ii = ii + 1
+		continue
 
-        lonmin, lonmax, latmin, latmax = IASITools.GetLatLonRegion(fn)
-        print  lonmin, lonmax, latmin, latmax, fnpng
-        if latmin < -75:
-            Projection='spstere'
-        elif latmax > 75:
-            Projection='npstere'
-        else:
-            Projection='lcc'
+# Add this clause since file may have been deleted since the launch of the program
+	try:
+		print "processing file: ", fn 
+		lonmin, lonmax, latmin, latmax = IASITools.GetLatLonRegion(fn, verbose=args.verbose)
+		if args.verbose:
+			print  lonmin, lonmax, latmin, latmax, fnpng
+		if latmin < -75:
+		    Projection='spstere'
+		elif latmax > 75:
+		    Projection='npstere'
+		else:
+		    Projection='lcc'
 
-        ReadAndPlotSingleIASIFile(fn, Projection=Projection, PlotFile=fnpng)
+		Scene = ReadAndPlotSingleIASIFile(fn, Projection=Projection, PlotFile=fnpng, ncfFile=fnncf, verbose=args.verbose)
+		ii = ii + 1
+		processed_count = processed_count + 1
+		if Scene == None:
+			print "Invalid Scene"
+			continue
+	except:
+                ii = ii + 1
+		continue
 
 #        exit()
+# Remove lock file
+# DEV: this will be removed by atexit?
+#    os.unlink(lock_file)