import copy
from logging import debug
from pathlib import Path
import numpy as np
import xarray as xr
from ....utils.check.errclass import CifError, CifValueError, IllegalArgumentError
from .utils import find_coord, get_bounds, get_centers
[docs]
def read_vcoord(domain):
"""Read vertical coordinate information from a NetCDF file into the domain"""
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 = (
Path(domain.dir_vcoord, domain.file_vcoord)
if hasattr(domain, 'file_vcoord')
else Path(domain.dir, domain.file)
)
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 CifError(
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 CifValueError(
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 CifValueError(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.vflip = False
elif np.all(np.diff(pressures)) > 0:
# Wrong order flip upside down
sigma_a = np.flip(sigma_a)
sigma_b = np.flip(sigma_b)
domain.vflip = True
else:
raise CifValueError(
"Levels are not in a stricly increasing or decreasing order. "
f"Please check your vertical coordinates 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], heights[:, 1]])
else:
domain.heights = heights
domain.nlev = len(domain.heights) - 1
else:
raise CifValueError(
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}"
)
[docs]
def read_grid(domain, **kwargs):
"""Read horizontal and vertical grid coordinates from a NetCDF file"""
file_path = Path(domain.dir, domain.file)
delta_lat = getattr(domain, 'delta_lat', None)
delta_lon = getattr(domain, 'delta_lon', None)
kw = {"decode_times": False}
if hasattr(domain, "drop_variables"):
kw["drop_variables"] = domain.drop_variables
# Check if drop_variables is needed or need to be uptated by the user
try:
with xr.open_dataset(file_path, **kw) as ds:
pass
# FIXME: This is unsafe xarra.core is not an exposed part of the xarray API.
# everything in xarray.core is subject to change without notice.
except xr.core.variable.MissingDimensionsError as e:
raise CifError(
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 to find out which variable is the problem."
) from e
# Process the reference file
with xr.open_dataset(file_path, **kw) 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.lon_cyclic = (lon[0] % 360.0) == (lon[-1] % 360.0) or (lonc[0] % 360.0) == (
lonc[-1] % 360.0
)
# Unwrap longitudes to avoid big gaps
if domain.lon_cyclic:
lonc = np.unwrap(lonc, period=360)
lon = np.unwrap(lon, period=360)
# Makes mesh
domain.zlon, domain.zlat = np.meshgrid(lon, lat)
domain.zlonc, domain.zlatc = np.meshgrid(lonc, latc)
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
read_vcoord(domain)
