我有两个数据集:一个是建模的,一个是观测的。我需要使观测数据集达到与模型相同的分辨率。
目前,模型数据有100个经度和100个纬度点,因此每个数据点都是1.8度乘3.6度。
我尝试了以下操作,但数据点不太匹配,所以我无法连接数据集。
import xarray as xa
import numpy as np
import cmocean.cm as cm
import matplotlib.patches as mpatches
from matplotlib.patches import Patch
from matplotlib.lines import Line2D
import matplotlib.pyplot as plt
import seaborn as sns
import pandas as pd
import cartopy.crs as ccrs
import cartopy.feature as cfeat
# =============================================================================
# Step 1: Get data
# =============================================================================
data_model = xa.open_dataset("PD_tavg_out_diss_5.nc",decode_times = False)
sal_obs_data = xa.open_dataset("sal_obs_NS_all.nc",decode_times = False)
temp_obs_data = xa.open_dataset("temp_obs_NS_all.nc",decode_times = False)
data_obs = xa.merge([sal_obs_data,temp_obs_data])
data_obs = data_obs.rename(lon = "longitude",
lat = "latitude")
data_extract = data_model[["O_cons_temp","O_abs_sal"]]
long1 = data_extract.longitude.values
long1[long1>180]-=360
data_extract["longitude"] = long1
data_sorted = data_extract.sortby("longitude")
long2 = data_obs.longitude.values
long2[long2>180]-=360
data_obs["longitude"] = long2
data_sorted_obs = data_obs.sortby("longitude")
long_max = 1.8
long_min = -1.8
lat_max = 60.3
lat_min = 54.9
dep_max = 100
dep_min = 0
tim_max = 35422.0
tim_min = 35421.0
def extract_shelf_sea(long_max, long_min,
lat_max, lat_min,
dep_max, dep_min,
tim_max, tim_min):
# =============================================================================
# Step 3: Extract data
# =============================================================================
extract_model_data = data_sorted.sel(longitude = slice(long_min,long_max),
latitude = slice(lat_min,lat_max),
depth = slice(dep_min,dep_max),
time = slice(tim_min,tim_max))
extract_obs_data = data_obs.sel(time = data_obs.time,
longitude = slice(long_min,long_max),
latitude = slice(lat_min,lat_max))
new_lon = np.linspace(extract_obs_data.longitude[0],extract_obs_data.longitude[7],extract_model_data.sizes['longitude'])
new_lat = np.linspace(extract_obs_data.latitude[0],extract_obs_data.latitude[10],extract_model_data.sizes['latitude'])
obs_interpolated = extract_obs_data.interp(latitude = new_lat, longitude = new_lon)
extract_obs_depth = obs_interpolated.sel(depth = extract_model_data.depth, method="nearest")
观察输出:
extract_obs_depth
Out[178]:
<xarray.Dataset>
Dimensions: (depth: 2, latitude: 4, longitude: 2, time: 12)
Coordinates:
* depth (depth) float64 15.07 82.92
* time (time) float32 480.5 481.5 482.5 483.5 ... 489.5 490.5 491.5
* latitude (latitude) float64 54.25 56.42 58.58 60.75
* longitude (longitude) float64 -1.75 1.75
Data variables:
salt (time, depth, latitude, longitude) float64 nan 34.56 ... 35.38
temp (time, depth, latitude, longitude) float64 nan 6.586 ... 8.907
型号输出:
extract_model_data
Out[179]:
<xarray.Dataset>
Dimensions: (depth: 2, latitude: 4, longitude: 2, time: 12)
Coordinates:
* time (time) float64 3.542e+04 3.542e+04 ... 3.542e+04 3.542e+04
* longitude (longitude) float64 -1.8 1.8
* latitude (latitude) float64 54.9 56.7 58.5 60.3
* depth (depth) float64 17.5 82.5
Data variables:
O_cons_temp (time, depth, latitude, longitude) float64 5.615 ... 7.437
O_abs_sal (time, depth, latitude, longitude) float64 33.28 ... 35.21
任何关于如何使纬度和经度完美匹配的建议都将不胜感激。提前谢谢。
顺便说一句,不幸的是,我无法上传数据集,因为它们太大了。
将观测数据集插值到建模数据集的坐标怎么样?
特别是,您可以使用Scipy模块中的RectSphereBivariateSpline
,它直接在球面坐标中插值结构化数据(在规则网格中定义的数据,就像您的情况一样((我在这里假设您的其他维度-时间和深度-数据集之间匹配;如果不匹配,您将不得不执行多次插值(。
看看我上面链接的文档中的例子,如果你有任何疑问,请随时询问;我自己也用过一两次这个函数,如果我没有记错的话,文档中的一个坐标参数是错误的(我认为这是v
参数,它不是像文档所说的那样在[-pi, pi]
区间中定义的,而是实际上必须在[0, 2pi]
中——不过要谨慎对待(。