我更改了脚本，并尝试使用底图显示它。
# Copyright (c) 2016 MetPy Developers.
# Distributed under the terms of the BSD 3-Clause License.
# SPDX-License-Identifier: BSD-3-Clause
"""
[Point_Interpolation](https://unidata.github.io/MetPy/latest/examples/gridding/Point_Interpolation.html?highlight=basic_map)
"""
导入一些库：
import cartopy.crs as ccrs
import cartopy.feature as cfeature
from matplotlib.colors import BoundaryNorm
import matplotlib.pyplot as plt
from matplotlib import cm as CM
from mpl_toolkits.basemap import Basemap
from matplotlib.patches import Polygon
from metpy.cbook import get_test_data
from metpy.gridding.gridding_functions import (interpolate, remove_nan_observations,
remove_repeat_coordinates)
import pandas as pd
import numpy as np
import numpy.ma as ma
这里有一些功能：
def basic_map(proj):
"""Make our basic default map for plotting"""
fig = plt.figure(figsize=(15, 10))
#add_metpy_logo(fig, 0, 80, size='large')
view = fig.add_axes([0, 0, 1, 1], projection=proj)
view.set_extent([100, 130, 20, 50])
view.add_feature(cfeature.STATES.with_scale('50m'))
view.add_feature(cfeature.OCEAN)
view.add_feature(cfeature.COASTLINE)
view.add_feature(cfeature.BORDERS, linestyle=':')
return fig, view
获取station_test_data：
def station_test_data(variable_names, proj_from=None, proj_to=None):
#with get_test_data('station_data.txt') as f:
all_data = np.loadtxt('2016072713.000', skiprows=1,
# all_data = np.loadtxt(f, skiprows=1,delimiter=',',
usecols=(0,1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11),
dtype=np.dtype([('stid', '5S'), ('lat', 'f'), ('lon', 'f'),
('aqi', 'f'), ('grade', 'f'),
('pm25', 'f'), ('pm10', 'f'),
('co', 'f'),('no2','f'),('o3','f'),
('o38h', 'f'), ('so2', '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
def remove_inf_x(x, y, z):
x_ = x[~np.isinf(x)]
y_ = y[~np.isinf(x)]
z_ = z[~np.isinf(x)]
return x_, y_, z_
def remove_inf_y(x, y, z):
x_ = x[~np.isinf(y)]
y_ = y[~np.isinf(y)]
z_ = z[~np.isinf(y)]
return x_, y_, z_
添加普洛伊根图
def plot_rec(map,lower_left,upper_left,upper_right,lower_right,text):
lon_poly = np.array([lower_left[0], upper_left[0],upper_right[0], lower_right[0]])
lat_poly = np.array([lower_left[1], upper_left[1],upper_right[1], lower_right[1]])
X, Y = map(lon_poly, lat_poly)
xy = np.vstack([X,Y]).T
poly = Polygon(xy, closed=True,
facecolor='None',edgecolor='red',
linewidth=1.
)
ax, ay = map(lower_left[0],lower_left[1])
plt.text(ax, ay,text,fontsize=6,fontweight='bold',
ha='left',va='bottom',color='k')
plt.gca().add_patch(poly)
获取数据：
from_proj = ccrs.Geodetic()
to_proj = ccrs.AlbersEqualArea(central_longitude=110.0000, central_latitude=32.0000)
x, y, temp = station_test_data('pm25', from_proj, to_proj)
x, y, temp = remove_nan_observations(x, y, temp)
x, y, temp = remove_inf_x(x, y, temp)
x, y, temp = remove_inf_y(x, y, temp)
x, y, temp = remove_repeat_coordinates(x, y, temp)
###########################################
巴恩斯插值：search_radius = 100km min_neighbors = 3
gx, gy, img1 = interpolate(x, y, temp, interp_type='barnes', hres=75000, search_radius=100000)
img1 = np.ma.masked_where(np.isnan(img1), img1)
然后我做剧情：
fig = plt.figure(figsize=(8, 8))
#fig, view = basic_map(to_proj)
m = Basemap(width=8000000,height=7000000,
resolution='l',projection='aea',
lat_1=0.,lat_2=40,lon_0=110,lat_0=20)
#m.shadedrelief()
m.drawparallels(np.arange(20.,40,2.5),linewidth=1, dashes=[4, 2], labels=[1,0,0,0], color= 'gray',zorder=0, fontsize=10)
m.drawmeridians(np.arange(100.,125.,2.),linewidth=1, dashes=[4, 2], labels=[0,0,0,1], color= 'gray',zorder=0, fontsize=10)
m.readshapefile('dijishi_2004','state',color='gray')
levels = list(range(0, 200, 10))
cmap = plt.get_cmap('viridis')
norm = BoundaryNorm(levels, ncolors=cmap.N, clip=True)
mmb = m.pcolormesh(gx, gy, img1, cmap=cmap, norm=norm)
plt.colorbar(label=r'$PM_{2.5}(mu g/m^3 $)')
plt.show()
看来位置不对！
也许我使用Point_Interpolation，它改变了坐标。
这是我像这样使用的导入数据：
99306 31.1654 121.412 NaN NaN NaN 37.000 0.875 10.000 141.000 NaN NaN
99299 31.2036 121.478 NaN NaN 49.000 0.420 18.000 157.000 NaN 16.000
99302 31.1907 121.703 NaN NaN 75.000 112.000 1.571 54.000 167.000 NaN 34.000
99300 31.3008 121.467 NaN 53.000 NaN 0.414 21.000 128.000 NaN 10.000
99304 31.2071 121.577 NaN NaN 20.000 66.000 NaN 20.000 192.000 NaN 28.000
99305 31.0935 120.978 NaN NaN 5.000 0.717 23.000 140.000 NaN 13.000
我的猜测是你需要调整你的网格位置或你的投影。底图将 (0， 0( 定义为地图的左下角，而数据看起来假定 (0， 0( 是中心。