Adapt the reading of the flexpart footprints to my case. Hard-coded fix for...

Adapt the reading of the flexpart footprints to my case. Hard-coded fix for now, needs to make it flexible.

pycif/plugins/models/flexpart/utils/mod_flexpart.f90

@@ -44,7 +44,7 @@ subroutine read_header(filename, numpoint, ibdate, ibtime, releases, &
 
 
 !  character(len=256)                           :: flexversion
-  character(len=29)                           :: flexversion
+  character(len=13)                           :: flexversion
   logical                                      :: lexist
   integer                                      :: nageclass, imethod, nspec
   integer                                      :: loutstep, loutaver, loutsample

pycif/plugins/models/flexpart/utils/read.py

@@ -8,6 +8,7 @@ import numpy as np
 import glob
 import re
 from datetime import datetime, timedelta
+from netCDF4 import Dataset
 
 from . import flexpart_header
 from . import mod_flexpart
@@ -75,25 +76,33 @@ def read_flexpart_grid(subdir, file_name, fp_header, **kwargs):
     file_dates = os.path.join(subdir, 'dates')
     grid = np.ndarray((fp_header.numx, fp_header.numy, fp_header.maxngrid))
     path_file = os.path.join(subdir, file_name)
+
+
+    #HJM
+    if False:
+        grid_fp, ngrid, gtime = mod_flexpart.read_grid(
+            path_file, file_dates, fp_header.numx, fp_header.numy, \
+            fp_header.maxngrid, fp_header.xshift,  fp_header.ndgrid, fp_header.trajdays)
+
+        # Convert grid times to datetime format
+        gtime_dt = []
+
+        # Convert FLEXPART julian days to datetime
+        for i in range(len(gtime)):
+            if gtime[i] == 0.:
+                break
+            gtime_dt.append(j2d(gtime[i]))
+    else:
+        print('HJM', path_file, file_dates, fp_header.numx, fp_header.numy,
+              fp_header.maxngrid, fp_header.xshift,  fp_header.ndgrid, fp_header.trajdays)
+        grid_fp, ngrid, gtime_dt = read_grid_nc(path_file, file_dates, fp_header)
     
-    grid_fp, ngrid, gtime = mod_flexpart.read_grid(
-        path_file, file_dates, fp_header.numx, fp_header.numy, \
-        fp_header.maxngrid, fp_header.xshift,  fp_header.ndgrid, fp_header.trajdays)
 
     # Convert from ppt to ppmv or ppbv
 
     # Convert s.m3/kg to s.m2/kg and apply numerical scaling
     grid_fp = grid_fp/(fp_header.outheight[0]*numscale)
 
-    # Convert grid times to datetime format
-    gtime_dt = []
-
-
-    # Convert FLEXPART julian days to datetime
-    for i in range(len(gtime)):
-        if gtime[i] == 0.:
-            break
-        gtime_dt.append(j2d(gtime[i]))
 
     # Generate the footprint dates as a series (for debugging), adding one hour
     # gtime_dt.append(datetime.strptime(re.findall(r'\d+', file_name)[0], '%Y%m%d%H%M%S'))
@@ -108,6 +117,31 @@ def read_flexpart_grid(subdir, file_name, fp_header, **kwargs):
 
     return grid_fp, gtime_dt, ngrid
 
+def read_grid_nc(path_file, file_dates, fp_header):
+    print('HJM reading nc')
+    # grid_fp, ngrid, gtime = mod_flexpart.read_grid(
+    #     path_file, file_dates, fp_header.numx, fp_header.numy, \
+    #     fp_header.maxngrid, fp_header.xshift,  fp_header.ndgrid, fp_header.trajdays)
+    with Dataset(path_file) as nc:
+        #spec001 has dimension nageclass, numpoint, time, level, lat, lon
+        # Numpoint determines the station 
+        #TODO add station pick
+        # For now, 0
+        grid = nc['spec001'][0,0,:,0,:,:] # time, lat, lon              
+        # swap axes to get it to lon, lat, time, 
+        grid_fp = np.swapaxes(grid,0,-1)
+
+        #iedate = "20120111"
+        times = np.sort(nc['time'][:])
+        print(times)
+        enddate = datetime.strptime(nc.getncattr('iedate'),'%Y%m%d')
+        gtime = []
+        for t in times:
+            gtime.append(enddate + timedelta(seconds=int(t)))
+    print('HJM done reading')
+
+    return grid_fp, len(gtime), gtime
+
 
 def read_init(subdir, file_name, nx, ny, nz, **kwargs):
     """ Reads FLEXPART grid_initial file.