我有两组卫星数据。对于这两个集合，我都有像素几何（像素每个角的纬度和经度）。我想将一组重新设置为另一组。因此，我的目标是从不规则网格到另一个不规则网格的区域加权重新网格化。我知道xESMF，但不确定这是否是这项工作的最佳工具。也许虹膜区域加权重新网格是合适的？
问问题
161 次
1 回答
1
我过去也遇到过类似的事情。我在 Windows 上，xEMSF 对我来说并不是一个真正的选择。
我已经编写了这个包，并添加了一些计算网格到网格权重的方法： https ://github.com/Deltares/numba_celltree
（你可以点安装它。）
数据结构可以处理完全非结构化的 2D 网格，并期望采用这种格式的数据。请参阅下面的代码。
您将需要进行一些更改：您的坐标很可能不会命名为 x 和 y。您还需要ugrid2d_topology稍微更新函数，因为我在这里假设规则的四边形网格（但在彼此的坐标系中看到它们是不规则的）。
它仍然非常简单，只需确保您有 2D 顶点数组，以及为每个面映射其四个顶点face_node_connectivity的形状数组。(n_cell, 4)有关更多背景信息，请参阅此文档：
https ://ugrid-conventions.github.io/ugrid-conventions/
import numpy as np
import pandas as pd
import pyproj
import xarray as xr
from numba_celltree import CellTree2d
FloatArray = np.ndarray
IntArray = np.ndarray
def _coord(da, dim):
"""
Transform N xarray midpoints into N + 1 vertex edges
"""
delta_dim = "d" + dim # e.g. dx, dy, dz, etc.
# If empty array, return empty
if da[dim].size == 0:
return np.array(())
if delta_dim in da.coords: # equidistant or non-equidistant
dx = da[delta_dim].values
if dx.shape == () or dx.shape == (1,): # scalar -> equidistant
dxs = np.full(da[dim].size, dx)
else: # array -> non-equidistant
dxs = dx
_check_monotonic(dxs, dim)
else: # undefined -> equidistant
if da[dim].size == 1:
raise ValueError(
f"DataArray has size 1 along {dim}, so cellsize must be provided"
" as a coordinate."
)
dxs = np.diff(da[dim].values)
dx = dxs[0]
atolx = abs(1.0e-4 * dx)
if not np.allclose(dxs, dx, atolx):
raise ValueError(
f"DataArray has to be equidistant along {dim}, or cellsizes"
" must be provided as a coordinate."
)
dxs = np.full(da[dim].size, dx)
dxs = np.abs(dxs)
x = da[dim].values
if not da.indexes[dim].is_monotonic_increasing:
x = x[::-1]
dxs = dxs[::-1]
# This assumes the coordinate to be monotonic increasing
x0 = x[0] - 0.5 * dxs[0]
x = np.full(dxs.size + 1, x0)
x[1:] += np.cumsum(dxs)
return x
def _ugrid2d_dataset(
node_x: FloatArray,
node_y: FloatArray,
face_x: FloatArray,
face_y: FloatArray,
face_nodes: IntArray,
) -> xr.Dataset:
ds = xr.Dataset()
ds["mesh2d"] = xr.DataArray(
data=0,
attrs={
"cf_role": "mesh_topology",
"long_name": "Topology data of 2D mesh",
"topology_dimension": 2,
"node_coordinates": "node_x node_y",
"face_node_connectivity": "face_nodes",
"edge_node_connectivity": "edge_nodes",
},
)
ds = ds.assign_coords(
node_x=xr.DataArray(
data=node_x,
dims=["node"],
)
)
ds = ds.assign_coords(
node_y=xr.DataArray(
data=node_y,
dims=["node"],
)
)
ds["face_nodes"] = xr.DataArray(
data=face_nodes,
coords={
"face_x": ("face", face_x),
"face_y": ("face", face_y),
},
dims=["face", "nmax_face"],
attrs={
"cf_role": "face_node_connectivity",
"long_name": "Vertex nodes of mesh faces (counterclockwise)",
"start_index": 0,
"_FillValue": -1,
},
)
ds.attrs = {"Conventions": "CF-1.8 UGRID-1.0"}
return ds
def ugrid2d_topology(data: Union[xr.DataArray, xr.Dataset]) -> xr.Dataset:
"""
Derive the 2D-UGRID quadrilateral mesh topology from a structured DataArray
or Dataset, with (2D-dimensions) "y" and "x".
Parameters
----------
data: Union[xr.DataArray, xr.Dataset]
Structured data from which the "x" and "y" coordinate will be used to
define the UGRID-2D topology.
Returns
-------
ugrid_topology: xr.Dataset
Dataset with the required arrays describing 2D unstructured topology:
node_x, node_y, face_x, face_y, face_nodes (connectivity).
"""
# Transform midpoints into vertices
# These are always returned monotonically increasing
x = data["x"].values
xcoord = _coord(data, "x")
if not data.indexes["x"].is_monotonic_increasing:
xcoord = xcoord[::-1]
y = data["y"].values
ycoord = _coord(data, "y")
if not data.indexes["y"].is_monotonic_increasing:
ycoord = ycoord[::-1]
# Compute all vertices, these are the ugrid nodes
node_y, node_x = (a.ravel() for a in np.meshgrid(ycoord, xcoord, indexing="ij"))
face_y, face_x = (a.ravel() for a in np.meshgrid(y, x, indexing="ij"))
linear_index = np.arange(node_x.size, dtype=np.int32).reshape(
ycoord.size, xcoord.size
)
# Allocate face_node_connectivity
nfaces = (ycoord.size - 1) * (xcoord.size - 1)
face_nodes = np.empty((nfaces, 4))
# Set connectivity in counterclockwise manner
face_nodes[:, 0] = linear_index[:-1, 1:].ravel() # upper right
face_nodes[:, 1] = linear_index[:-1, :-1].ravel() # upper left
face_nodes[:, 2] = linear_index[1:, :-1].ravel() # lower left
face_nodes[:, 3] = linear_index[1:, 1:].ravel() # lower right
# Tie it together
ds = _ugrid2d_dataset(node_x, node_y, face_x, face_y, face_nodes)
return ds
def area_weighted_mean(
da: xr.DataArray,
destination_index: np.ndarray,
source_index: np.ndarray,
weights: np.ndarray,
):
"""
Area weighted mean.
Parameters
----------
da: xr.DataArray
Contains source data.
destination_index: np.ndarray
In which destination the overlap is located.
source_index: np.ndarray
In which source cell the overlap is located.
weights: np.ndarray
Area of each overlap.
Returns
-------
destination_index: np.ndarray
values: np.ndarray
"""
values = da.data.ravel()[source_index]
df = pd.DataFrame(
{"dst": destination_index, "area": weights, "av": weights * values}
)
aggregated = df.groupby("dst").sum("sum", min_count=1)
out = aggregated["av"] / aggregated["area"]
return out.index.values, out.values
class Regridder:
"""
Regridder to reproject and/or regrid rasters. When no ``crs_source`` and
``crs_destination`` are provided, it is assumed that ``source`` and
``destination`` share the same coordinate system.
Note that an area weighted regridding method only makes sense for projected
(Cartesian!) coordinate systems.
Parameters
----------
source: xr.DataArray
Source example. Must have dimensions ("y", "x").
destination: xr.DataArray
Destination example. Must have dimensions ("y", "x").
crs_source: optional, default: None
crs_destination: optional, default: None
"""
def __init__(
self,
source: xr.DataArray,
destination: xr.DataArray,
crs_source=None,
crs_destination=None,
):
src = ugrid2d_topology(source)
dst = ugrid2d_topology(destination)
src_yy = src["node_y"].values
src_xx = src["node_x"].values
if crs_source and crs_destination:
transformer = pyproj.Transformer.from_crs(
crs_from=crs_source, crs_to=crs_destination, always_xy=True
)
src_xx, src_yy = transformer.transform(xx=src_xx, yy=src_yy)
elif crs_source ^ crs_destination:
raise ValueError("Received only one of (crs_source, crs_destination)")
src_vertices = np.column_stack([src_xx, src_yy])
src_faces = src["face_nodes"].values.astype(int)
dst_vertices = np.column_stack((dst["node_x"].values, dst["node_y"].values))
dst_faces = dst["face_nodes"].values
celltree = CellTree2d(src_vertices, src_faces, fill_value=-1)
self.source = source.copy()
self.destination = destination.copy()
(
self.destination_index,
self.source_index,
self.weights,
) = celltree.intersect_faces(
dst_vertices,
dst_faces,
fill_value=-1,
)
def regrid(self, da: xr.DataArray, fill_value=np.nan):
"""
Parameters
----------
da: xr.DataArray
Data to regrid.
fill_value: optional, default: np.nan
Default value of the output grid, e.g. where no overlap occurs.
Returns
-------
regridded: xr.DataArray
Data of da, regridded using an area weighted mean.
"""
src = self.source
if not (np.allclose(da["y"], src["y"]) and np.allclose(da["x"], src["x"])):
raise ValueError("da does not match source")
index, values = area_weighted_mean(
da,
self.destination_index,
self.source_index,
self.weights,
)
data = np.full(self.destination.shape, fill_value)
data.ravel()[index] = values
out = self.destination.copy(data=data)
out.name = da.name
return out
# Example use
da = xr.open_dataarray("gw_abstraction_sum.nc")
like = xr.open_dataarray("example.nc")
regridder = Regridder(
source=da, destination=like, crs_source=4326, crs_destination=3035
)
result = regridder.regrid(da)
result.to_netcdf("area-weighted_sum.nc")
于 2021-10-23T11:22:51.057 回答