气象绘图 |
您所在的位置：网站首页 › xlim在r › 气象绘图 |
气象绘图
|
本节提要：通过collection功能的开发实现图形的迁移。 在前面推送中我们提到了通过collection功能而在3D地图中添加地图的方法，也短暂提到了栅格与填色两种图形样式的降维方法。但是从matplotlib这两个函数的底层有一定的局限性，比如下面这两张图的侧面填色就无法绘出：   前一张图只能画最上面的等值线填色和地图，下面这张的栅格也是无法绘制出来的，只能画地图。所以通过相同的collection办法，我们来实现图形的迁移。 一、Axes子图平面pcolormesh的迁移 import matplotlib.pyplot as plt import cartopy.crs as ccrs import cartopy.io.shapereader as spder import numpy as np import itertools from matplotlib.collections import LineCollection, PolyCollection plt.rcParams['font.sans-serif']=['SimHei']#显示中文 plt.rcParams['axes.unicode_minus']=False shp=r'E:\家园\利川市地图\利川.shp' fig=plt.figure(figsize=(4,2)) fig,ax=plt.subplots(ncols=2,subplot_kw={'projection':ccrs.PlateCarree()},dpi=500) for i in range(2): ax[i].add_geometries(spder.Reader(shp).geometries(),crs=ccrs.PlateCarree(),facecolor='none',edgecolor='k') ax[i].set_extent([108.3,109.35,29.7,30.7]) x=np.arange(108.3,109.35,0.05) y=np.arange(29.7,30.7,0.05) lon,lat=np.meshgrid(x,y) data=(lon-(109.35-108.3))**2+(lat-(29.7-30.7))**2 ap=ax[0].pcolormesh(lon,lat,data,cmap='Spectral_r') fc=fig.colorbar(ap,ax=[ax[0],ax[1]],shrink=0.5) fc.ax.tick_params(labelsize=4) datas=data.flatten() paths=ap.get_paths() concat = lambda iterable: list(itertools.chain.from_iterable(iterable)) polys = concat(path.to_polygons() for path in paths) lc=PolyCollection(polys,cmap='Spectral_r',edgecolor='none',closed=False,array=datas) ax[1].add_collection(lc) plt.show() ap=ax[0].pcolormesh(lon,lat,data,cmap='Spectral_r')这一句，在第一张子图上绘制了一个pcolormesh，并将它返回的图形几何省称为ap。 datas=data.flatten() paths=ap.get_paths()前一句对data降维，后一句获取ap的path。 concat = lambda iterable: list(itertools.chain.from_iterable(iterable)) polys = concat(path.to_polygons() for path in paths) lc=PolyCollection(polys, cmap='Spectral_r', edgecolor='none', closed=False, array=datas)变path为polygon再收集为collection和地图的添加步骤类似，唯array这一个参数是用来给polygon上色的。 二、跨越Axes与Axes3D进行collection的迁移 import itertools import pandas as pd from mpl_toolkits.mplot3d import Axes3D import matplotlib.pyplot as plt from matplotlib.collections import LineCollection, PolyCollection import numpy as np from cartopy.io.shapereader import Reader import cartopy.feature from cartopy.mpl.patch import geos_to_path import cartopy.crs as ccrs import numpy as np ########################################################### proj= ccrs.PlateCarree() plt.rcParams['font.sans-serif']=['FangSong']#正常显示中文 fig = plt.figure(figsize=(8,5),dpi=700) ax = Axes3D(fig, xlim=[108.3,109.35], ylim=[29.7,30.7]) ax.set_zlim(0,50) target_projection = ccrs.PlateCarree() ###################################################################### reader = Reader(r'E:\家园\利川市地图\利川.shp') provinces = cartopy.feature.ShapelyFeature(reader.geometries(), proj, edgecolor='k', facecolor='none') #################################################################### geoms = provinces.geometries() geoms = [target_projection.project_geometry(geom, provinces.crs) for geom in geoms] paths = list(itertools.chain.from_iterable(geos_to_path(geom) for geom in geoms)) segments = [] for path in paths: vertices = [vertex for vertex, _ in path.iter_segments()] vertices = np.asarray(vertices) segments.append(vertices) lc = LineCollection(segments, color='black',lw=1) ax.add_collection3d(lc,zs=50) ax.set_xlabel('经度 (E)') ax.set_ylabel('纬度 (N)') ax.set_zlabel('层次') ax.view_init(elev=35,azim=290)#改变绘制图像的视角,即相机的位置,azim沿着z轴旋转，elev沿着y轴 plt.title('pcolormesh迁移',fontsize=20) ax1=fig.add_axes([0,0,0.5,0.5],projection=ccrs.PlateCarree()) ax1.add_geometries(Reader(r'E:\家园\利川市地图\利川.shp').geometries(),crs=ccrs.PlateCarree(),facecolor='none',edgecolor='k') ax1.set_extent([108.3,109.35,29.7,30.7]) x=np.arange(108.3,109.35,0.05) y=np.arange(29.7,30.7,0.05) lon,lat=np.meshgrid(x,y) data=(lon-(109.35-108.3))**2+(lat-(29.7-30.7))**2 ap=ax1.pcolormesh(lon,lat,data,cmap='Spectral_r') datas=data.flatten() paths=ap.get_paths() concat = lambda iterable: list(itertools.chain.from_iterable(iterable)) polys = concat(path.to_polygons() for path in paths) lc=PolyCollection(polys,cmap='Spectral_r',edgecolor='none',closed=False,array=datas) lc1=PolyCollection(polys,cmap='BuGn_r',edgecolor='none',closed=False,array=datas) lc2=PolyCollection(polys,cmap='viridis',edgecolor='none',closed=False,array=datas) ax.add_collection3d(lc,zs=50) ax.add_collection3d(lc1,zs=25) ax.add_collection3d(lc2,zs=0) ax1.set_visible(False)#解除ax1的显示原本生成的图：  解除掉用于生成平面pcolormesh的子图的显示后的图：  三、多层contourf的叠加 import itertools import pandas as pd from scipy.interpolate import Rbf from mpl_toolkits.mplot3d import Axes3D import matplotlib.pyplot as plt from matplotlib.collections import LineCollection, PolyCollection import numpy as np from cartopy.io.shapereader import Reader import cartopy.feature from cartopy.mpl.patch import geos_to_path import cartopy.crs as ccrs import numpy as np ########################################################### proj= ccrs.PlateCarree() plt.rcParams['font.sans-serif']=['FangSong']#正常显示中文 fig = plt.figure(figsize=(8,5),dpi=700) ax = Axes3D(fig, xlim=[108.3,109.35], ylim=[29.7,30.7]) ax.set_zlim(0,50) target_projection = ccrs.PlateCarree() ###################################################################### reader = Reader(r'E:\家园\利川市地图\利川.shp') provinces = cartopy.feature.ShapelyFeature(reader.geometries(), proj, edgecolor='k', facecolor='none') #################################################################### geoms = provinces.geometries() geoms = [target_projection.project_geometry(geom, provinces.crs) for geom in geoms] paths = list(itertools.chain.from_iterable(geos_to_path(geom) for geom in geoms)) segments = [] for path in paths: vertices = [vertex for vertex, _ in path.iter_segments()] vertices = np.asarray(vertices) segments.append(vertices) lc = LineCollection(segments, color='black',lw=1) ax.add_collection3d(lc,zs=0) ax.set_xlabel('经度 (E)') ax.set_ylabel('纬度 (N)') ax.set_zlabel('层次') ax.view_init(elev=35,azim=290)#改变绘制图像的视角,即相机的位置,azim沿着z轴旋转，elev沿着y轴 plt.title('三维contourf平面化',fontsize=20) ax1=fig.add_axes([0,0,0.5,0.5],projection=ccrs.PlateCarree()) ax1.add_geometries(Reader(r'E:\家园\利川市地图\利川.shp').geometries(),crs=ccrs.PlateCarree(),facecolor='none',edgecolor='k') ax1.set_extent([108.3,109.35,29.7,30.7]) x=np.arange(108.3,109.35,0.05) y=np.arange(29.7,30.7,0.05) lon,lat=np.meshgrid(x,y) data=(lon-(109.35-108.3))**2+(lat-(29.7-30.7))**2 ax.contourf(lon,lat,data,cmap='Spectral_r',offset=50) ax.contourf(lon,lat,data,cmap='BuGn_r',offset=25) ax.contourf(lon,lat,data,cmap='viridis',offset=0) ax1.set_visible(False) 四、Axes的pcolormesh的多面迁移 import itertools import pandas as pd from scipy.interpolate import Rbf from mpl_toolkits.mplot3d import Axes3D import matplotlib.pyplot as plt from matplotlib.collections import LineCollection, PolyCollection import numpy as np from cartopy.io.shapereader import Reader import cartopy.feature from cartopy.mpl.patch import geos_to_path import cartopy.crs as ccrs import numpy as np ########################################################### proj= ccrs.PlateCarree() concat = lambda iterable: list(itertools.chain.from_iterable(iterable)) plt.rcParams['font.sans-serif']=['FangSong']#正常显示中文 plt.rcParams['axes.unicode_minus']=False fig = plt.figure(figsize=(8,5),dpi=700) ax = Axes3D(fig, xlim=[-50,50], ylim=[-50,50]) ax.set_zlim(-50,50) ax.set(xticks=np.arange(-50,60,25),yticks=np.arange(-50,60,25),zticks=np.arange(-50,60,25)) ax2=fig.add_axes([0,0,0.5,0.5]) x=np.arange(-50,55,5) y=np.arange(-50,55,5) x,y=np.meshgrid(x,y) z=x**2+y**2 a2=ax2.pcolormesh(x,y,z,shading='auto') ################################################################## zs=z.flatten() paths2=a2.get_paths() polys2 =concat(path.to_polygons() for path in paths2) lc02=PolyCollection(polys2,cmap='Spectral_r',edgecolor='none',closed=False,array=zs) lc03=PolyCollection(polys2,cmap='Greys',edgecolor='none',closed=False,array=zs) lc04=PolyCollection(polys2,cmap='BuGn',edgecolor='none',closed=False,array=zs) ax.add_collection3d(lc02,zdir='y',zs=-51) ax.add_collection3d(lc03,zdir='x',zs=51) ax.add_collection3d(lc04,zdir='z',zs=51) ax2.set_visible(False) ax.set_title('立方pcolocmesh',fontsize=17,pad=-100000000) 通过修改zdir来实现各个面的贴瓷砖。 五、Axes的contourf多面迁移 import itertools import pandas as pd from scipy.interpolate import Rbf from mpl_toolkits.mplot3d import Axes3D import matplotlib.pyplot as plt from matplotlib.collections import LineCollection, PolyCollection from mpl_toolkits.mplot3d.art3d import Poly3DCollection import numpy as np from cartopy.io.shapereader import Reader import cartopy.feature from cartopy.mpl.patch import geos_to_path import cartopy.crs as ccrs import numpy as np ########################################################### proj= ccrs.PlateCarree() concat = lambda iterable: list(itertools.chain.from_iterable(iterable)) plt.rcParams['font.sans-serif']=['FangSong']#正常显示中文 plt.rcParams['axes.unicode_minus']=False fig = plt.figure(figsize=(8,5),dpi=700) ax = Axes3D(fig, xlim=[-50,50], ylim=[-50,50]) ax.set_zlim(-50,50) ax.set(xticks=np.arange(-50,60,25),yticks=np.arange(-50,60,25),zticks=np.arange(-50,60,25)) ax2=fig.add_axes([0,0,0.5,0.5]) x=np.arange(-50,55,5) y=np.arange(-50,55,5) x,y=np.meshgrid(x,y) z=x+np.sin(x)-y-np.tan(y) cs=ax2.contourf(x,y,z) ################################################################## for level,collection in zip(cs.levels,cs.collections): paths=collection.get_paths() polys=concat(path.to_polygons() for path in paths) pc=PolyCollection(polys,facecolor=collection.get_facecolor(),edgecolor='none') ax.add_collection3d(pc,zdir='z',zs=50) for level,collection in zip(cs.levels,cs.collections): paths=collection.get_paths() polys=concat(path.to_polygons() for path in paths) pc=PolyCollection(polys,facecolor=collection.get_facecolor(),edgecolor='none') ax.add_collection3d(pc,zdir='x',zs=50) for level,collection in zip(cs.levels,cs.collections): paths=collection.get_paths() polys=concat(path.to_polygons() for path in paths) pc=PolyCollection(polys,facecolor=collection.get_facecolor(),edgecolor='none') ax.add_collection3d(pc,zdir='y',zs=-50) ax2.set_visible(False) 但是这个办法有个问题，如果有封闭的polygon则会添加不出来，还没有找到适合的方法来解决这个问题,如果有大神希望能指导一二：  六、仿制一张图  import itertools
import pandas as pd
from scipy.interpolate import Rbf
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt
from matplotlib.collections import LineCollection,PolyCollection
import numpy as np
from cartopy.io.shapereader import Reader
import cartopy.feature
from cartopy.mpl.patch import geos_to_path
import cartopy.crs as ccrs
from metpy.interpolate import interpolate_to_grid
concat = lambda iterable: list(itertools.chain.from_iterable(iterable))
from matplotlib.path import Path
from cartopy.mpl.patch import geos_to_path
shp_data=Reader(r'E:\家园\利川市地图\利川.shp')
records=shp_data.records()
for record in records:
path=Path.make_compound_path(*geos_to_path([record.geometry]))
###########################################################
proj= ccrs.PlateCarree()
plt.rcParams['font.sans-serif']=['FangSong']#正常显示中文
fig = plt.figure(figsize=(5,10),dpi=700)
[75,105,25,40]
ax=Axes3D(fig, xlim=[60,105], ylim=[25,48])
ax.set_zlim(1,5)
######################################################################
reader=Reader(r'E:\tibetan\tibetan.shp')
feature=cartopy.feature.ShapelyFeature(reader.geometries(), proj, edgecolor='k', facecolor='none')
filename=r'C:\Users\lenovo\Desktop\利川.xlsx'
df=pd.read_excel(filename)
data=r'C:\Users\lenovo\Desktop\sh1.csv'
datash=pd.read_csv(data)
lon=datash['站点经度']
lat=datash['站点纬度']
data1984=datash['1984']
###########这一步为cressman插值#####################
olon=np.linspace(75,105,60)
olat=np.linspace(25,40,30)
olon,olat=np.meshgrid(olon,olat)
func=Rbf(lon,lat,data1984,function='linear')
data_new=func(olon,olat)
####################预先设置地图的参数######################################
reader = Reader(r'E:\tibetan\tibetan.shp')
provinces = cartopy.feature.ShapelyFeature(reader.geometries(), proj, edgecolor='k', facecolor='none')
target_projection = proj
proj_ax = plt.figure().add_axes([0, 0, 1, 1], projection=proj)
####################################################################
geoms = provinces.geometries()
geoms = [target_projection.project_geometry(geom, provinces.crs)
for geom in geoms]
paths = list(itertools.chain.from_iterable(geos_to_path(geom) for geom in geoms))
segments = []
for path in paths:
vertices = [vertex for vertex, _ in path.iter_segments()]
vertices = np.asarray(vertices)
segments.append(vertices)
lc = LineCollection(segments, color='black',lw=1)
lc1 = LineCollection(segments, color='black',lw=1)
ax.add_collection3d(lc,zs=1)
ax.add_collection3d(lc1,zs=5)
#############################################################################################
records=Reader(r'E:\tibetan\tibetan.shp').records()
for record in records:
path=Path.make_compound_path(*geos_to_path([record.geometry]))
patch=path.to_polygons()
#############################################################################################
ax1=fig.add_axes([0,0,0.5,0.5],projection=ccrs.PlateCarree())
ax1.add_geometries(Reader(r'E:\tibetan\tibetan.shp').geometries(),crs=ccrs.PlateCarree(),facecolor='none',edgecolor='k')
ax1.set_extent([60,105,25,48])
a2=ax1.pcolormesh(olon,olat,data_new,cmap='hot')
zs=data_new.flatten()
paths=a2.get_paths()
polys=concat(path.to_polygons() for path in paths)
lc01=PolyCollection(polys,cmap='hot',edgecolor='none',closed=False,array=zs)
lc02=PolyCollection(polys,cmap='hot',edgecolor='none',closed=False,array=zs)
ax.add_collection3d(lc01,zdir='z',zs=1)
ax.add_collection3d(lc02,zdir='z',zs=5)
lw=0.5
ax.plot([65,65],[24,24],[1,5],c='k',lw=lw)
ax.plot([105,105],[24,24],[1,5],c='k',zorder=10,lw=lw)
ax.plot([105,105],[48,48],[1,5],c='k',zorder=10,lw=lw)
ax.plot([65,65],[48,48],[1,5],c='k',lw=lw)
ax.plot([65,105],[24,24],[1,1],c='k',lw=lw)
ax.plot([65,105],[48,48],[1,1],c='k',lw=lw)
ax.plot([65,105],[24,24],[5,5],c='k',lw=lw)
ax.plot([65,105],[48,48],[5,5],c='k',lw=lw)
ax.plot([65,65],[24,48],[1,1],c='k',lw=lw)
ax.plot([65,65],[24,48],[5,5],c='k',lw=lw)
ax.plot([105,105],[24,48],[1,1],c='k',zorder=10,lw=lw)
ax.plot([105,105],[24,48],[5,5],c='k',zorder=10,lw=lw)
ax.axis('off')
ax1.set_visible(False)
proj_ax.set_visible(False)
ax.text(63,48,5.5,'立方图',fontsize=20)
plt.show()
|
【本文地址】