MetPy interpolate_to_grid 函数基于数据域返回意外的 nan
MetPy interpolate_to_grid function returning unexpected nan based on data domain
我一直在使用 MetPy Mondays #154 中的“Gridding METAR Observations”示例代码一段时间，没有任何问题。直到最近，我都没有限制地传递了整个数据集（除了南极附近的站点，它们破坏了兰伯特共形变换。）
最近，我试图将我处理的 METAR 数据的域限制在北美。此时，MetPy 的 interpolate_to_grid 函数似乎在返回 nan 之前没有返回。由于我感兴趣的区域远离数据集的边界，因此我预计不会对从插值数据导出的轮廓产生影响；相反，这会产生深远的影响（请参见下面的示例。）我尝试使用 SciPy 的 interp2d 对缺失数据区域 (nan) 进行插值 函数，但是有太多 nan 无法用那个“创可贴步骤”克服。
问题：这是 interpolate_to_grid 的预期行为，还是我使用不当？我总是可以继续使用整个数据集，但这确实会减慢速度。感谢您帮助理解这一点。
在下面的示例中，我使用来自 https://tgftp.nws.noaa.gov/data/observations/metar/cycles/ 的 00Z.TXT 文件，但我从其他来源的 METAR 数据中看到这一点。
import cartopy.crs as ccrs
import cartopy.feature as cfeature
import matplotlib.pyplot as plt
import numpy as np
from metpy.io import parse_metar_file
from metpy.interpolate import interpolate_to_grid, remove_nan_observations
%matplotlib inline
# If we are getting data from the filesystem:
month=10
fp = '00Z.TXT'
df0 = parse_metar_file(fp, month=month)
# To avoid the Lambert Conformal transformation from blowing up
# at the South Pole if reports from Antarctica are present:
q = df0.loc[df0['latitude'].values>=-30]
df = q
# Set up the map projection
mapcrs = ccrs.LambertConformal(central_longitude=-100, central_latitude=35,standard_parallels=(30,60))
datacrs= ccrs.PlateCarree()
# 1) Remove NaN
df1=df.dropna(subset=['latitude','longitude','air_temperature'])
lon1=df1['longitude'].values
lat1=df1['latitude'].values
xp1 , yp1 , _ = mapcrs.transform_points(datacrs, lon1, lat1).T
# Interpolate observation data onto grid.
xm1, ym1, tmp = remove_nan_observations(xp1, yp1, df1['air_temperature'].values)
Tgridx, Tgridy, Temp = interpolate_to_grid(xm1, ym1, tmp, hres = 20000, interp_type='cressman')
fig = plt.figure(figsize=(20,15))
ax = fig.add_subplot(1,1,1,projection=mapcrs)
ax.set_extent([-105, -95, 32, 40],datacrs)
ax.add_feature(cfeature.COASTLINE.with_scale('50m'))
ax.add_feature(cfeature.STATES.with_scale('50m'))
c = ax.contour(Tgridx, Tgridy, Temp ,levels=50)
输出是我期望看到的：
当我采用相同的数据集并限制其域时，我们得到不连续的轮廓：
# Limit to ~ North America
q = df0.loc[(df0['latitude'].values>=20) & (df0['latitude'].values<=70) &
(df0['longitude'].values>=-150) & (df0['longitude'].values<=-60)]
df = q
# 1) Remove NaN
df2=df.dropna(subset=['latitude','longitude','air_temperature'])
lon2=df2['longitude'].values
lat2=df2['latitude'].values
xp2 , yp2 , _ = mapcrs.transform_points(datacrs, lon2, lat2).T
# Interpolate observation data onto grid.
xm2, ym2, tmp2 = remove_nan_observations(xp2, yp2, df2['air_temperature'].values)
Tgridx2, Tgridy2, Temp2 = interpolate_to_grid(xm2, ym2, tmp2, hres = 20000, interp_type='cressman')
fig2 = plt.figure(figsize=(20,15))
ax2 = fig2.add_subplot(1,1,1,projection=mapcrs)
ax2.set_extent([-105, -95, 32, 40],datacrs)
ax2.add_feature(cfeature.COASTLINE.with_scale('50m'))
ax2.add_feature(cfeature.STATES.with_scale('50m'))
c2 = ax2.contour(Tgridx2, Tgridy2, Temp2 ,levels=50)
我确实确认站点位置实际上在北美（未显示）。然后我在插值数据中检查了 nan 的位置，发现等高线在区域中断充满了nan。作为最终图，我绘制了 nan（蓝色）的位置、站点位置（绿色）以及断开的等高线。
xn , yn = np.where(np.isnan(Temp2))
fig4 = plt.figure(figsize=(20,15))
ax4 = fig4.add_subplot(1,1,1,projection=mapcrs)
ax4.set_extent([-105, -95, 32, 40],datacrs)
ax4.add_feature(cfeature.COASTLINE.with_scale('50m'))
ax4.add_feature(cfeature.STATES.with_scale('50m'))
c4 = ax4.contour(Tgridx2, Tgridy2, Temp2 ,levels=50)
plt.scatter(df2['longitude'],df2['latitude'],transform=datacrs,color='lightgreen')
plt.scatter(Tgridx2[xn,yn], Tgridy2[xn,yn])
plt.show()
这是一个很难弄清楚的问题。发生的事情是，MetPy 的 Cressman（和 Barnes）插值实现使用距离加权平均值中包含的点的最大搜索半径。如果您不指定此最大搜索半径，它将使用站点之间平均最小间距的 5 倍。
通过将您的数据子集化到大约北美，您构建了站点更靠近的数据子集；这导致搜索半径更小（从 ~150km 减少到 ~66km）。这显然会为您的数据集产生次优结果，我认为部分原因是站点有限。我在这里使用 66 公里的搜索半径在结果之上绘制了站点位置：
你可以看到在站与站之间有相当大的差距的地方出现了漏失。这里最好的解决方案是手动将 search_radius 参数指定为 120km 之类的东西，这似乎给出了合理的结果：
Tgridx, Tgridy, Temp = interpolate_to_grid(xp1, yp1, df1['air_temperature'].values,
hres=20000, interp_type='cressman',
search_radius=120000)
虽然 search_radius 的适当值实际上取决于您要分析的特征以及您愿意将值从观察点分散到多远。我会注意到，您还可以在 interpolate_to_grid.
中使用 Cressman 的其他参数调整其中一些参数
我一直在使用 MetPy Mondays #154 中的“Gridding METAR Observations”示例代码一段时间，没有任何问题。直到最近，我都没有限制地传递了整个数据集（除了南极附近的站点，它们破坏了兰伯特共形变换。）
最近，我试图将我处理的 METAR 数据的域限制在北美。此时，MetPy 的 interpolate_to_grid 函数似乎在返回 nan 之前没有返回。由于我感兴趣的区域远离数据集的边界，因此我预计不会对从插值数据导出的轮廓产生影响；相反，这会产生深远的影响（请参见下面的示例。）我尝试使用 SciPy 的 interp2d 对缺失数据区域 (nan) 进行插值 函数，但是有太多 nan 无法用那个“创可贴步骤”克服。
问题：这是 interpolate_to_grid 的预期行为，还是我使用不当？我总是可以继续使用整个数据集，但这确实会减慢速度。感谢您帮助理解这一点。
在下面的示例中，我使用来自 https://tgftp.nws.noaa.gov/data/observations/metar/cycles/ 的 00Z.TXT 文件，但我从其他来源的 METAR 数据中看到这一点。
import cartopy.crs as ccrs
import cartopy.feature as cfeature
import matplotlib.pyplot as plt
import numpy as np
from metpy.io import parse_metar_file
from metpy.interpolate import interpolate_to_grid, remove_nan_observations
%matplotlib inline
# If we are getting data from the filesystem:
month=10
fp = '00Z.TXT'
df0 = parse_metar_file(fp, month=month)
# To avoid the Lambert Conformal transformation from blowing up
# at the South Pole if reports from Antarctica are present:
q = df0.loc[df0['latitude'].values>=-30]
df = q
# Set up the map projection
mapcrs = ccrs.LambertConformal(central_longitude=-100, central_latitude=35,standard_parallels=(30,60))
datacrs= ccrs.PlateCarree()
# 1) Remove NaN
df1=df.dropna(subset=['latitude','longitude','air_temperature'])
lon1=df1['longitude'].values
lat1=df1['latitude'].values
xp1 , yp1 , _ = mapcrs.transform_points(datacrs, lon1, lat1).T
# Interpolate observation data onto grid.
xm1, ym1, tmp = remove_nan_observations(xp1, yp1, df1['air_temperature'].values)
Tgridx, Tgridy, Temp = interpolate_to_grid(xm1, ym1, tmp, hres = 20000, interp_type='cressman')
fig = plt.figure(figsize=(20,15))
ax = fig.add_subplot(1,1,1,projection=mapcrs)
ax.set_extent([-105, -95, 32, 40],datacrs)
ax.add_feature(cfeature.COASTLINE.with_scale('50m'))
ax.add_feature(cfeature.STATES.with_scale('50m'))
c = ax.contour(Tgridx, Tgridy, Temp ,levels=50)
输出是我期望看到的：
当我采用相同的数据集并限制其域时，我们得到不连续的轮廓：
# Limit to ~ North America
q = df0.loc[(df0['latitude'].values>=20) & (df0['latitude'].values<=70) &
(df0['longitude'].values>=-150) & (df0['longitude'].values<=-60)]
df = q
# 1) Remove NaN
df2=df.dropna(subset=['latitude','longitude','air_temperature'])
lon2=df2['longitude'].values
lat2=df2['latitude'].values
xp2 , yp2 , _ = mapcrs.transform_points(datacrs, lon2, lat2).T
# Interpolate observation data onto grid.
xm2, ym2, tmp2 = remove_nan_observations(xp2, yp2, df2['air_temperature'].values)
Tgridx2, Tgridy2, Temp2 = interpolate_to_grid(xm2, ym2, tmp2, hres = 20000, interp_type='cressman')
fig2 = plt.figure(figsize=(20,15))
ax2 = fig2.add_subplot(1,1,1,projection=mapcrs)
ax2.set_extent([-105, -95, 32, 40],datacrs)
ax2.add_feature(cfeature.COASTLINE.with_scale('50m'))
ax2.add_feature(cfeature.STATES.with_scale('50m'))
c2 = ax2.contour(Tgridx2, Tgridy2, Temp2 ,levels=50)
我确实确认站点位置实际上在北美（未显示）。然后我在插值数据中检查了 nan 的位置，发现等高线在区域中断充满了nan。作为最终图，我绘制了 nan（蓝色）的位置、站点位置（绿色）以及断开的等高线。
xn , yn = np.where(np.isnan(Temp2))
fig4 = plt.figure(figsize=(20,15))
ax4 = fig4.add_subplot(1,1,1,projection=mapcrs)
ax4.set_extent([-105, -95, 32, 40],datacrs)
ax4.add_feature(cfeature.COASTLINE.with_scale('50m'))
ax4.add_feature(cfeature.STATES.with_scale('50m'))
c4 = ax4.contour(Tgridx2, Tgridy2, Temp2 ,levels=50)
plt.scatter(df2['longitude'],df2['latitude'],transform=datacrs,color='lightgreen')
plt.scatter(Tgridx2[xn,yn], Tgridy2[xn,yn])
plt.show()
这是一个很难弄清楚的问题。发生的事情是，MetPy 的 Cressman（和 Barnes）插值实现使用距离加权平均值中包含的点的最大搜索半径。如果您不指定此最大搜索半径，它将使用站点之间平均最小间距的 5 倍。
通过将您的数据子集化到大约北美，您构建了站点更靠近的数据子集；这导致搜索半径更小（从 ~150km 减少到 ~66km）。这显然会为您的数据集产生次优结果，我认为部分原因是站点有限。我在这里使用 66 公里的搜索半径在结果之上绘制了站点位置：
你可以看到在站与站之间有相当大的差距的地方出现了漏失。这里最好的解决方案是手动将 search_radius 参数指定为 120km 之类的东西，这似乎给出了合理的结果：
Tgridx, Tgridy, Temp = interpolate_to_grid(xp1, yp1, df1['air_temperature'].values,
hres=20000, interp_type='cressman',
search_radius=120000)
虽然 search_radius 的适当值实际上取决于您要分析的特征以及您愿意将值从观察点分散到多远。我会注意到，您还可以在 interpolate_to_grid.