Я изменяю сценарий и пытаюсь показать его с помощью базовой карты.
# 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
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_
добавить плойгеновский участок
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 = 100 км 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