import copy
import os
from logging import debug
import numpy as np
import xarray as xr
from ....utils.check.errclass import IllegalArgumentError
from .utils import get_bounds, get_centers, find_coord
[docs]
def read_grid(domain, **kwargs):
"""Reads a grid from a NetCDF file
Args:
domain (Plugin): dictionary defining the domain
Return:
Grid dictionary with meshgrids for center lon/lat and corner lon/lat
"""
file_path = os.path.join(domain.dir, domain.file)
delta_lat = getattr(domain, 'delta_lat', None)
delta_lon = getattr(domain, 'delta_lon', None)
# Decode times or not for very badly formatted files...
try:
with xr.open_dataset(file_path) as ds:
decode_times = True
except ValueError:
decode_times = False
# Filter MissingDimensionErrors
try:
with xr.open_dataset(
file_path,
decode_times=decode_times,
drop_variables=domain.drop_variables
) as ds:
pass
except xr.core.variable.MissingDimensionsError as e:
raise IOError(
f"Could not open {file_path} with xarray. "
"This happens when one dimension variable has multiple dimensions, "
"with a name common with one of its own dimensions. "
"Please consider reformatting your NetCDF properly. "
"It is possible to use the argument `drop_variables` in your yaml "
"to avoid reading the erroneous variables. \n"
"Please check the error message above to find out which variable is the problem."
) from e
# Process the reference file
with xr.open_dataset(
file_path,
decode_times=decode_times,
drop_variables=domain.drop_variables
) as ds:
# Getting grid coordinate names
lat_coord_name = find_coord(ds, domain.latitude_varname, 'lat')
lon_coord_name = find_coord(ds, domain.longitude_varname, 'lon')
# Getting coordinates values
lat = ds[lat_coord_name].values
lon = ds[lon_coord_name].values
# Flipping coordinates if they are not strictly increasing
flip_lat = not np.all(np.diff(lat) > 0)
flip_lon = not np.all(np.diff(lon) > 0)
lat = np.flip(lat) if flip_lat else lat
lon = np.flip(lon) if flip_lon else lon
# Asserting coordinates are strictly increasing
assert np.all(np.diff(lat) > 0), \
"latitude coordinate is not strictly increasing"
assert np.all(np.diff(lon) > 0), \
"longitude coordinate is not strictly increasing"
# Getting bounds
if not domain.use_corners:
if delta_lat is None:
latc = get_bounds(ds, lat_coord_name, flip_lat,
domain.lat_min, domain.lat_max,
domain.regular_lat)
else:
if "lat_max" in domain.is_default_value or "lat_min" in domain.is_default_value:
raise IllegalArgumentError(
"Using specified `delta_lat` from the configuration "
"while `lat_max` and/or `lat_min` are not explicitly determined. \n"
"Please include them in the Yaml consistently to your domain extension."
)
latc = np.arange(
domain.lat_min, domain.lat_max + delta_lat,
delta_lat)
lat = latc[:-1] + 0.5 * delta_lat
debug(
"Generated fixed latitude corners as followed:\n"
f"\t- lat_min = {domain.lat_min}\n"
f"\t- lat_max = {domain.lat_max}\n"
f"\t- delta_lat = {domain.delta_lat}\n"
f"\t- latc = {latc}\n"
)
if delta_lon is None:
lonc = get_bounds(ds, lon_coord_name, flip_lon,
domain.lon_min, domain.lon_max,
domain.regular_lon)
else:
if "lon_max" in domain.is_default_value or "lon_min" in domain.is_default_value:
raise IllegalArgumentError(
"Using specified `delta_lon` from the configuration "
"while `lon_max` and/or `lon_min` are not explicitly determined. \n"
"Please include them in the Yaml consistently to your domain extension."
)
lonc = np.arange(
domain.lon_min, domain.lon_max + delta_lon,
delta_lon)
else:
latc = copy.deepcopy(lat)
lonc = copy.deepcopy(lon)
if domain.extend_lat:
latc = np.append(latc, 2 * latc[-1] - latc[-2])
if domain.extend_lon:
lonc = np.append(lonc, 2 * lonc[-1] - lonc[-2])
# TODO: Isn't there one excess flip there ??
lat = get_centers(latc, flip_lat)
lon = get_centers(lonc, flip_lon)
if domain.force_centered_lat:
lat = 0.5 * (latc[1:] + latc[:-1])
if domain.force_centered_lon:
lon = 0.5 * (lonc[1:] + lonc[:-1])
# Putting back lat and lon in their original order if necessary
if not domain.sort_lat and flip_lat:
lat = np.flip(lat)
latc = np.flip(latc)
if not domain.sort_lon and flip_lon:
lon = np.flip(lon)
lonc = np.flip(lonc)
domain.nlat = len(lat)
domain.nlon = len(lon)
domain.zlon, domain.zlat = np.meshgrid(lon, lat)
domain.zlonc, domain.zlatc = np.meshgrid(lonc, latc)
domain.lon_cyclic = (
(lon[0] % 360.0) == (lon[-1] % 360.0)
or (lonc[0] % 360.0) == (lonc[-1] % 360.0)
)
debug("Domain latitudes:"
f"\n domain.nlat={domain.nlat}"
f"\n lat={lat}"
f"\n latc={latc}")
debug("Domain longitudes:"
f"\n domain.nlon={domain.nlon}"
f"\n lon={lon}"
f"\n lonc={lonc}"
f"\n domain.lon_cyclic={domain.lon_cyclic}")
# Now deal with vertical coordinates
if not hasattr(domain, 'vertical_coord'):
# Only one level
domain.pressure_unit = 'Pa'
domain.sigma_a_mid = np.array([0])
domain.sigma_b_mid = np.array([1])
domain.nlev = len(domain.sigma_b_mid)
else:
# Getting vertical levels file path
file_vcoord_path = os.path.join(domain.dir_vcoord, domain.file_vcoord) \
if hasattr(domain, 'file_vcoord') else file_path
with xr.open_dataset(file_vcoord_path) as ds:
if domain.vertical_dim_type == "sigma_pressure":
# Reading levels
sigma_a = ds[domain.sigma_a_var_name].values
sigma_b = ds[domain.sigma_b_var_name].values
# Reading units
units = getattr(domain, 'pressure_unit',
ds[domain.sigma_a_var_name].attrs.get('units'))
elif domain.vertical_dim_type == "heights":
heights = ds[domain.vertical_dim_name].values
units = "m"
else:
raise Exception(
f'Type {domain.vertical_dim_type} is not recognized '
'as a valid `vertical_dim_type`. Please check your yaml.'
)
if domain.reverse_vcoord:
if domain.vertical_dim_type == "sigma_pressure":
sigma_a = np.flip(sigma_a)
sigma_b = np.flip(sigma_b)
else:
heights = np.flip(heights)
if units is None:
raise ValueError(
f"no pressure units found in '{domain.sigma_a_var_name}' "
"variable attributes, please use the argument 'pressure_unit' "
"to provide one.")
elif units not in ['Pa', 'hPa', 'm']:
raise ValueError(f"unsupported pressure units {units}")
domain.pressure_unit = units
# Re-order vertical levels if needed
if domain.vertical_dim_type == "sigma_pressure":
psurf = 101300 if domain.pressure_unit == "Pa" else 1013
pressures = psurf * sigma_b + sigma_a
if np.all(np.diff(pressures)) < 0:
# Correct order, do nothing
domain.vshifted = False
elif np.all(np.diff(pressures)) > 0:
# Wrong order shift upside down
sigma_a = sigma_a[::-1]
sigma_b = sigma_b[::-1]
domain.vshifted = True
else:
raise Exception(
"Levels are not in a monotonous order, "
f"please check your file {file_vcoord_path}"
)
# Computing values at level mid points
if domain.vertical_coord == 'mids':
if domain.vertical_dim_type == "sigma_pressure":
domain.sigma_a_mid = sigma_a
domain.sigma_b_mid = sigma_b
domain.nlev = len(domain.sigma_b_mid)
else:
domain.heights_mid = heights
domain.nlev = len(domain.heights_mid)
elif domain.vertical_coord == 'bounds':
if domain.vertical_dim_type == "sigma_pressure":
if sigma_a.ndim == 2:
domain.sigma_a = np.concatenate(
[sigma_a[[0], 0], sigma_a[:, 1]])
domain.sigma_b = np.concatenate(
[sigma_b[[0], 0], sigma_b[:, 1]])
else:
domain.sigma_a = sigma_a
domain.sigma_b = sigma_b
domain.nlev = len(domain.sigma_b) - 1
else:
if heights.ndim == 2:
domain.heights = np.concatenate(
[heights[[0], 0], sheightsigma_a[:, 1]])
else:
domain.heights = heights
domain.nlev = len(domain.heights) - 1
else:
raise ValueError(f"unexpected value '{domain.vertical_coord}' for "
"argument 'vertical_coord'")
if domain.vertical_dim_type == "sigma_pressure":
if hasattr(domain, "sigma_a_mid"):
debug("Domain levels:"
f"\n domain.nlev={domain.nlev}"
f"\n domain.sigma_a_mid={domain.sigma_a_mid}"
f"\n domain.sigma_b_mid={domain.sigma_b_mid}"
f"\n domain.pressure_unit={domain.pressure_unit}")
else:
debug("Domain levels:"
f"\n domain.nlev={domain.nlev}"
f"\n domain.sigma_a={domain.sigma_a}"
f"\n domain.sigma_b={domain.sigma_b}"
f"\n domain.pressure_unit={domain.pressure_unit}")
else:
if hasattr(domain, "heights_mid"):
debug("Domain levels:"
f"\n domain.nlev={domain.nlev}"
f"\n domain.heights_mid={domain.heights_mid}"
f"\n domain.pressure_unit={domain.pressure_unit}")
else:
debug("Domain levels:"
f"\n domain.nlev={domain.nlev}"
f"\n domain.heights={domain.heights}"
f"\n domain.pressure_unit={domain.pressure_unit}")
