Skip to content

GitLab

Commit 3dd554b5 authored by Espen Sollum's avatar Espen Sollum
Browse files

Merged in cganges from branch satdata

parent f63de215
Hide whitespace changes
Inline Side-by-side
bin/IASITools.py
View file @ 3dd554b5
......@@ -43,7 +43,6 @@ import matplotlib.pyplot as plt
from matplotlib.legend import rcParams
from matplotlib.ticker import MultipleLocator, FormatStrFormatter
import math
import pdb
# Note: in readbufr,
# set "old_iasi" to False if reading "older" IASI files (file names of form
......@@ -51,6 +50,11 @@ import pdb
# Newer files look like
# W_XX-EUMETSAT-Darmstadt,SOUNDING+SATELLITE,METOPB+IASI_C_EUMP_20191116020259_37157_eps_o_l1.bin
# NOTE #2: There is a hardcoded relative path (was ../data/SEVIRI_spectral_response/) that
# was changed to absolute. You need to manually edit this path if running the program
# from a location where this path is unavailable
SEVIRI_SR_PATH = "/viper/nadir/home/satdata/ash-iasi/git/ash-iasi/data/SEVIRI_spectral_response/"
### Map resolution for plots
map_res='l'
......@@ -120,7 +124,7 @@ class Scene(object):
fo.close()
return
def ReadBufr(self,fn,test=False,verbose=False, old_iasi=None):
def ReadBufr(self,fn,test=False,verbose=False, old_iasi=False):
#
# IASI data are given for latitude/longitude where latitude goes from
# -90$^{\circ}$ (S) to +90$^{\circ}$ (N), and longitude covers
......@@ -251,7 +255,10 @@ class Scene(object):
while True:
if verbose:
print("message: %s" % cnt)
bfr = bufr.codes_bufr_new_from_file(f)
# The program will usually hang here if bufr file is corrupt
try:
bfr = bufr.codes_bufr_new_from_file(f)
if bfr is None:
break
bufr.codes_set(bfr,'unpack',1)
......@@ -290,8 +297,11 @@ class Scene(object):
swath_tmp = (bufr.codes_get_array(bfr,"scaledIasiRadiance"))
swath.append(swath_tmp[0:self.Nwvls*self.Nspectra])
cnt += 1
# delete handle
# delete handle
bufr.codes_release(bfr)
except Exception as e:
print("ERROR reading file ", fn)
return -1
......@@ -386,7 +396,6 @@ class Scene(object):
return 0
def ReadFlexpartNetcdf(self,fn,test=False,verbose=False):
if verbose:
......@@ -428,11 +437,16 @@ class Scene(object):
def CalculateSEVIRIChannels(self, Channels=['10.8','12.0']):
self.SEVIRI = np.zeros((self.Nscanlines*self.Nspectra, len(Channels)))
print(self.Rad.shape, self.SEVIRI.shape)
## print(self.Rad.shape, self.SEVIRI.shape)
self.Channels=Channels
ic=0
for Channel in Channels:
fns='../data/SEVIRI_spectral_response/SEVIRI_'+Channel+'_wvl.dat'
# fns='../data/SEVIRI_spectral_response/SEVIRI_'+Channel+'_wvl.dat'
SEVIRI_SR_PATH
fns=SEVIRI_SR_PATH + 'SEVIRI_' +Channel+ '_wvl.dat'
if not os.path.isfile(fns):
print("ERROR: file not found: ", fns)
sys.exit(1)
SpectralResponseData=np.loadtxt(fns)
SpectralResponseWvn=SpectralResponseData[:,0]
SpectralResponse =SpectralResponseData[:,1]
......@@ -449,7 +463,7 @@ class Scene(object):
def CalculateAVHRRChannels(self, Channels=['4','5']):
self.AVHRR = np.zeros((self.Nscanlines*self.Nspectra, len(Channels)))
print(self.Rad.shape, self.SEVIRI.shape)
## print(self.Rad.shape, self.SEVIRI.shape)
self.Channels=Channels
ic=0
for Channel in Channels:
......@@ -709,7 +723,7 @@ class Scene(object):
ndig=-1
ymax=round(math.ceil(max(nxs)/div)*div,base)
ymaL=round(ymax/8,ndig)
print(ymax, ymaL)
##print(ymax, ymaL)
ymiL=ymaL/10
ymin=0
xlabel='Ash mass loading (g/m$^2$)'
......@@ -752,7 +766,7 @@ class Scene(object):
fig=plt.figure(figsize=(fact*xsize, fact*ysize))
ax = plt.subplot(111)
rcParams.update({'font.size':fsize})
print(XType, YType)
##print(XType, YType)
if XType=='IWC' and YType=='AshML':
ind = np.where(self.typenr==3)
data = self.IceIWC
......@@ -983,7 +997,8 @@ class Scene(object):
m.drawparallels(parallels,labels=[1,0,0,0],ax=ax)
m.drawmeridians(meridians,labels=[0,0,0,1],ax=ax)
x,y = m(self.lons,self.lats)
cmap = 'jet'#'gist_ncar_r'#
# cmap = 'jet'#'gist_ncar_r'#
cmap = plt.get_cmap('jet') #'gist_ncar_r'#
if PlotAshContour:
# dBT = self.BT[:,self.ind1231_5]-self.BT[:,self.ind1168]
......@@ -1022,7 +1037,8 @@ class Scene(object):
# vmax= 0.05*1000 # 2 # 20 #data.max()*0.5#
step=0.08 # vmax/5. #0.01
cblabel='Ash Mass Loading (g/m$^2$)'
cmap = 'gist_ncar_r'#
# cmap = 'gist_ncar_r'#
cmap = plt.get_cmap('gist_ncar_r')
print("GABBA", data.shape, self.lats.shape, self.lons.shape)
elif Type=='AshRe':
data = self.AshRe[:]
......@@ -1031,6 +1047,7 @@ class Scene(object):
step=1.0
cblabel='Ash Effective Radius ($\mu$m)'
cmap = 'gist_ncar_r'#
cmap = plt.get_cmap('gist_ncar_r')
print("GABBA", data.shape, self.lats.shape, self.lons.shape)
elif Type=='IceIWC':
data = self.IceIWC[:]
......@@ -1039,12 +1056,13 @@ class Scene(object):
step=0.08 #vmax/5 #0.005
cblabel='Ice water content (g/m$^3$)'
cmap = 'gist_ncar_r'#
cmap = plt.get_cmap('gist_ncar_r')
elif Type=='IceCloudTop':
data = self.IceCloudTop[:]
ind = np.where( self.IceIWC > 0.01)
print(data.shape)
##print(data.shape)
data = data[ind]
print(data.shape)
##print(data.shape)
self.lats = self.lats[ind]
self.lons = self.lons[ind]
self.ascan = self.ascan[ind]
......@@ -1056,6 +1074,8 @@ class Scene(object):
cblabel='Ice cloud top (km)'
cmap = 'jet'#
cmap = 'gist_ncar_r'#
cmap = plt.get_cmap('gist_ncar_r')
elif Type=='Pearson':
data = self.corr[:]
vmin=0.85
......@@ -1155,6 +1175,7 @@ class Scene(object):
norm = mpl.colors.Normalize(vmin=vmin, vmax=vmax)
# cmap = cm.jet
cmap = plt.get_cmap('jet')
ms = cm.ScalarMappable(norm=norm, cmap=cmap)
ip=0
Ntot, = data.shape #self.Nspectra*self.Nscanlines
......@@ -1242,6 +1263,7 @@ class Scene(object):
rcParams.update({'font.size':fsize})
if PlotFile==None:
print("No PlotFile given")
if PlotShow:
# plt.show()
sys.exit()
......@@ -1251,6 +1273,8 @@ class Scene(object):
plt.savefig(PlotFile)
plt.close('all')
print("Writing plot to ", PlotFile)
def ReadNetCDF(self,fn,test=False,verbose=False, AVHRRinfo=False):
"""
......@@ -1271,61 +1295,63 @@ class Scene(object):
if verbose:
print("ReadNetCDF: fn:", fn)
if (type(fn) != list):
ncfile = Dataset(fn,'r')
if verbose:
print(ncfile.variables)
self.lats = ncfile.variables['lat'][:]
self.lons = ncfile.variables['lon'][:]
self.ascan = ncfile.variables['ascan'][:]
self.atrack = ncfile.variables['atrack'][:]
self.dBT = ncfile.variables['dBT'][:]
self.npoints = len(self.lats)
# These variables are possibly missing
try:
self.SEVIRI_dBT = ncfile.variables['SEVIRI_dBT'][:]
self.AshML = ncfile.variables['ML'][:]
self.AshRe = ncfile.variables['Re'][:]
self.corr = ncfile.variables['corr'][:]
self.rms = ncfile.variables['rms'][:]
self.IceIWC = ncfile.variables['IceIWC'][:]
self.IceCloudTop = ncfile.variables['IceCloudTop'][:]
self.typenr = ncfile.variables['typenr'][:]
except:
pass
if AVHRRinfo:
self.AVHRRBT4std = ncfile.variables['AVHRRBT4std'][:]
ncfile.close()
else:
for i,f in enumerate(fn):
ncfile = Dataset(f, 'r')
if i == 0:
self.lats = ncfile.variables['lat'][:]
self.lons = ncfile.variables['lon'][:]
self.ascan = ncfile.variables['ascan'][:]
self.atrack = ncfile.variables['atrack'][:]
self.dBT = ncfile.variables['dBT'][:]
else:
self.lats = np.concatenate((self.lats, ncfile.variables['lat'][:]))
self.lons = np.concatenate((self.lons, ncfile.variables['lon'][:]))
self.ascan = np.concatenate((self.ascan, ncfile.variables['ascan'][:]))
self.atrack = np.concatenate((self.atrack, ncfile.variables['atrack'][:]))
self.dBT = np.concatenate((self.dBT,ncfile.variables['dBT'][:]))
ncfile.close()
self.npoints = len(self.lats)
return
def WriteNetcdf(self,fn, AVHRRinfo=False ):
try:
if (type(fn) != list):
ncfile = Dataset(fn,'r')
if verbose:
print(ncfile.variables)
self.lats = ncfile.variables['lat'][:]
self.lons = ncfile.variables['lon'][:]
self.ascan = ncfile.variables['ascan'][:]
self.atrack = ncfile.variables['atrack'][:]
self.dBT = ncfile.variables['dBT'][:]
self.npoints = len(self.lats)
# These variables are possibly missing
try:
self.SEVIRI_dBT = ncfile.variables['SEVIRI_dBT'][:]
self.AshML = ncfile.variables['ML'][:]
self.AshRe = ncfile.variables['Re'][:]
self.corr = ncfile.variables['corr'][:]
self.rms = ncfile.variables['rms'][:]
self.IceIWC = ncfile.variables['IceIWC'][:]
self.IceCloudTop = ncfile.variables['IceCloudTop'][:]
self.typenr = ncfile.variables['typenr'][:]
except:
pass
if AVHRRinfo:
self.AVHRRBT4std = ncfile.variables['AVHRRBT4std'][:]
ncfile.close()
else:
for i,f in enumerate(fn):
ncfile = Dataset(f, 'r')
if i == 0:
self.lats = ncfile.variables['lat'][:]
self.lons = ncfile.variables['lon'][:]
self.ascan = ncfile.variables['ascan'][:]
self.atrack = ncfile.variables['atrack'][:]
self.dBT = ncfile.variables['dBT'][:]
else:
self.lats = np.concatenate((self.lats, ncfile.variables['lat'][:]))
self.lons = np.concatenate((self.lons, ncfile.variables['lon'][:]))
self.ascan = np.concatenate((self.ascan, ncfile.variables['ascan'][:]))
self.atrack = np.concatenate((self.atrack, ncfile.variables['atrack'][:]))
self.dBT = np.concatenate((self.dBT,ncfile.variables['dBT'][:]))
ncfile.close()
self.npoints = len(self.lats)
except Exception as e:
print("Exception: ", e)
return
def WriteNetcdf(self,fn, AVHRRinfo=False, verbose=False ):
# According to CF all names should be lower case
ncf = Dataset(fn,'w')
......@@ -1526,6 +1552,8 @@ class Scene(object):
data[np.abs(data) > 999] = -999.9
SEVIRI_dBT[:] = data
if verbose:
print("Write netcdf file", fn)
ncf.close()
......
bin/PlotSelectedNCF.py
View file @ 3dd554b5
#!/usr/bin/env python3
#!/usr/bin/env python3.7
"""
Plot ash detections for selected IASI (processed) netcdf files, selected
......@@ -26,6 +26,18 @@ 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,
......@@ -38,7 +50,10 @@ def PlotIASIFilteredGeo(nc_directory,
Projection=None,
start_date=None,
end_date=None,
link_directory=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
......@@ -70,7 +85,7 @@ def PlotIASIFilteredGeo(nc_directory,
if start_date==None:
# search 1 day back in time by default
start_date=dt.datetime.utcnow() - timedelta(days=1)
start_date=dt.datetime.utcnow() - timedelta(days=days_to_process)
else:
start_date = dt.datetime.strptime(start_date, '%Y%m%d%H%M%S')
......@@ -130,6 +145,11 @@ def PlotIASIFilteredGeo(nc_directory,
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',
......@@ -169,11 +189,11 @@ if __name__ == "__main__":
dest='plot_multi',
help="Concatenate a list of netcdf files and combine in one plot")
parser.add_argument(
'--netcdf_directory', default=None,
'--netcdf_directory', default="/viper/nadir/nvap/IASI/",
dest='netcdf_directory',
help="")
parser.add_argument(
'--plot_directory', default=None,
'--plot_directory', default="/viper/nadir/nvap/IASI/plots/",
dest='plot_directory',
help="")
parser.add_argument(
......@@ -200,6 +220,10 @@ if __name__ == "__main__":
'--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',
......@@ -216,6 +240,15 @@ if __name__ == "__main__":
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
......@@ -249,7 +282,9 @@ if __name__ == "__main__":
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)
old_iasi=args.old_iasi,link_directory=args.link_directory,
verbose=args.verbose,
days_to_process=int(args.days_to_process))
......
bin/ReadAndPlotIASI.py
View file @ 3dd554b5
#!/usr/bin/env python3
#!/usr/bin/env /usr/local/bin/python3.7
"""
Produce netCDF files of IASI ash detections for met.no
......@@ -17,6 +17,7 @@ import argparse
import atexit
import numpy as np
import importlib
import datetime as dt
importlib.reload(IASITools)
......@@ -31,14 +32,17 @@ importlib.reload(IASITools)
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/"
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/"
# 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):
......@@ -54,7 +58,7 @@ def at_exit(file_name):
def ReadAndPlotSingleIASIFile(
fn, Projection='lcc', PlotFile=None, ncfFile=None, verbose=None, old_iasi=None):
fn, Projection='lcc', PlotFile=None, ncfFile=None, verbose=None, old_iasi=False):
Scene = IASITools.Scene()
ret = Scene.ReadBufr(fn, test=False, verbose=verbose, old_iasi=old_iasi)#True)#
......@@ -84,7 +88,7 @@ def ReadAndPlotSingleIASIFile(
PlotAshContour=False)
Scene.WriteNetcdf(ncfFile)
Scene.WriteNetcdf(ncfFile, verbose=verbose)
return Scene
......@@ -143,10 +147,10 @@ if __name__ == "__main__":
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)
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)
......@@ -173,8 +177,8 @@ if __name__ == "__main__":
'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
......@@ -191,10 +195,7 @@ if __name__ == "__main__":
if args.write_png is True:
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')
else:
fnpng = None
......@@ -213,16 +214,25 @@ if __name__ == "__main__":
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?)
# Tka into account that the file may have been deleted since start of processing
try:
if os.path.getsize(fn) < 2.9*1.e7:
print("Skipping file (likely corrupted due to file size):", fn)
ii = ii + 1
continue
# except FileNotFoundError as f_err:
except Exception as e:
continue
# Add this clause since file may have been deleted since the launch of the program
try:
print("processing file: ", fn)
if args.verbose:
print(dt.datetime.utcnow().strftime("%Y-%m-%d %H:%M"),": processing file ", fn)
lonmin, lonmax, latmin, latmax = IASITools.GetLatLonRegion(fn, verbose=args.verbose)
if args.verbose:
print(lonmin, lonmax, latmin, latmax, fnpng)
......@@ -239,11 +249,12 @@ if __name__ == "__main__":
if Scene == None:
print("Invalid Scene")
continue
except:
except Exception as e:
print("Error", e)
ii = ii + 1
continue
# exit()
# Remove lock file
# DEV: this will be removed by atexit?
# if os.path.isfile(lock_file):
# os.unlink(lock_file)
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment