Commit a525511f authored by Espen Sollum's avatar Espen Sollum
Browse files

Merge branch 'python3-esolap'

parents 63947c7e 77063ce1
This diff is collapsed.
#!/usr/bin/env python3.7
"""
Plot ash detections for selected IASI (processed) netcdf files, selected
by geographical bounding box and start/end dates
TODO
- add command line options
"""
import IASITools as IASITools
from glob import glob
import os
import sys
import argparse
import atexit
import numpy as np
import re
import shutil
import datetime as dt
from datetime import timedelta
import pdb
import importlib
importlib.reload(IASITools)
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)
# TEST
#--verbose --netcdf_directory=/viper/nadir/nvap/IASI/ --plot_region --lon_min=-30. --lon_max=30. --lat_min=40. --lat_max=80. --projection=PS_NE
def PlotIASIFilteredGeo(nc_directory,
lon_w,
lon_e,
lat_s,
lat_n,
plot_directory,
old_iasi=None,
Projection=None,
start_date=None,
end_date=None,
link_directory=None,
verbose=None,
days_to_process=1):
"""
From a directory of processed IASI netcdf files, make plots
of the ones that have a swath center within the boundaries
specified by input coordinates
INPUT ARGUMENTS
start_date, end_date -- dates on form yyyymmddHHMMSS
TODO:
* add date filter
* take care where longitudes change from negative to positive
"""
lon_w = float(lon_w)
lon_e = float(lon_e)
lat_s = float(lat_s)
lat_n = float(lat_n)
### start_date/end_date converted to datetime objects
if end_date==None:
end_date=dt.datetime.utcnow()
else:
# TODO add check for correct date format
end_date = dt.datetime.strptime(end_date, '%Y%m%d%H%M%S')
if start_date==None:
# search 1 day back in time by default
start_date=dt.datetime.utcnow() - timedelta(days=days_to_process)
else:
start_date = dt.datetime.strptime(start_date, '%Y%m%d%H%M%S')
fs = nc_directory + "*.nc"
fns = sorted(glob(fs))
# Loop over files
for f in fns:
s = IASITools.Scene()
s.ReadNetCDF(f)
s.CenterGeo()
# Extract date info to check for date limits, and use for plot title
fnbits = f.split('/')
fnpng = fnbits[len(fnbits)-1]
fnpng = plot_directory + fnpng.replace(f.split('.')[-1],'png')
if old_iasi:
ydm = re.findall("[0-9]{8}",fnpng)[0]
hms = re.findall("[0-9]_[0-9]{6}",fnpng)[0]
date_str = ydm + hms[2::]
else:
date_str = re.findall("[0-9]{14}",fnpng)[0]
file_date = dt.datetime.strptime(date_str, '%Y%m%d%H%M%S')
####################################################################
### Check if swath within specified time limits
####################################################################
if (file_date > end_date) or (file_date < start_date):
del(s)
continue
####################################################################
### Check if swath within geographic bounds
####################################################################
# TODO correct longitudes if scene crosses greenwich
# l_min = np.min(s.lons[:])
# l_max = np.max(s.lons[:])
# if l_max > 270.0 and l_min < 90:
# s.lons[np.where(s.lons > 180.0)] = s.lons[np.where(s.lons > 180.0)] - 360.0
# print "s.center_lon, s.center_lat, lon_w, lon_e, lat_s, lat_n\n",\
# s.center_lon, s.center_lat, lon_w, lon_e, lat_s, lat_n
#s.center_lon, s.center_lat, lon_w, lon_e, lat_s, lat_n
#26.8290200103 59.9700449311 -30. 30. 40. 80.
#False True False True
# print s.center_lon>lon_w, s.center_lon<lon_e, s.center_lat>lat_s, s.center_lat<lat_n
if (s.center_lon>lon_w) and (s.center_lon<lon_e) and (s.center_lat>lat_s) and (s.center_lat<lat_n):
# Plot (TODO: add date etc to title)
# print("PLOTTING")
man_title = date_str[0:4] +'-'+date_str[4:6]+'-'+date_str[6:8]+' '+date_str[8:10]+':'+date_str[10:12] + ' UTC'
# pdb.set_trace()
if verbose:
print("Plotting file ", fnpng)
s.PlotLatLonVals(Title=False,
ManualTitle=man_title,
Type='dBT_already_calculated',
test=True,
PlotAshContour=False,
PlotFile=fnpng,
Projection=Projection)
# Copy/link the file to a separate directory
if link_directory is not None:
if os.path.isdir(link_directory):
shutil.copy(f, link_directory)
else:
print("File operation failed")
del(s)
else:
del(s)
continue
if __name__ == "__main__":
parser = argparse.ArgumentParser(
prog="PlotSelectedNetCDF",
description="Plot IASI ash detections from processed netcdf files")
parser.add_argument(
'-V', '--verbose', action='store_true', default=False,
dest='verbose',
help="Print extra program output to stdout")
parser.add_argument(
'--plot_region', action='store_true', default=False,
dest='plot_region',
help="Plot all files within specified region")
parser.add_argument(
'--plot_multi', action='store_true', default=False,
dest='plot_multi',
help="Concatenate a list of netcdf files and combine in one plot")
parser.add_argument(
'--netcdf_directory', default="/viper/nadir/nvap/IASI/",
dest='netcdf_directory',
help="")
parser.add_argument(
'--plot_directory', default="/viper/nadir/nvap/IASI/plots/",
dest='plot_directory',
help="")
parser.add_argument(
'--lon_min', default=None,
dest='lon_min',
help="")
parser.add_argument(
'--lon_max', default=None,
dest='lon_max',
help="")
parser.add_argument(
'--lat_min', default=None,
dest='lat_min',
help="")
parser.add_argument(
'--lat_max', default=None,
dest='lat_max',
help="")
parser.add_argument(
'--start_date', default=None,
dest='start_date',
help="String on format YYYYmmddHHMMSS")
parser.add_argument(
'--end_date', default=None,
dest='end_date',
help="String on format YYYYmmddHHMMSS")
parser.add_argument(
'--days_to_process', default=1,
dest='days_to_process',
help="For continuous processing, how many days back in time to look for data")
parser.add_argument(
'--link_directory', default=None,
dest='link_directory',
help="Copy (or link) the selected files to this directory.")
parser.add_argument(
'--projection', default='PS_NE', #'lcc',
dest='projection',
help="Map projection for plotting. Using 'PS_NE' will set "
"a polar stereographic projection for northern europe "
"as desired by met.no" )
parser.add_argument(
'--old_iasi', action='store_true', default=False,
dest='old_iasi',
help="Reading old IASI format")
# Lock file handling (avoid more than 1 running instance)
pid = str(os.getpid())
lock_file = os.path.join('/var/lock/', 'PlotSelectedNCF.lock')
if os.path.isfile(lock_file):
print("%s already exists, exiting" % lock_file)
sys.exit(0)
open(lock_file, 'w').write(pid)
atexit.register(at_exit, lock_file)
args = parser.parse_args()
start_date = args.start_date
end_date = args.end_date
lon_min = args.lon_min
lon_max = args.lon_max
lat_min = args.lat_min
lat_max = args.lat_max
# nc_directory = config.NETCDF_OUT
if args.netcdf_directory == None:
nc_directory = "/homevip/espen/tmp/IASI/" # for testing
else:
nc_directory = args.netcdf_directory
# plot_directory = "/viper/nadir/nvap/IASI/plots/"
plot_directory = args.plot_directory
if args.plot_region:
### Plot all files in netcdf directory that have swath center within coordinate
### bounds
if lon_min is None:
lon_min = -30.
if lon_max is None:
lon_max = 30.
if lat_min is None:
lat_min = 40.
if lat_max is None:
lat_max = 80.
PlotIASIFilteredGeo(nc_directory, lon_min, lon_max, lat_min, lat_max,
plot_directory, start_date=start_date,
end_date=end_date, Projection=args.projection,
old_iasi=args.old_iasi,link_directory=args.link_directory,
verbose=args.verbose,
days_to_process=int(args.days_to_process))
if args.plot_multi:
fn1 = nc_directory + "W_XX-EUMETSAT-Darmstadt,SOUNDING+SATELLITE,METOPA+IASI_C_EUMC_20180110224455_58266_eps_o_l1.nc"
# fn2 = nc_directory + "W_XX-EUMETSAT-Darmstadt,SOUNDING+SATELLITE,METOPA+IASI_C_EUMC_20180110224159_58266_eps_o_l1.nc"
# fn1 = nc_directory + "W_XX-EUMETSAT-Darmstadt,SOUNDING+SATELLITE,METOPA+IASI_C_EUMC_20180110085655_58258_eps_o_l1.nc"
fn2 = nc_directory + "W_XX-EUMETSAT-Darmstadt,SOUNDING+SATELLITE,METOPA+IASI_C_EUMC_20180110085359_58258_eps_o_l1.nc"
fn = [fn1, fn2]
s = IASITools.Scene()
s.ReadNetCDF(fn)
s.PlotLatLonVals(Title=False,
Type='dBT_already_calculated',
test=True,
# PlotFile=plot_directory + '/0853-0856.png',
# ManualTitle='08:53 UTC + 08:56 UTC',
# PlotFile=plot_directory + '/2241-2244.png',
# ManualTitle='22:41 UTC + 22:44 UTC',
PlotFile=plot_directory + '/0853-2244.png',
ManualTitle='08:53 UTC + 22:44 UTC',
NScans=1)
#!/usr/bin/env python
#!/usr/bin/env /usr/local/bin/python3.7
"""
Produce netCDF files of IASI ash detections for met.no
......@@ -16,7 +16,9 @@ import sys
import argparse
import atexit
import numpy as np
reload(IASITools)
import importlib
import datetime as dt
importlib.reload(IASITools)
### Location of input satellite data
......@@ -24,203 +26,235 @@ reload(IASITools)
#IASI_PATH = "/sat_data/eumetcast/EPS-11/"
# at betty with nsfv4:
IASI_PATH = "/sat_data_from_fred/eumetcast/EPS-11/"
#IASI_PATH = "/sat_data_from_fred/eumetcast/EPS-11/"
# local files / 2010 Eyja
IASI_PATH = "/home/eso/repos/ash-iasi/data/0416/"
### Paths to output
# locally for testing
### locally for testing @betty
#NETCDF_OUT = "/viper/nadir/home/satdata/ash-iasi/netCDF/"
#PNG_OUT = "/homevip/espen/repos/ash-iasi/Plots/"
# for met.no
# PNG_OUT = "/viper/nadir/home/satdata/ash-iasi/netCDF/"
### locally for testing at esolap
# NETCDF_OUT = "/home/eso/repos/ash-iasi/ncf/"
# PNG_OUT = "/home/eso/repos/ash-iasi/Plots/"
#### for met.no
NETCDF_OUT = "/viper/nadir/nvap/IASI/"
PNG_OUT = "/viper/nadir/nvap/IASI/"
args = 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
print("at_exit: error removing lock file: ", any_err)
def ReadAndPlotSingleIASIFile(fn, Projection='lcc', PlotFile=None, ncfFile=None, verbose=None):
def ReadAndPlotSingleIASIFile(
fn, Projection='lcc', PlotFile=None, ncfFile=None, verbose=None, old_iasi=False):
Scene = IASITools.Scene()
ret = Scene.ReadBufr(fn, test=False, verbose=verbose)#True)#
if ret != 0:
return None
ret = Scene.ReadBufr(fn, test=False, verbose=verbose, old_iasi=old_iasi)#True)#
if ret != 0:
return None
Scene.IfCrossDateLineCorrectLongitude()
Scene.CalculateFootprint()
Scene.ConvertRadiance2BT()
# Scene.CalculateSEVIRIChannels()
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,
# llcrnrlon=181,
# llcrnrlat=50,
# urcrnrlon=120,
# urcrnrlat=70,
verbose=verbose,
Title=True,
test=True,
PlotFile=PlotFile,
PlotAshContour=False)
# if getattr(args, 'write_png', False):
if PlotFile is not None:
# Type is plot type
# dBT is brightness temperature difference between two IASI channels
Scene.PlotLatLonVals(Type='dBT', #'BT', # 'SEVIRIdBT',
Projection=Projection,
# llcrnrlon=181,
# llcrnrlat=50,
# urcrnrlon=120,
# urcrnrlat=70,
verbose=verbose,
Title=True,
test=True,
PlotFile=PlotFile,
PlotAshContour=False)
Scene.WriteNetcdf(fnncf)
Scene.WriteNetcdf(ncfFile, verbose=verbose)
return Scene
#############################################################################################
if __name__ == "__main__":
if __name__ == "__main__":
## Program options
parser = argparse.ArgumentParser(prog="ReadAndPlotIASI",
description=":TODO:",
formatter_class=argparse.
ArgumentDefaultsHelpFormatter)
description=":DESCROPTION:",
formatter_class=argparse.
ArgumentDefaultsHelpFormatter)
parser.add_argument('-v', '--verbose',
default=False,
action='store_true',
dest='verbose',
help="verbose output")
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")
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")
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)")
default=False,
action='store_true',
dest='force_run',
help="If true, force run (reprocessing old files)")
parser.add_argument('--old',
default=False,
action='store_true',
dest='old_iasi',
help="If true, reads old IASI format.")
parser.add_argument('--max_files',
default=150,
dest='max_files',
help="Maximum number of files to process in one run.")
args = parser.parse_args()
# 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)
print("Nothing to do (give at least one of -p, -n)")
# sys.exit(0)
# 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)
print("%s already exists, exiting" % lock_file)
sys.exit(0)
# file(lock_file, 'w').write(pid)
# atexit.register(at_exit, lock_file)
open(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 = IASI_PATH
if args.old_iasi:
fs = path+"iasi_*bufr"
else:
fs = path+"W_XX-*_eps_o_l1.bin"
# 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]]
# 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']
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
if args.verbose:
print("Reading BUFR files", fns)
ii=0
processed_count = 0
# Set upper limit of files to process
max_process = 150
max_process = np.minimum(len(fns), max_process)
max_files = int(args.max_files)
max_files = np.minimum(len(fns), max_files)
# 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')
# while processed_count < max_files:
while processed_count < max_files and ii < len(fns):
fn = fns[ii]
# if fn.split('/')[-1] in invalid_files:
# continue
# Plot filename
if args.write_png is True:
fnbits = fn.split('/')
fnpng = fnbits[len(fnbits)-1]
fnpng = PNG_OUT + fnpng.replace(fn.split('.')[-1],'png')
else:
fnpng = None
# ncf filename
fnbits = fn.split('/')
fnncf = fnbits[len(fnbits)-1]
fnncf = NETCDF_OUT + fnncf.replace(fn.split('.')[-1],'nc')
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 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
# Skip if input file less limit (some files seems to cantain only 1 scan line
# Skip if input file less limit (some files seems to contain 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