I change the script, and try to show it by using basemap.
# 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 some library:
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
Here some fuction:
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
get 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_
add ploygen plot
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)
Get the data:
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)
###########################################
Barnes Interpolation: 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)
Then I make the plot:
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()
It seems the location is wrong!
Maybe I use the Point_Interpolation, which has changed the coordinate.
Here is the import data which I use like this:
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
