Python - 温度插值

Question

我目前正在尝试使用此网站上提供的代码(https://unidata.github.io/MetPy/latest/examples/gridding/Point_Interpolation.html#sphx-glr-examples-gridding-point-interpolation-py(来创建台湾地图，并在Jupyter笔记本上对数据进行线性插值。

我的数据采用以下形式：

17070123, lat, lon, tem
C0A92, 25.27, 121.56, 29.3
C0AD0, 25.26, 121.49, 28.2
C0A94, 25.23, 121.64, 26.2
46691, 25.19, 121.52, 23.4
46690, 25.17, 121.44, 27.3
46693, 25.17, 121.54, 22.5
C0AD1, 25.15, 121.4, 28.5
46694, 25.13, 121.73, 28.6
C0A95, 25.13, 121.92, -999
C0A9B, 25.12, 121.51, 26.8
C0A9C, 25.12, 121.53, 28.3
C0A66, 25.11, 121.79, 27.8
C0A98, 25.11, 121.46, 29.6
C0A68, 25.09, 121.43, -999

并且也以这种形式：

#17070123   lat lon T
C0A92   25.27   121.56  29.3
C0AD0   25.26   121.49  28.2
C0A94   25.23   121.64  26.2
46691   25.19   121.52  23.4
46690   25.17   121.44  27.3
46693   25.17   121.54  22.5
C0AD1   25.15   121.4   28.5
46694   25.13   121.73  28.6
C0A95   25.13   121.92  -999
C0A9B   25.12   121.51  26.8
C0A9C   25.12   121.53  28.3
C0A66   25.11   121.79  27.8
C0A98   25.11   121.46  29.6
C0A68   25.09   121.43  -999

我的代码如下所示：

# In[1]:
import cartopy
import cartopy.crs as ccrs
from matplotlib.colors import BoundaryNorm
import matplotlib.pyplot as plt
import numpy as np

# In[2]:
from metpy.cbook import get_test_data
from metpy.gridding.gridding_functions import (interpolate, 
remove_nan_observations,
remove_repeat_coordinates)

# In[3]:
def basic_map(map_proj):
"""Make our basic default map for plotting"""
fig = plt.figure(figsize=(15, 10))
view = fig.add_axes([0, 0, 1, 1], projection=to_proj)
view.set_extent([120.5, 122.5, 24.5, 25.5])
view.add_feature(cartopy.feature.NaturalEarthFeature(category='cultural',
name='admin_1_states_provinces_lakes',
scale='50m', 
facecolor='none'))
view.add_feature(cartopy.feature.OCEAN)
view.add_feature(cartopy.feature.COASTLINE)
view.add_feature(cartopy.feature.BORDERS, linestyle=':')
return view

# In[4]:
def station_test_data(variable_names, proj_from=None, proj_to=None):
f = ('temp.txt')
all_data = np.loadtxt(f, skiprows=0, delimiter='t',
usecols=(0, 1, 2, 3),
dtype=np.dtype([('stid', '5S'), ('lat', 'f'), ('lon', 
'f'), ('T', 'f')]))
all_stids = [s.decode('ascii') for s in all_data['stid']]
data = np.concatenate([all_data[all_stids.index(site)].reshape(1, ) for 
site in all_stids])
value = data[variable_names]
lon = data['lon']
lat = data['lat']
if proj_from is not None and proj_to is not None:
try:
proj_points = proj_to.transform_points(proj_from, lon, lat)
return proj_points[:, 0], proj_points[:, 1], value
except Exception as e:
print(e)
return None
return lon, lat, value

# In[5]:
from_proj = ccrs.Geodetic()
to_proj = ccrs.AlbersEqualArea(central_longitude=120.0000, 
central_latitude=25.0000)

# In[6]:
levels = list(range(20, 30, 1))
cmap = plt.get_cmap('magma')
norm = BoundaryNorm(levels, ncolors=cmap.N, clip=True)

# In[7]:
x, y, temp = station_test_data('T', from_proj, to_proj)

# In[8]:
x, y, temp = remove_nan_observations(x, y, temp)
x, y, temp = remove_repeat_coordinates(x, y, temp)

# In[9]:
gx, gy, img = interpolate(x, y, temp, interp_type='linear', hres=75000)
img = np.ma.masked_where(np.isnan(img), img)
view = basic_map(to_proj)
mmb = view.pcolormesh(gx, gy, img, cmap=cmap, norm=norm)
plt.colorbar(mmb, shrink=.4, pad=0, boundaries=levels)

# In[10]:
#Show map of TW with interpolated temps
plt.title("Interpolated Temperatures 17070100")
plt.show()

代码运行没有错误，但我最终得到了一张空的台湾地图。

我超级绝望，任何帮助将不胜感激!!

Answer 1

每当向卡地平面添加数据时，请务必记住定义坐标系。在这种情况下，由于您的数据位于纬度/纬度中，我将首先执行以下操作：

view.pcolormesh(..., transform=ccrs.PlateCarree())

您可能还有兴趣在使用卡托米在其他投影中绘制投影数据中认真使用转换关键字。