resample.py 3.98 KiB
import numpy as np
from scipy.interpolate import RegularGridInterpolator, griddata
from cartopy.crs import LambertAzimuthalEqualArea
import cartopy.crs as ccrs
from netCDF4 import Dataset
import xarray as xr
# resample methods:
linear = 'linear'
cubic = 'cubic'
nearest = 'nearest'
def fill_missing(fld_a, fld_b, mask):
min_val = min(np.nanmin(fld_a), np.nanmin(fld_b))
max_val = max(np.nanmax(fld_a), np.nanmax(fld_b))
random_values_a = np.random.uniform(min_val, max_val, size=fld_a.shape)
random_values_b = np.random.uniform(min_val, max_val, size=fld_a.shape)
fld_a[mask] = random_values_a[mask]
fld_b[mask] = random_values_b[mask]
def get_projection(cartopy_map_name, cen_lat, cen_lon):
if cartopy_map_name == "LambertAzimuthalEqualArea":
projection = ccrs.LambertAzimuthalEqualArea(
central_longitude=cen_lon, central_latitude=cen_lat
)
elif cartopy_map_name == "AlbersEqualArea":
projection = ccrs.AlbersEqualArea(
central_longitude=cen_lon, central_latitude=cen_lat
)
elif cartopy_map_name == "Sinusoidal":
projection = ccrs.Sinusoidal(central_longitude=cen_lon)
else:
raise ValueError("Projection: " + cartopy_map_name + " is not supported")
return projection
def resample_reg_grid(scalar_field, y, x, y_s, x_s, method='linear'):
intrp = RegularGridInterpolator((y, x), scalar_field, method=method, bounds_error=False)
xg, yg = np.meshgrid(x_s, y_s, indexing='xy')
yg, xg = yg.flatten(), xg.flatten()
pts = np.array([yg, xg])
t_pts = np.transpose(pts)
return np.reshape(intrp(t_pts), (y_s.shape[0], x_s.shape[0]))
def resample(scalar_field, y_d, x_d, y_t, x_t, method='linear'):
# 2D target locations shape
t_shape = y_t.shape
# reproject scalar fields
fld_repro = griddata((y_d.flatten(), x_d.flatten()), scalar_field.flatten(),(y_t.flatten(), x_t.flatten()), method=method)
fld_repro = fld_repro.reshape(t_shape)
return fld_repro
def reproject(fld_2d, lat_2d, lon_2d, proj, target_grid=None, grid_spacing=15000, method=linear):
"""
:param fld_2d: the 2D scalar field to reproject
:param lat_2d: 2D latitude of the scalar field domain
:param lon_2d: 2D longitude of the scalar field domain
:param proj: the map projection (Cartopy). Default: LambertAzimuthalEqualArea
:param region_grid: the larger region grid that we pull the target grid from
:param target_grid: the resampling target (y_map, x_map) where y_map and x_map are 2D. If None, the grid is created
automatically. The target grid is always returned.
:param grid_spacing: distance between the target grid points (in meters)
:param method: resampling method: 'linear', 'nearest', 'cubic'
:return: reprojected 2D scalar field, the target grid (will be 2D if rotate=True)
"""
data_xy = proj.transform_points(ccrs.PlateCarree(), lon_2d, lat_2d)[..., :2]
# Generate a regular 2d grid extending the min and max of the xy dimensions with grid_spacing
if target_grid is None:
x_min, y_min = np.amin(data_xy, axis=(0, 1))
x_max, y_max = np.amax(data_xy, axis=(0, 1))
x_map = np.arange(x_min, x_max, grid_spacing)
y_map = np.arange(y_min, y_max, grid_spacing)
x_map_2d, y_map_2d = np.meshgrid(x_map, y_map)
else:
y_map_2d, x_map_2d = target_grid
fld_reproj = resample(fld_2d, data_xy[..., 1], data_xy[..., 0], y_map_2d, x_map_2d, method=method)
return fld_reproj, (y_map_2d, x_map_2d)
def bisect_great_circle(lon_a, lat_a, lon_b, lat_b):
lon_a = np.radians(lon_a)
lat_a = np.radians(lat_a)
lon_b = np.radians(lon_b)
lat_b = np.radians(lat_b)
dlon = lon_b - lon_a
Bx = np.cos(lat_b) * np.cos(dlon)
By = np.cos(lat_b) * np.sin(dlon)
lat_c = np.arctan2(np.sin(lat_a) + np.sin(lat_b), np.sqrt((np.cos(lat_a) + Bx) ** 2 + By ** 2))
lon_c = lon_a + np.arctan2(By, np.cos(lat_a) + Bx)
lon_c = np.degrees(lon_c)
lat_c = np.degrees(lat_c)
return lon_c, lat_c