Python可视化 | 三维图形迁移

郭好奇同学

发布于 2021-05-28 17:22:14

1.8K0

发布于 2021-05-28 17:22:14

文章被收录于专栏：好奇心Log

本节提要：通过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()

本文参与腾讯云自媒体同步曝光计划，分享自微信公众号。

原始发表：2021-05-17，如有侵权请联系 cloudcommunity@tencent.com 删除

迁移