前往小程序,Get更优阅读体验!
立即前往
首页
学习
活动
专区
工具
TVP
发布
社区首页 >专栏 >Python可视化 | 三维图形迁移

Python可视化 | 三维图形迁移

作者头像
郭好奇同学
发布2021-05-28 17:22:14
1.7K0
发布2021-05-28 17:22:14
举报
文章被收录于专栏:好奇心Log好奇心Log

本节提要:通过collection功能的开发实现图形的迁移。



在前面推送中我们提到了通过collection功能而在3D地图中添加地图的方法,也短暂提到了栅格与填色两种图形样式的降维方法。但是从matplotlib这两个函数的底层有一定的局限性,比如下面这两张图的侧面填色就无法绘出:

前一张图只能画最上面的等值线填色和地图,下面这张的栅格也是无法绘制出来的,只能画地图。所以通过相同的collection办法,我们来实现图形的迁移。

一、Axes子图平面pcolormesh的迁移

代码语言:javascript
复制
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()
代码语言:javascript
复制
ap=ax[0].pcolormesh(lon,lat,data,cmap='Spectral_r')

这一句,在第一张子图上绘制了一个pcolormesh,并将它返回的图形几何省称为ap。

代码语言:javascript
复制
datas=data.flatten()
paths=ap.get_paths()

前一句对data降维,后一句获取ap的path。

代码语言:javascript
复制
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的迁移

代码语言:javascript
复制
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的叠加

代码语言:javascript
复制
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的多面迁移

代码语言:javascript
复制
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多面迁移

代码语言:javascript
复制
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则会添加不出来,还没有找到适合的方法来解决这个问题,如果有大神希望能指导一二:

六、仿制一张图

代码语言:javascript
复制
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 删除

本文分享自 好奇心Log 微信公众号,前往查看

如有侵权,请联系 cloudcommunity@tencent.com 删除。

本文参与 腾讯云自媒体同步曝光计划  ,欢迎热爱写作的你一起参与!

评论
登录后参与评论
0 条评论
热度
最新
推荐阅读
相关产品与服务
图像处理
图像处理基于腾讯云深度学习等人工智能技术,提供综合性的图像优化处理服务,包括图像质量评估、图像清晰度增强、图像智能裁剪等。
领券
问题归档专栏文章快讯文章归档关键词归档开发者手册归档开发者手册 Section 归档