登录后查看更多精彩内容~
您需要 登录 才可以下载或查看,没有帐号?立即注册
x
本帖最后由 QCD 于 2018-7-7 22:28 编辑
更新:借助Google找到了解决方法
StackOverFlow上已经有人提过类似的问题(https://stackoverflow.com/questions/45042722/cartopy-matplotlib-contourf-map-overriding-data),是由于两幅图都是同一个类型,使用的zorder是相同的,所以需要手动制定zorder,可以参考matplotlib关于zorder的demo(https://matplotlib.org/examples/pylab_examples/zorder_demo.html)
ax.add_feature(adm2_shapes,zorder=1)
axis = ax.contourf(olon,olat,rain_data_new, levels=levels, cmap='viridis',transform=ccrs.PlateCarree(),zorder=5)
最终结果
最近在学习使用Python进行色斑图,在学习了论坛里多位大神的经验后(如非对称的帖子1和他在Github上的MicapsDraw;平流层的萝卜的maskout方法,也参考过chiqu296的code)成功使用Basemap画出的想要的图。
但是考虑到Basemap势必逐渐被CartoPy取代,MetPy和Iris等的例子里面也多使用CartoPy,因此决定尝试CartoPy,在学习了po_po1的方法后实现了成功剪切,如
但是随后在叠加地图的时候遇到了麻烦,如果使用ax.add_feature来加入用于裁剪的shp文件,就会发生只有地图而没有色斑图的现象,如下图所示
最后发现如果将读入的地图数据的alpha值降低后就可以将2者都显示出来,但是色斑图也被叠加了alpha值处理,如下图所示
不知道社区的大神们有没有遇到过类似的问题,或者有没有其他使用cartopy来画色斑图的成功经验?谢谢!
附上使用的代码:
- import numpy as np
- from scipy.interpolate import Rbf
- import pandas as pd
- import matplotlib.pyplot as plt
- from matplotlib.pyplot import colormaps as cmap
- import cartopy.crs as ccrs
- import cartopy.feature as cfeat
- from cartopy.io.shapereader import Reader
- from cartopy.feature import ShapelyFeature
- from cartopy.mpl.patch import geos_to_path
- from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter
复制代码- <font face="微软雅黑">data = pd.read_csv('../rain05.csv')
- # 插值
- lon = data['Lon']
- lat = data['Lat']
- rain_data = data['PRE_Time_0808']
- # [103,110,24.5,29.3]
- olon = np.linspace(103,110,88)
- olat = np.linspace(24.5,29.3,88)
- olon,olat = np.meshgrid(olon,olat)
- # 插值处理
- func = Rbf(lon, lat, rain_data,function='linear')
- rain_data_new = func(olon, olat)
- levels = [0,10,25,50,100,150,200,300,450]
- fname = '../Micaps-shapefiles/Province.shp'
- fname2 = '../Micaps-shapefiles/Province.shp'
- adm1_shapes = ShapelyFeature(Reader(fname).geometries(),ccrs.PlateCarree(), facecolor='White',edgecolor='black')
- adm2_shapes = ShapelyFeature(Reader(fname2).geometries(),ccrs.PlateCarree(),facecolor='White',edgecolor='black',alpha=0.1,linewidth=2)
- #设定裁剪区域
- reader = Reader(fname)
- provinces = reader.records()
- #Micaps地图中Province.shp中可以使用PAC区分不同省份,同理City.shp和County.shp也可使用该关键字
- us_multipoly, = [province.geometry for province in provinces if province.attributes['PAC'] == 520000]
- main_us_geom = sorted(us_multipoly.geoms, key=lambda geom: geom.area)[0]
- us_path, = geos_to_path(main_us_geom)
- #画图部分
- fig = plt.figure(figsize=(12,9))
- ax = fig.add_subplot(1, 1, 1, projection=ccrs.PlateCarree())
- plate_carre_data_transform = ccrs.PlateCarree()._as_mpl_transform(ax)
- axis = ax.contourf(olon,olat,rain_data_new, levels=levels, cmap='viridis',transform=ccrs.PlateCarree(),zorder=5)
- for collection in axis.collections:
- collection.set_clip_path(us_path, plate_carre_data_transform)#设置显示区域
- ax.add_feature(adm2_shapes,zorder=1)
- # ax.set_extent([103,110,24.5,29.3])</font>
复制代码
|