Added Python script for combining swaths and adding netcdf flags

bin/IASITools.py

@@ -386,6 +386,7 @@ class Scene(object):
 
         return 0
 
+    
     def ReadFlexpartNetcdf(self,fn,test=False,verbose=False):
 
         if verbose:

bin/PlotSelectedNCF.py

@@ -26,6 +26,8 @@ import importlib
 importlib.reload(IASITools)
 
 
+# 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, 
@@ -96,8 +98,6 @@ def PlotIASIFilteredGeo(nc_directory,
 
         file_date = dt.datetime.strptime(date_str, '%Y%m%d%H%M%S')
 
-#        pdb.set_trace()
-
 ####################################################################
 ### Check if swath within specified time limits
 ####################################################################
@@ -124,15 +124,12 @@ def PlotIASIFilteredGeo(nc_directory,
 
  #       print s.center_lon>lon_w, s.center_lon<lon_e, s.center_lat>lat_s, s.center_lat<lat_n
         
-#        pdb.set_trace()
-
         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()
             s.PlotLatLonVals(Title=False,
                              ManualTitle=man_title,
                              Type='dBT_already_calculated', 
@@ -241,13 +238,13 @@ if __name__ == "__main__":
 ### Plot all files in netcdf directory that have swath center within coordinate
 ### bounds
         if lon_min is None:
-            lon_min = 155.
+            lon_min = -30.
         if lon_max is None:            
-            lon_max = 180.
+            lon_max = 30.
         if lat_min is None:
-            lat_min = 50.
+            lat_min = 40.
         if lat_max is None:
-            lat_max = 65.
+            lat_max = 80.
         
         PlotIASIFilteredGeo(nc_directory, lon_min, lon_max, lat_min, lat_max,
                             plot_directory, start_date=start_date,

bin/ReadAndPlotIASI.py

@@ -151,7 +151,7 @@ if __name__ == "__main__":
             
 
 # Start by clearing the netcdf output directory for files older that 14 days
-    IASITools.remove_old_files(NETCDF_OUT, 10000, False, dry_run=False)
+    IASITools.remove_old_files(NETCDF_OUT, 14, False, dry_run=False)
 
     path = IASI_PATH
 
@@ -165,7 +165,7 @@ if __name__ == "__main__":
     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 (..."

bin/TestIASI.py

@@ -3,7 +3,7 @@
 import ReadAndPlotIASI
 
 # New file format
-#ReadAndPlotIASI.ReadAndPlotSingleIASIFile('W_XX-EUMETSAT-Darmstadt,SOUNDING+SATELLITE,METOPB+IASI_C_EUMP_20191116020259_37157_eps_o_l1.bin', ncfFile='test.nc')
+ReadAndPlotIASI.ReadAndPlotSingleIASIFile('W_XX-EUMETSAT-Darmstadt,SOUNDING+SATELLITE,METOPB+IASI_C_EUMP_20191116020259_37157_eps_o_l1.bin', ncfFile='test.nc')
 
 # Older file format, tends to not work...
-ReadAndPlotIASI.ReadAndPlotSingleIASIFile('/xnilu_wrk/projects/sat_data/events/2010-04_Eyjafjallajokull/eumetcast/EPS-11/iasi_20100417_214753_metopa_18131_eps_o.l1_bufr', ncfFile='test_old.nc', old_iasi=True, PlotFile='test.png')
+#ReadAndPlotIASI.ReadAndPlotSingleIASIFile('/xnilu_wrk/projects/sat_data/events/2010-04_Eyjafjallajokull/eumetcast/EPS-11/iasi_20100417_214753_metopa_18131_eps_o.l1_bufr', ncfFile='test_old.nc', old_iasi=True, PlotFile='test.png')

bin/makeCombined.py

+#!/usr/bin/env python3
+"""
+DESCRIPTION
+    Combine IASI netcdf files (created by `PlotSelectedNCF.py´) into 
+    6 hour aggregate files.
+
+    Add ash flags, and set/edit various netCDF attributes
+
+NOTE
+    Uses the "non-standard" pynco library, install under your user directory with
+    `python3 -m pip install --user nco´
+
+TODO
+    A cron job deletes IASI files older than n (=14?) days. Take care in this script 
+    of the possibility that files may disappear while the script is running.
+
+    Should probably have a "running" mode where last 2-3 days are (re)processed,
+    and a "force" mode to reprocess everything
+
+
+"""
+
+import os
+import argparse
+import glob
+import sys
+import re
+import shutil
+import numpy as np
+from datetime import datetime, timedelta
+from nco import Nco
+
+# Parameters, paths etc.
+# TMPDIR = "/tmp/"
+TMPDIR = "./tmp/"
+
+# Directory where input files are found
+INDIR = "." # For test
+
+# Directory where merged files will be written
+OUTDIR = TMPDIR
+
+# Time interval for merged files (6h). Should probably be integer divisor of 24h
+DT = timedelta(hours=6)
+
+
+
+def run_test(idir):
+    """
+    Run a test, using two files
+
+    """
+
+    fn1 = "iasi_20100417_213857_metopa_18131_eps_o.nc"
+    fn2 = "iasi_20100417_213553_metopa_18131_eps_o.nc"
+
+    print("Running test")
+
+    file_list = [fn1, fn2]
+    time_stamp = datetime.utcnow().strftime("%Y%m%d_%H%M")
+    
+    add_parameter(file_list, time_stamp, idir)
+
+
+def make_file_list(in_directory):
+    """
+    Return list of IASI nc files to be merged
+    
+    Finds dates of oldest/newest IASI files, return dictionary 
+    of the form {time inteval : [file_list]}
+    
+
+    """
+    
+    in_files = [f for f in os.listdir(in_directory) if (f.endswith('.nc') and
+                                                        f.startswith('W_XX'))]
+
+    # Note: for older files (e.g. 2010) use this:
+                                                        # f.startswith('iasi_'))]
+
+    # Make a dictionary with time as key and file name as value
+    files_dict = {}
+
+    for f in in_files:
+        files_dict.update({datetime_iasi(f) : f})
+
+    return files_dict
+
+
+def make_interval(files_dict, days_to_process=None):
+    """
+    Find which files to be used for each 6hr interval.
+    Stored as dict with time interval (str) as key, files as list of values
+
+    """
+
+    earliest_date = min(files_dict)
+    latest_date = max(files_dict)
+
+
+    # Start processing intervals beginning at hour 0 of oldest day
+    start_day = earliest_date.replace(hour=0, minute=0)
+    end_day = latest_date.replace(hour=0, minute=0) + timedelta(days=1)
+
+    # If days_to_process is given, only process the given days back in time
+    if days_to_process is not None:
+        start_day = end_day - timedelta(days=days_to_process)
+
+    start_interval = start_day
+
+    int_dict = {}
+
+    # Loop over all files and put files in correct time interval
+    while start_interval < end_day:
+        end_interval = start_interval + DT
+
+        int_files = []
+
+        for key, value in files_dict.items():
+            if key > start_interval and key < end_interval:
+                int_files.append(value)
+                
+        int_dict.update({dt_to_name(start_interval, end_interval) : int_files})
+
+        start_interval += DT
+
+    return int_dict
+
+
+def dt_to_name(start, end):
+    """
+    Return a string from two given datetimes.
+    Used for dict keys (and for file names)
+
+    """
+
+    return start.strftime("%Y%m%d_%H%M") + end.strftime("-%H%M") 
+
+
+def datetime_iasi(file_name):
+    """
+    Extract date from IASI file name. 
+    Return as datetime object
+
+    """
+
+    # For older IASI files:
+    # DT_FMT = "%Y%m%d_%H%M"
+    # date_str = re.findall(r"\d{8}_\d{4}", file_name)
+
+    DT_FMT = "%Y%m%d%H%M"
+    date_str = re.findall(r"\d{12}", file_name)
+        
+
+
+    if len(date_str) != 1:
+        print ("<datetime_iasi> Error: could not extract date from filename")
+        raise Exception
+
+    return datetime.strptime(date_str[0], DT_FMT)
+    
+
+def add_parameter(file_list, time_stamp, idir):
+    """
+    Add ncf flags.
+    Merge files over time period
+    
+    """
+
+    # Check input args
+    if len(file_list) == 0:
+        print("<add_parameter> Warning: empty file list")
+        return
+
+    # nco instance
+    nco = Nco()
+
+    # Delete tmp files, if existing
+    files_in = glob.glob(os.path.join(TMPDIR, "IASI_in*.nc"))
+    files_out = glob.glob(os.path.join(TMPDIR, "IASI_out*.nc"))
+
+    for f in files_in:
+        try:
+            os.remove(f)
+        except:
+            print("Error removing file", f)
+
+    for f in files_out:
+        try:
+            os.remove(f)
+        except:
+            print("Error removing file", f)
+                         
+    for i, f in enumerate(file_list):
+        in_tmp = os.path.join(TMPDIR, "IASI_in_" + str(i) + ".nc")
+        out_tmp = os.path.join(TMPDIR, "IASI_out_" + str(i) + ".nc")
+
+        try:
+            shutil.copyfile(os.path.join(args.idir, f), in_tmp)
+        except FileNotFoundError as f_err:
+            continue
+        
+        nco.ncks(input=in_tmp, output=out_tmp, options=["-O", "--mk_rec_dmn", "npoints" ])
+
+        nco.ncap2(input=out_tmp, output=in_tmp, options=["-O", "-s", "'time=int(dBT*0+060000);time@long_name=\"Time(hhmmss)\";time@standard_name=\"hhmmss\";time@units=\"hhmmss\"'"])
+
+        nco.ncap2(input=in_tmp, output=out_tmp, options=["-O", "-s", "'flag=int(dBT*0+1);flag@long_name=\"flag\";flag@standard_name=\"flag\";flag@units=\"2=ash,1=uncertain,0=no ash\" '"])
+
+        nco.ncap2(input=out_tmp, output=in_tmp, options=["-O", "-s", "'where(dBT<2.0) flag=0;'"])
+
+        nco.ncap2(input=in_tmp, output=out_tmp, options=["-O", "-s", "'where(dBT>3.0) flag=2;'"])
+                 
+    # Merge files
+    in_files = [os.path.join(TMPDIR, "IASI_out_" + str(n) +".nc") for n in range(i+1)]
+
+    # Make output file name
+    file_out = "Combined_IASI_" + time_stamp + ".nc"
+    nco.ncrcat(input=in_files, output=os.path.join(TMPDIR, file_out), options=["-O", "-h"])
+
+
+    # Delete tmp files on exit
+    files_in = glob.glob(os.path.join(TMPDIR, "IASI_in*.nc"))
+    files_out = glob.glob(os.path.join(TMPDIR, "IASI_out*.nc"))
+
+    for f in files_in:
+        try:
+            os.remove(f)
+        except:
+            print("Error removing file", f)
+
+    for f in files_out:
+        try:
+            os.remove(f)
+        except:
+            print("Error removing file", f)
+                  
+
+def run_program(in_dir, out_dir, days_tp):
+    """
+    Process all files in in_dir, wrirte output to out_dir
+
+    """
+
+    # Get all dates and file names
+    files_dict = make_file_list(in_dir)
+    
+    # Find files within each 6h interval
+    interval_files = make_interval(files_dict, days_to_process=int(days_tp))
+
+    # Merge files and add ncf parameters
+    for key in interval_files:
+        if len(interval_files[key]) > 0:
+            add_parameter(interval_files[key], key, in_dir)
+            
+
+        
+if __name__ == "__main__":      
+# Program options
+    parser = argparse.ArgumentParser(prog="makeCombined", 
+                                     description=":DESCRIPTION:",
+                                     formatter_class=argparse.
+                                     ArgumentDefaultsHelpFormatter)
+        
+    parser.add_argument('-v', '--verbose', 
+                        default=False,
+                        action='store_true',
+                        dest='verbose',
+                        help="verbose output")
+
+    parser.add_argument('--test', 
+                        default=False,
+                        action='store_true',
+                        dest='test',
+                        help="If true, run a test")
+
+    parser.add_argument('--run', 
+                        default=False,
+                        action='store_true',
+                        dest='run',
+                        help="If true, run for all files in input directory")
+
+    parser.add_argument('--idir', 
+                        default=INDIR,
+                        dest='idir',
+                        help="Location of IASI netcdf file to process. Default: " + INDIR)
+
+    parser.add_argument('--odir', 
+                        default=OUTDIR,
+                        dest='odir',
+                        help="Where to write processed files. Default: " + OUTDIR)
+
+    parser.add_argument('--days_to_process', 
+                        default=None,
+                        dest='days_to_process',
+                        help="Days back in time to process. If not given, process all "
+                        "files found on disk (typically 14 days)")
+
+
+    args = parser.parse_args()
+
+
+    # Check input arguments
+
+    # Lock file handling (avoid more than 1 running instance)
+
+
+    if args.test:
+        run_test(args.idir)
+    elif args.run:
+        run_program(args.idir, args.odir, args.days_to_process)
+    else:
+        print("Nothing to do")
+        sys.exit(0)
+    
+    

bin/make_IASI_ash_12h.sh

+#!/bin/bash
+
+#file1='IASI_'
+#procdir=/viper/nadir/tmp/IASI/
+
+#file1='W_XX-EUMETSAT-Darmstadt,SOUNDING+SATELLITE,METOP?+IASI_C_EUM?_'
+#procdir=/xnilu_wrk/users/eso/IASI/netcdf/
+procdir=/home/eso/repos/ash-iasi/data/links/
+file1='iasi_'
+
+###### Loop over year, month, day and hour ######
+
+for year in $(seq -f "%04g" 2010 2010); do
+for month in $(seq -f "%02g" 04 04); do
+for dd in $(seq -f "%02g" 16 17); do
+ls $procdir$file1$year$month$dd* > file.txt
+for hind in 00 12; do
+
+let H2=$hind+12
+H2=$(printf %02g $H2)
+
+new_file=Combined_IASI_$year$month$dd'_'$hind'00-'$H2'00.nc'
+touch $new_file
+
+i=0
+while read line
+do
+  echo 
+  echo $line
+
+#IASI directory
+#  HHMMname=${line:92:6}
+#  HH=${line:92:2}
+#Espens home
+  HHMMname=${line:47:6}
+  HH=${line:47:2}
+  echo $HHMMname $HH
+  echo $hind
+ 
+  calc=0 
+  if [ $hind -eq 00 ] && [ $HH -lt 12 ]; then
+     calc=1
+  elif [ $hind -eq 12 ] && [ $HH -ge 12 ] && [ $HH -le 24 ]; then
+     calc=1
+  fi  
+  echo $calc
+
+
+  if [ $calc -eq 1 ]; then
+
+  let i=i+1 
+
+  #Add dimension (NB! ncks do not accept file names with comma)
+  cp $line tmp2.nc
+  ncks -O --mk_rec_dmn npoints tmp2.nc tmp1.nc
+
+  echo $i
+
+  #Add parameter:
+  ncap2 -O -s 'time=int(dBT*0+'$HHMMname');time@long_name="Time(hhmmss)";time@standard_name="hhmmss";time@units="hhmmss"' tmp1.nc tmp2.nc
+  ncap2 -O -s 'flag=int(dBT*0+1);flag@long_name="flag";flag@standard_name="flag";flag@units="2=ash,1=uncertain,0=no ash"' tmp2.nc tmp1.nc
+  ncap2 -O -s 'where(dBT<2.0) flag=0;' tmp1.nc tmp2.nc
+  ncap2 -O -s 'where(dBT>3.0) flag=2;' tmp2.nc FILE_$i.nc
+
+  ## End file within time frame (calc=1)
+  fi
+
+done < file.txt
+
+#Merge files, delete long list of global attributes, rename
+echo COMBINE FILE: $new_file
+echo Number of files: $i
+
+if [ $i -gt 0 ]; then
+  ncrcat -O -h FILE_?.nc final.nc 
+  if [ $i -ge 10 ]; then
+     mv final.nc tmp1.nc
+     ncrcat -O -h FILE_??.nc tmp2.nc
+     ncrcat -O -h tmp1.nc tmp2.nc final.nc 
+  fi
+  ncatted -O -h -a ash_if_dBT_above,global,c,c,3.0 final.nc
+  ncatted -O -h -a no_ash_if_dBT_below,global,c,c,2.0 final.nc
+  ncatted -O -h -a history,global,d,, final.nc $new_file
+  rm -f tmp1.nc tmp2.nc final.nc FILE_*
+fi
+
+#Done: hour, Day, Month, Year
+done
+done
+done
+done

bin/plotCombined.py

@@ -20,33 +20,36 @@ import pdb
 import importlib
 importlib.reload(IASITools)
 
-# above is pasted from plot
-#import PlotSelectedNC
 
-
-def plotCombined(file_name):
+def plotCombined(file_name, ignore=False, Type=None):
     ## :TODO: make input arg
 #    Type='dBT_already_calculated'
-    Type='SEVIRI_dBT_already_calculated'
-    
-    if Type == 'SEVIRI_dBT_already_calculated':
+    if Type is None or Type == "SEVIRI_dBT":
+        plotmethod='SEVIRI_dBT_already_calculated'
         str1 = 'SEVIRI_dBT'
-    elif Type == 'dBT_already_calculated':
+    elif Type == "IASI_dBT":
+        plotmethod = 'SEVIRI_dBT_already_calculated'
         str1 = 'IASI_dBT'
     
+    if ignore:
+        date_str = dt.datetime.utcnow().strftime('%Y%m%d%H%M%S')
+        man_title = date_str
+        fnpng = os.path.basename(file_name).split('.nc')[0] + '_' + Type + ".png"
+    else:
+        _tmp = os.path.basename(file_name).split('_')
+        date_str = _tmp[2][0:4] + '-' + _tmp[2][4:6] + '-' + _tmp[2][6:8] + '-' + _tmp[3][0:3]
+        man_title = date_str
+        fnpng = os.path.basename(file_name).split('.nc')[0] + '_' + Type + ".png"
     
-    _tmp = os.path.basename(file_name).split('_')
-    date_str = _tmp[2][0:4] + '-' + _tmp[2][4:6] + '-' + _tmp[2][6:8] + '-' + _tmp[3][0:3]
-
     s = IASITools.Scene()
     s.ReadNetCDF(file_name)
 
-    man_title = date_str
-    fnpng = os.path.basename(file_name).split('.nc')[0] + '_' + str1 + ".png"
+
+
 
     s.PlotLatLonVals(Title=False,
                      ManualTitle=man_title,
-                     Type=Type, 
+                     Type=plotmethod, 
                      test=True,
                      PlotAshContour=True,
                      PlotFile=fnpng,
@@ -67,12 +70,24 @@ if __name__ == "__main__":
                         dest='path_name',
                         help="Path to IASI Netcdf files to plot")
 
+    parser.add_argument('-t', '--type', 
+                        dest='Type',
+                        default ='SEVIRI_dBT',
+                        help="plot type [IASI_dBT, SEVIRI_dBT]")
+
+    parser.add_argument('-i', '--ignore_date', 
+                        dest='ignore_date', action="store_true",
+                        default=False,
+                        help="Just plot the file, ignoring date and time info from file name")
+
+    
+
     args = parser.parse_args()
 
     # Plot file (skip if size 0)
     full_file_name = os.path.join(args.path_name, args.file_name) 
     if os.stat(full_file_name).st_size > 0:
-        plotCombined(full_file_name)
+        plotCombined(full_file_name, ignore=args.ignore_date, Type=args.Type)