-1
我目前正在尝试使用本网站上提供的代码（https://unidata.github.io/MetPy/latest/examples/gridding/Point_Interpolation.html#sphx-glr-examples-gridding-point-interpolation-py）创建一个台湾地图，在Jupyter笔记本上使用线性插值数据。Python - 温度插值
我的数据是这样的形式：
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
，并以这种形式：
#1707lat 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()
的代码运行没有错误但我最终得到一张空白的台湾地图。
我超级绝望，任何帮助将不胜感激！
我不知道这个库，但它看起来像你代表地图本身之外的颜色映射。代码'view.set_extent（[120.5,122.5,24.5,25.5]）'看起来像是在设置坐标。然而，地图的中心'to_proj = ccrs.AlbersEqualArea（central_longitude = 120.0000， central_latitude = 25.0000）'在该框之外...... – agastalver