""" Routines for reading Lagrangian footprints
"""
import os
import numpy as np
import pandas as pd
from scipy.io import FortranFile
import datetime
from logging import info, warning
from netCDF4 import Dataset
[docs]
def read_grid(path_file, release_date, fp_header, scaleconc):
"""Read a FLEXPART binary footprint grid file.
Parses the Fortran-unformatted binary format (``grid_time_*_001`` or
``grid_initial_*_001``) and reconstructs the sparse footprint array as
a dense ``(nlon, nlat, ntimes)`` array.
Args:
path_file (str): path to the binary footprint file.
release_date (pd.Timestamp): observation release date.
fp_header: FLEXPART domain header object (carries grid metadata).
scaleconc (float): unit scaling factor applied to the footprint.
Returns:
tuple: ``(grid_fp, gtime, ngrid, valid_file)`` where *grid_fp* is
the dense footprint array, *gtime* is the list of footprint dates,
*ngrid* is the time dimension, and *valid_file* is ``True`` if the
file was successfully read.
"""
days = []
times = []
counts_i = []
counts_r = []
sparse_i = []
sparse_r = []
valid_file = False
with FortranFile(path_file, 'r') as f:
# Looping until the end of the binary file
while True:
try:
yyyymmdd = f.read_ints('i4')[0].astype(str)
hhmmss = f"{f.read_ints('i4')[0]:06d}"
i = f.read_ints('i4')
spi = f.read_ints('i4')
r = f.read_ints('i4')
spr = f.read_reals('f4')
days.append(yyyymmdd)
times.append(hhmmss)
counts_i.extend(i)
sparse_i.extend(spi)
counts_r.extend(r)
sparse_r.extend(spr)
except TypeError:
break
gtime_dt = [datetime.datetime.strptime(f'{d}{h}', '%Y%m%d%H%M%S')
for d, h in zip(days, times)][::-1]
ngrid = len(gtime_dt)
grid_fp = np.zeros((fp_header.numx, fp_header.numy, ngrid + 2))
if len(sparse_r) > 0:
sign = np.sign(sparse_r)
cum_counts = np.unique(np.cumsum(counts_r) % sign.size)
signchange = ((np.roll(sign, 1) - sign) != 0)
signchange[cum_counts] = True
inds_out = np.zeros((len(sparse_r)), dtype=int)
inds_out[signchange] = sparse_i
mask = inds_out == 0
idx = np.where(~mask, np.arange(mask.size), 0)
np.maximum.accumulate(idx, out=idx)
inds_out = inds_out[idx] + np.arange(mask.size) \
- idx - fp_header.numx * fp_header.numy
try:
jy, jx = np.unravel_index(
inds_out, (fp_header.numy, fp_header.numx))
jt = np.zeros((len(sparse_r)), dtype=int)
jt[cum_counts] = np.arange(ngrid)[np.array(counts_r) > 0]
np.maximum.accumulate(jt, out=jt)
jt = ngrid - jt - 1
grid_fp[jx, jy, jt] = np.abs(sparse_r) * scaleconc
valid_file = True
except ValueError:
warning(f"Could not properly read file {path_file}. Returning zero footprints")
pass
grid_fp = np.roll(grid_fp, fp_header.xshift, axis=0)
return grid_fp, ngrid, gtime_dt, valid_file
[docs]
def read_grid_nc(path_file, release_date, fp_header):
"""Read a FLEXPART NetCDF footprint file.
Reads the ``spec001`` variable from the NetCDF footprint file and
returns a dense ``(nlon, nlat, ntimes)`` footprint array.
Args:
path_file (str): path to the NetCDF footprint file.
release_date (pd.Timestamp): observation release date.
fp_header: FLEXPART domain header object.
Returns:
tuple: ``(grid_fp, gtime, ngrid, valid_file)``.
"""
valid_file = False
with Dataset(path_file) as nc:
# spec001 has dimension nageclass, numpoint, time, level, lat, lon
# Assume that there is 1 numpoint/station per file
shape = nc['spec001'].shape
if shape[0] > 1 or shape[1] > 1:
info('WARNING: There are more than 1 station in the flexpart file.')
if shape[3] > 1:
info('INFO: Select the bottom layer of the grid')
grid = nc['spec001'][0, 0, :, 0, :, :].data.astype(np.float64) # time,
# lat, lon
# swap axes to get it to lon, lat, time
grid_fp = np.swapaxes(grid, 0, -1)
times = np.sort(nc['time'][:])
enddate = datetime.datetime.strptime(nc.getncattr('iedate'), '%Y%m%d')
gtime = []
for t in times:
gtime.append(enddate + datetime.timedelta(seconds=int(t)))
valid_file = True
grid_fp *= 1e12 # This is applied in mod_flexpart
return grid_fp, len(gtime), gtime, valid_file
[docs]
def read_grid_stilt(self, path_file, release_date, fp_header):
"""Read a STILT footprint CSV file.
Reads the STILT footprint format (one CSV per observation) and
reconstructs a dense ``(nlon, nlat, nhours)`` footprint array matched
to the CIF domain grid.
Args:
self: Lagrangian model plugin instance (carries ``backward_trajdays``
and ``domain``).
path_file (str): path to the STILT footprint CSV file.
release_date (pd.Timestamp): observation release date.
fp_header: FLEXPART header (for grid metadata; only ``numpoint``
used).
Returns:
tuple: ``(grid_fp, gtime, ngrid, valid_file)``.
"""
nhours = int(pd.to_timedelta(
self.backward_trajdays).total_seconds() / 3600)
grid_fp = np.zeros((self.domain.nlon, self.domain.nlat, nhours))
gtime = release_date.to_pydatetime() \
- np.arange(nhours) * datetime.timedelta(hours=1)
valid_file = False
# Try several times in case there is a hick-up in data
trials = 0
while not valid_file and trials < 3:
try:
with Dataset(path_file, "r") as f:
if release_date.strftime("ftp_%Y%m%d%H_lon") not in f.variables:
info(f"No footprint is available for date {release_date}")
return grid_fp, len(gtime), gtime, valid_file
lon = f.variables[release_date.strftime("ftp_%Y%m%d%H_lon")][:]
lat = f.variables[release_date.strftime("ftp_%Y%m%d%H_lat")][:]
val = f.variables[release_date.strftime("ftp_%Y%m%d%H_val")][:]
ind = f.variables[release_date.strftime(
"ftp_%Y%m%d%H_indptr")][:]
valid_file = True
except RuntimeError:
trials += 1
nb_indexes = np.diff(ind)
hours2process = np.where(nb_indexes != 0)[0]
domain = self.domain
xmin = domain.xmin
xmax = domain.xmax
nlon = domain.nlon
dx = (xmax - xmin) / nlon
ymin = domain.ymin
ymax = domain.ymax
nlat = domain.nlat
dy = (ymax - ymin) / nlat
ind_lon = ((lon - xmin) / dx).astype(int)
ind_lat = ((lat - ymin) / dy).astype(int)
ind_hour = np.zeros(len(lon))
# np.append(0, np.cumsum(nb_indexes))
ind_hour[
np.minimum(len(lon) - 1,
np.append(0, np.cumsum(nb_indexes[hours2process])[:-1]))
] = np.arange(nhours)[hours2process]
ind_hour = np.maximum.accumulate(ind_hour).astype(int) + 1
np.add.at(grid_fp,
(ind_lon, ind_lat, ind_hour),
val)
return grid_fp, len(gtime), gtime, valid_file
[docs]
def read_flexpart_gridinit(subdir, filename, fp_header,
scaleconc=1, **kwargs):
"""Read a FLEXPART initial-condition sensitivity binary file.
Parses ``grid_initial_*_001`` Fortran-binary files and reconstructs the
sparse initial-condition sensitivity as a dense ``(nlon, nlat, nlev)``
array.
Args:
subdir (str): directory containing the file.
filename (str): file name (without directory).
fp_header: FLEXPART domain header object.
scaleconc (float): unit scaling factor.
**kwargs: unused.
Returns:
tuple: ``(grid_fp, gtime, ngrid, valid_file)`` where *grid_fp* is
the dense sensitivity array, *gtime* is the sensitivity date list,
*ngrid* is the time dimension, and *valid_file* is ``True`` on
success.
"""
path_file = os.path.join(subdir, filename)
with FortranFile(path_file, 'r') as f:
yyyymmdd = f.read_ints('i4')[0].astype(str)
hhmmss = f"{f.read_ints('i4')[0]:06d}"
sp_count_i = f.read_ints('i4')[0]
sparse_i = f.read_ints('i4')
sp_count_r = f.read_ints('i4')[0]
sparse_r = f.read_reals('f4')
nx = fp_header.numx
ny = fp_header.numy
nz = (fp_header.outheight != 0.).sum()
grid_init = np.zeros((nx, ny, nz))
if len(sparse_r) > 0:
sign = np.sign(sparse_r)
cum_counts = np.unique(np.cumsum(sp_count_r) % sign.size)
signchange = ((np.roll(sign, 1) - sign) != 0)
signchange[cum_counts] = True
inds_out = np.zeros((len(sparse_r)), dtype=int)
inds_out[signchange] = sparse_i
mask = inds_out == 0
idx = np.where(~mask, np.arange(mask.size), 0)
np.maximum.accumulate(idx, out=idx)
inds_out = inds_out[idx] + np.arange(mask.size) \
- idx - fp_header.numx * fp_header.numy
jz, jy, jx = np.unravel_index(inds_out,
(nz, fp_header.numy, fp_header.numx))
grid_init[jx, jy, jz] = np.abs(sparse_r) * scaleconc
grid_init = np.roll(grid_init, fp_header.xshift, axis=0)
return grid_init
[docs]
def get_spec(subdir, **kwargs):
""" Get species name from simulation output
"""
with open(os.path.join(subdir, 'header_txt'), 'r') as f:
txt = f.readlines()
spec = txt[20].split()[1]
info("species", spec)
return spec
