python - 使用 Cartopy 插入城市地区的温度数据

Question

我正在尝试对在由 5 个地点组成的市区观察到的温度数据进行插值。我正在使用 cartopy 来插值和绘制 map ，但是，当我运行脚本时，温度插值没有显示，我只得到带有调色板的市区图层。有人可以帮我解决这个错误吗？ shapefile 的链接是

https://www.dropbox.com/s/0u76k3yegvr09sx/LimiteAMG.shp?dl=0

https://www.dropbox.com/s/yxsmm3v2ey3ngsp/LimiteAMG.cpg?dl=0

https://www.dropbox.com/s/yx05n31dfkggbb6/LimiteAMG.dbf?dl=0

https://www.dropbox.com/s/a6nk0xczgjeen2d/LimiteAMG.prj?dl=0

https://www.dropbox.com/s/royw7s51n2f0a6x/LimiteAMG.qpj?dl=0

https://www.dropbox.com/s/7k44dcl1k5891qc/LimiteAMG.shx?dl=0

数据

        Lat   Lon       tmax
 0   20.8208 -103.4434  22.8
 1   20.7019 -103.4728  17.7
 2   20.6833 -103.3500  24.9
 3   20.6280 -103.4261   NaN
 4   20.7205 -103.3172  26.4
 5   20.7355 -103.3782  25.7
 6   20.6593 -103.4136   NaN
 7   20.6740 -103.3842  25.8
 8   20.7585 -103.3904   NaN
 9   20.6230 -103.4265   NaN
 10  20.6209 -103.5004   NaN
 11  20.6758 -103.6439  24.5
 12  20.7084 -103.3901  24.0
 13  20.6353 -103.3994  23.0
 14  20.5994 -103.4133  25.0
 15  20.6302 -103.3422   NaN
 16  20.7400 -103.3122  23.0
 17  20.6061 -103.3475   NaN
 18  20.6400 -103.2900  23.0
 19  20.7248 -103.5305  24.0
 20  20.6238 -103.2401   NaN
 21  20.4753 -103.4451   NaN

代码:

import cartopy
import cartopy.crs as ccrs
from matplotlib.colors import BoundaryNorm
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import cartopy.io.shapereader as shpreader

from metpy.calc import get_wind_components
from metpy.cbook import get_test_data
from metpy.gridding.gridding_functions import interpolate,      remove_nan_observation 
from metpy.plots import add_metpy_logo
from metpy.units import units

to_proj = ccrs.PlateCarree()


data=pd.read_csv('/home/borisvladimir/Documentos/Datos/EMAs/EstacionesZMG/RedZMG.csv',usecols=(1,2,3),names=['Lat','Lon','tmax'],na_values=-99999,header=0)


fname='/home/borisvladimir/Dropbox/Diversos/Shapes/LimiteAMG.shp'
adm1_shapes = list(shpreader.Reader(fname).geometries())


 lon = data['Lon'].values
 lat = data['Lat'].values
 xp, yp, _ = to_proj.transform_points(ccrs.Geodetic(), lon, lat).T


 x_masked, y_masked, t = remove_nan_observations(xp, yp, data['tmax'].values)

 #Interpola temp usando Cressman
tempx, tempy, temp = interpolate(x_masked, y_masked, t, interp_type='cressman', minimum_neighbors=3, search_radius=400000, hres=35000)
temp = np.ma.masked_where(np.isnan(temp), temp)

levels = list(range(-20, 20, 1))
cmap = plt.get_cmap('viridis')
norm = BoundaryNorm(levels, ncolors=cmap.N, clip=True)

fig = plt.figure(figsize=(15, 10))
view = fig.add_subplot(1, 1, 1, projection=to_proj)

view.add_geometries(adm1_shapes, ccrs.PlateCarree(),edgecolor='black',  facecolor='white', alpha=0.5)


view.set_extent([-103.8, -103, 20.3, 21.099 ], ccrs.PlateCarree())

ZapLon,ZapLat=-103.50,20.80
GuadLon,GuadLat=-103.33,20.68
TonaLon,TonaLat=-103.21,20.62
TlaqLon,TlaqLat=-103.34,20.59
TlajoLon,TlajoLat=-103.44,20.47

plt.text(ZapLon,ZapLat,'Zapopan',transform=ccrs.Geodetic())
plt.text(GuadLon,GuadLat,'Guadalajara',transform=ccrs.Geodetic())
plt.text(TonaLon,TonaLat,'Tonala',transform=ccrs.Geodetic())
plt.text(TlaqLon,TlaqLat,'Tlaquepaque',transform=ccrs.Geodetic())
plt.text(TlajoLon,TlajoLat,'Tlajomulco',transform=ccrs.Geodetic())

mmb = view.pcolormesh(tempx, tempy, temp,transform=ccrs.PlateCarree(),cmap=cmap, norm=norm)
plt.colorbar(mmb, shrink=.4, pad=0.02, boundaries=levels)
plt.show()

Answer 1

问题出在对 MetPy 的 interpolate 函数的调用中。使用 hres=35000 的设置，它生成一个间隔为 35km 的网格。但是，您的数据点似乎比这更紧密；一起，这导致生成的网格只有两个点，如下图红点所示(黑点是具有非屏蔽数据的原始站点):

结果是它只为网格创建了两个点，这两个点都在你的数据点的边界之外；因此这些点最终被掩盖了。相反，如果我们将 hres 设置为更低的值，比如 5km(即 5000)，那么会出现一个更合理的结果: