import numpy as np
from osgeo import gdal, ogr, osr
from logging import debug
from .....utils.path import init_dir
import matplotlib.pyplot as plt
import geopandas as gpd
import pandas as pd
from shapely.geometry import Polygon
from .....utils.check.errclass import CifError
new_pyproj = True
try:
from pyproj import Proj, Transformer
except ImportError:
from pyproj import Proj, transform
new_pyproj = False
[docs]
def make_hcoord(domain):
"""Build the CHIMERE horizontal coordinate arrays and write HCOORD files.
Supports three grid types controlled by ``domain.type``:
* ``'deg'`` — regular lat-lon grid defined by ``xmin``, ``xmax``,
``ymin``, ``ymax``, ``dx``, ``dy`` (in degrees). An optional buffer
region of coarser resolution can be added with ``add_buffer``.
* ``'km'`` — rotated conformal grid centred on (``xcenter``,
``ycenter``) with cell size ``dx`` × ``dy`` in kilometres, projected
from an oblique Mercator CRS to GPS coordinates via pyproj.
* ``'precomputed'`` — reads pre-existing ``COORD_{domid}`` and
``COORDcorner_{domid}`` text files from ``coord_precomputed_dir``.
After building the coordinate arrays the function:
1. Computes cell areas with :meth:`domain.calc_areas`.
2. Saves a PNG map of cell areas and grid mesh for visual inspection.
3. Writes ``COORD_{domid}`` (centres) and ``COORDcorner_{domid}``
(corners + areas) to ``{workdir}/domain/HCOORD/``.
Sets on *domain*: ``nlon``, ``nlat``, ``zlon``, ``zlat``, ``zlonc``,
``zlatc``.
Args:
domain: CHIMERE domain plugin instance with all geometry parameters
set.
Raises:
Exception: if the domain dimensions are not consistent with ``dx``/
``dy`` (``'deg'`` type) or if ``domain.type`` is unrecognised.
"""
# Different initialization depending on type of domain
domain_type = domain.type
dx = domain.dx
dy = domain.dy
# Regular in degree
if domain_type == "deg":
xmin = domain.xmin
xmax = domain.xmax
ymin = domain.ymin
ymax = domain.ymax
nzo = (xmax - xmin) / dx
nme = (ymax - ymin) / dy
if not np.isclose(nzo, round(nzo)) or not np.isclose(nme, round(nme)):
raise CifError(f"The domain was not well defined. Please specify a size compatible with the resolution. \nCurrent configuration: \n - xmin: {xmin}\n - xmax: {xmax}\n - ymin: {ymin}\n - ymax: {ymax}\n - dx: {dx}\n - dy: {dy}\n"
)
lonc = np.linspace(xmin, xmax, int(nzo + 1))
latc = np.linspace(ymin, ymax, int(nme + 1))
zlatc, zlonc = np.meshgrid(latc, lonc, indexing="ij")
zlon = zlonc[:-1, :-1] + dx / 2
zlat = zlatc[:-1, :-1] + dy / 2
add_buffer = domain.add_buffer
x_buffer_left = domain.x_buffer_left
x_buffer_right = domain.x_buffer_right
y_buffer_up = domain.y_buffer_up
y_buffer_down = domain.y_buffer_down
dx_buffer = domain.dx_buffer
dy_buffer = domain.dy_buffer
if add_buffer:
debug('buffer creation')
zlon, zlat, zlonc, zlatc, nzo, nme = add_buffer_region(xmax, xmin, ymax, ymin, nzo, nme, dx, dy, zlonc, zlatc,
x_buffer_left, x_buffer_right, y_buffer_down, y_buffer_up,
dx_buffer, dy_buffer)
domain.nlon = int(nzo)
domain.nlat = int(nme)
elif domain_type == "km":
xcenter = domain.xcenter
ycenter = domain.ycenter
nlat = domain.nlat
nlon = domain.nlon
stretching = domain.stretching
dx *= 1000
dy *= 1000
inProj = \
f'+proj=omerc +k_0=1.0 +lat_0={ycenter} +alpha=-55 +gamma=0 +lonc={xcenter} '\
'+ellps=WGS84 +x_0=0 +y_0=0 +no_u_off +no_v_off'
# Reference GPS SpatialReference
if new_pyproj:
outProj = 'epsg:4326'
else:
outProj = Proj(init='epsg:4326')
# Meshgrid generation and conversion
zlatc = dy * np.linspace(- nlat / 2., nlat / 2., nlat + 1)
zlat = dy * np.linspace(- (nlat - 1) / 2., (nlat - 1) / 2., nlat)
zlonc = dx * np.linspace(- nlon / 2., nlon / 2., nlon + 1)
zlon = dx * np.linspace(- (nlon - 1) / 2., (nlon - 1) / 2., nlon)
# Stretch longitudes and latitudes
zlonc = np.sign(zlonc) * np.abs(zlonc) \
* (1 + 0.5 * stretching * (np.abs(zlonc) - 1))
zlon = np.sign(zlon) * np.abs(zlon) \
* (1 + 0.5 * stretching * (np.abs(zlon) - 1))
zlatc = np.sign(zlatc) * np.abs(zlatc) \
* (1 + 0.5 * stretching * (np.abs(zlatc) - 1))
zlat = np.sign(zlat) * np.abs(zlat) \
* (1 + 0.5 * stretching * (np.abs(zlat) - 1))
zlatc, zlonc = np.meshgrid(zlatc, zlonc, indexing="ij")
zlat, zlon = np.meshgrid(zlat, zlon, indexing="ij")
if new_pyproj:
t = Transformer.from_crs(inProj, outProj, always_xy=True)
zlonc, zlatc = t.transform(zlonc, zlatc)
zlon, zlat = t.transform(zlon, zlat)
else:
zlonc, zlatc = transform(inProj, outProj, zlonc, zlatc)
zlon, zlat = transform(inProj, outProj, zlon, zlat)
elif domain_type == "precomputed":
centers = np.genfromtxt(
f"{domain.coord_precomputed_dir}/COORD_{domain.domid}")
corners = np.genfromtxt(
f"{domain.coord_precomputed_dir}/COORDcorner_{domain.domid}")
zlon = centers[:, 0].reshape((domain.nlat, domain.nlon))
zlat = centers[:, 1].reshape((domain.nlat, domain.nlon))
zlonc = corners[:, 0].reshape((domain.nlat + 1, domain.nlon + 1))
zlatc = corners[:, 1].reshape((domain.nlat + 1, domain.nlon + 1))
else:
raise CifError(
f"Can't recognized the domain type {domain_type}")
# Fill values in the domain object
domain.zlon = zlon
domain.zlonc = zlonc
domain.zlat = zlat
domain.zlatc = zlatc
# Get areas
domain.calc_areas()
areas = domain.areas
areas_corner = np.concatenate([areas, areas[:, [-1]]], axis=1)
areas_corner = np.concatenate([areas_corner, areas_corner[[-1]]], axis=0)
# Dump HCOORD file
init_dir(f"{domain.workdir}/domain/HCOORD/")
# Plot intermediate map in case
geom = create_2Dpolygons(domain.zlonc, domain.zlatc, domain.zlat.shape)
data = pd.DataFrame(dict(
lat=domain.zlat.flat,
lon=domain.zlon.flat,
area_m2=areas.flat))
gdf_dom = gpd.GeoDataFrame(data, geometry=geom, crs="EPSG:4326")
plt.figure(figsize=(21, 11))
gdf_dom.plot(column='area_m2', legend=True)
plt.savefig(f"{domain.workdir}/domain/HCOORD/map_areas.png")
plt.close()
plt.figure(figsize=(21, 11))
plt.plot(zlonc, zlatc) # use plot, not scatter
plt.plot(np.transpose(zlonc), np.transpose(zlatc)) # add this here
plt.savefig(f"{domain.workdir}/domain/HCOORD/map_mesh.png")
plt.close()
np.savetxt(
f"{domain.workdir}/domain/HCOORD/COORD_{domain.domid}",
np.concatenate([
zlon.reshape((-1, 1)),
zlat.reshape((-1, 1)),
areas.reshape((-1, 1))],
axis=1),
fmt="%10.5f"
)
np.savetxt(
f"{domain.workdir}/domain/HCOORD/COORDcorner_{domain.domid}",
np.concatenate([
zlonc.reshape((-1, 1)),
zlatc.reshape((-1, 1)),
areas_corner.reshape((-1, 1))],
axis=1),
fmt="%10.5f"
)
[docs]
def add_buffer_region(xmax, xmin, ymax, ymin, nzo, nme, dx, dy, zlonc, zlatc,
x_buffer_left, x_buffer_right, y_buffer_down, y_buffer_up,
dx_buffer, dy_buffer):
"""Extend a lat-lon domain with a coarser-resolution buffer frame.
Surrounds the inner high-resolution grid with buffer strips of width
``dx_buffer`` × ``dy_buffer`` on each side, then recomputes the full
cell-centre and corner arrays for the combined grid.
Args:
xmax (float): east longitude of the inner domain.
xmin (float): west longitude of the inner domain.
ymax (float): north latitude of the inner domain.
ymin (float): south latitude of the inner domain.
nzo (float): number of inner cells in the zonal direction.
nme (float): number of inner cells in the meridional direction.
dx (float): inner cell size in degrees (zonal).
dy (float): inner cell size in degrees (meridional).
zlonc (np.ndarray): corner longitudes of the inner domain.
zlatc (np.ndarray): corner latitudes of the inner domain.
x_buffer_left (float): buffer width on the west side (degrees).
x_buffer_right (float): buffer width on the east side (degrees).
y_buffer_down (float): buffer height on the south side (degrees).
y_buffer_up (float): buffer height on the north side (degrees).
dx_buffer (float): buffer cell size in degrees (zonal).
dy_buffer (float): buffer cell size in degrees (meridional).
Returns:
tuple: ``(zlon_buffer, zlat_buffer, zlonc_buffer, zlatc_buffer,
nzo, nme)`` for the combined inner + buffer domain.
"""
# define the number of lines and columns to add in each side
nzo_buffer_left = int(x_buffer_left / dx_buffer)
nzo_buffer_right = int(x_buffer_right / dx_buffer)
nme_buffer_up = int(y_buffer_up / dy_buffer)
nme_buffer_down = int(y_buffer_down / dy_buffer)
nzo = nzo + nzo_buffer_left + nzo_buffer_right
nme = nme + nme_buffer_up + nme_buffer_down
debug("Added lines :")
debug(f"left : {nzo_buffer_left} ; right : {nzo_buffer_right}")
debug(f"up : {nme_buffer_up} ; down : {nme_buffer_down}")
# define the new domain limits
# x_max_buffer = float(np.sign(xmax) * (abs(xmax) + x_buffer_right))
# x_min_buffer = float(np.sign(xmin) * (abs(xmin) + x_buffer_left))
# y_max_buffer = float(np.sign(ymax) * (abs(ymax) + y_buffer_up))
# y_min_buffer = float(np.sign(ymin) * (abs(ymin) + y_buffer_down))
x_max_buffer = float(xmax + x_buffer_right)
x_min_buffer = float(xmin - x_buffer_left)
y_max_buffer = float(ymax + y_buffer_up)
y_min_buffer = float(ymin - y_buffer_down)
debug("Domain limits :")
debug("without buffer :")
debug(f"xmin : {xmin} ; xmax : {xmax}")
debug(f"ymin : {ymin} ; ymax : {ymax}")
debug("with buffer :")
debug(f"xmin : {x_min_buffer} ; xmax : {x_max_buffer}")
debug(f"ymin : {y_min_buffer} ; ymax : {y_max_buffer}")
# preparing the lat lon domain
lonc_buffer = np.linspace(x_min_buffer, x_max_buffer, int(nzo + 1))
latc_buffer = np.linspace(y_min_buffer, y_max_buffer, int(nme + 1))
zlatc_buffer, zlonc_buffer = np.meshgrid(
latc_buffer, lonc_buffer, indexing="ij")
# Fill interior
zlonc_buffer[nme_buffer_down:-nme_buffer_up,
nzo_buffer_left:-nzo_buffer_right] = zlonc
# Fill left
zlonc_buffer[:, :nzo_buffer_left] = np.tile(
np.arange(x_min_buffer, xmin, dx_buffer), (zlonc_buffer.shape[0], 1))
# Fill right
zlonc_buffer[:, -nzo_buffer_right:] = np.tile(np.arange(
xmax + dx_buffer, x_max_buffer + dx_buffer, dx_buffer), (zlonc_buffer.shape[0], 1))
# Fill bottom
zlonc_buffer[:nme_buffer_down,
:] = zlonc_buffer[nme_buffer_down:2 * nme_buffer_down, :]
# Fill top
zlonc_buffer[-nme_buffer_up:, :] = zlonc_buffer[-2 *
nme_buffer_up:-nme_buffer_up, :]
# update buffer lat
# Fill interior
zlatc_buffer[nme_buffer_down:-nme_buffer_up,
nzo_buffer_left:-nzo_buffer_right] = zlatc
# Fill bottom
zlatc_buffer[:nme_buffer_down, :] = np.tile(
np.arange(y_min_buffer, ymin, dy_buffer), (zlatc_buffer.shape[1], 1)).T
# Fill top
zlatc_buffer[-nme_buffer_up:, :] = np.tile(np.arange(
ymax + dy_buffer, y_max_buffer + dy_buffer, dy_buffer), (zlatc_buffer.shape[1], 1)).T
# Fill left
zlatc_buffer[:, :nzo_buffer_left] = zlatc_buffer[:,
nzo_buffer_left:2 * nzo_buffer_left]
# Fill right
zlatc_buffer[:, -nzo_buffer_right:] = zlatc_buffer[:, -
2 * nzo_buffer_right:-nzo_buffer_right]
zlon_buffer = zlonc_buffer.copy()
zlon_buffer[nme_buffer_down:-nme_buffer_up, nzo_buffer_left:-nzo_buffer_right] = (
zlonc_buffer[nme_buffer_down:-nme_buffer_up, nzo_buffer_left:-nzo_buffer_right] + dx / 2)
zlon_buffer[:, :nzo_buffer_left] = (
zlonc_buffer[:, :nzo_buffer_left] + dx_buffer / 2)
zlon_buffer[:, -nzo_buffer_right:] = (
zlonc_buffer[:, -nzo_buffer_right:] - dx_buffer / 2)
zlon_buffer[:nme_buffer_down, :] = (
zlon_buffer[nme_buffer_down:2 * nme_buffer_down, :])
zlon_buffer[-nme_buffer_up:,
:] = (zlon_buffer[-2 * nme_buffer_up:-nme_buffer_up, :])
zlat_buffer = zlatc_buffer.copy()
zlat_buffer[nme_buffer_down:-nme_buffer_up, nzo_buffer_left:-nzo_buffer_right] = (
zlatc_buffer[nme_buffer_down:-nme_buffer_up, nzo_buffer_left:-nzo_buffer_right] + dy / 2)
zlat_buffer[:nme_buffer_down, :] = (
zlatc_buffer[:nme_buffer_down, :] + dy_buffer / 2)
zlat_buffer[-nme_buffer_up:,
:] = (zlatc_buffer[-nme_buffer_up:, :] - dy_buffer / 2)
zlat_buffer[:, :nzo_buffer_left] = (
zlat_buffer[:, nzo_buffer_left:2 * nzo_buffer_left])
zlat_buffer[:, -nzo_buffer_right:] = (
zlat_buffer[:, -2 * nzo_buffer_right:-nzo_buffer_right])
# supress the additional row and column
zlon_buffer = np.delete(zlon_buffer, -nzo_buffer_right - 1, axis=1)
zlon_buffer = np.delete(zlon_buffer, -nme_buffer_up - 1, axis=0)
zlat_buffer = np.delete(zlat_buffer, -nme_buffer_up - 1, axis=0)
zlat_buffer = np.delete(zlat_buffer, -nzo_buffer_right - 1, axis=1)
return zlon_buffer, zlat_buffer, zlonc_buffer, zlatc_buffer, nzo, nme
[docs]
def create_2Dpolygons(x_corner, y_corner, dim):
"""Build a list of Shapely Polygon objects from corner-coordinate arrays.
For each cell ``(i, j)`` in the grid, constructs a quadrilateral polygon
from the four surrounding corner points. Used to create a GeoDataFrame
for diagnostic area maps.
Args:
x_corner (np.ndarray): 2-D array of corner longitudes,
shape ``(nlat+1, nlon+1)``.
y_corner (np.ndarray): 2-D array of corner latitudes,
shape ``(nlat+1, nlon+1)``.
dim (tuple[int, int]): ``(nlat, nlon)`` — number of cells
(not corners) in each direction.
Returns:
list[shapely.geometry.Polygon]: one polygon per grid cell, in
row-major order.
"""
pixel_polygon_list = []
for index in np.ndindex(dim):
i = index[0]
j = index[1]
p1 = [x_corner[i, j], y_corner[i, j]]
p2 = [x_corner[i, j + 1], y_corner[i, j + 1]]
p3 = [x_corner[i + 1, j + 1], y_corner[i + 1, j + 1]]
p4 = [x_corner[i + 1, j], y_corner[i + 1, j]]
pixel_polygon_list.append(Polygon([p1, p2, p3, p4]))
return pixel_polygon_list
