我更改了脚本，并尝试使用底图显示它。
# 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_
添加 ploygen 图
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 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 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 NaN 5.000 0.717 23.000 140.000 NaN 13.000