6.Matplotlib绘制高级图像:Basemap和Mplot3D

Matplotlib绘制高级图像

本节将讲解的Matplotlib绘制高级图像主要包含两种图像:

1. 三维图像
2. 地理图像

Matplotlib绘制三维图像直接使用自带的mplot3d工具箱来绘制,Matplotlib绘制地理图像则要使用Basemap工具箱

由于Basemap工具箱不是自带的,因此需要额外下载

Matplotlib绘制三维图像

Matplotlib原本只能绘制二维图像,但是随着版本的迭代逐渐开始可以绘制三维图像了

绘制三维图像的工具箱是mplot3d,我们需要从mpl_toolkits导入这个工具箱,实际上所有基于Matplotlib开发的高级功能库都位于这个位置,所以称其为工具箱

我们首先需要导入这个工具箱

from mpl_toolkits import mplot3d

指定projection参数创建三维子图

和二维图像的情况一样,三维图像也需要三维子图承载,为此,我们导入mplot3d工具箱之后,在任何可以创建二维Axes对象的函数或者方法中指定projection参数为3d后即可创建一个三维子图(Axes3D对象)

from mpl_toolkits import mplot3d

Fig=plt.figure()
Axes3D=plt.axes(projection='3d')
print(type(Axes3D))
plt.show()
>>>
<class 'matplotlib.axes._subplots.Axes3DSubplot'>

这个网格图可以用鼠标来进行交互

Axes3D.plot3D()与Axes3D.scatter3D()绘制三维点与线

最基本的三维图像是(x,y,z)三维坐标点构成的线图与散点图.我们分别可以使用Axes3D对象的plot3D()和scatter3D()方法来绘制

from mpl_toolkits import mplot3d

z=np.linspace(start=0,stop=15,num=1000)
x=np.sin(z)
y=np.cos(z)

z_scatter=15*np.random.random(100)
x_scatter=np.sin(z_scatter)+0.1*np.random.randn(100)
y_scatter=np.cos(z_scatter)+0.1*np.random.randn(100)

Fig=plt.figure()
Axes3D=plt.axes(projection='3d')
Axes3D.plot3D(x,y,z,'gray')
Axes3D.scatter3D(x_scatter,y_scatter,z_scatter,c=z_scatter,cmap='Greens')
plt.show()

mplot3d自动的会将一些点设置为透明的,来表现整个图像的立体感,我们通过交互拖拽图像就能够得到透明的点的变化

Axes3D.contour3D()绘制三维等高线图

前面介绍过将三维图像投影到二维图像上绘制等高线图,我们其实可以通过mplot3d直接绘制三维等高线图

和Axes.contour()方法一样,Axes3D.contour3D()对象我们传入的xoy面的数据必须是二维网格数据

x=np.linspace(start=-6,stop=6,num=30)
y=np.linspace(start=-6,stop=6,num=30)

x_grid,y_grid=np.meshgrid(x,y)
z_grid=np.sin(np.sqrt(x_grid**2+y_grid**2))

Fig=plt.figure()
Axes3D=plt.axes(projection='3d')

Axes3D.contour3D(x_grid,y_grid,z_grid,50,cmap='RdGy')
Axes3D.set_xlabel('x')
Axes3D.set_ylabel('y')
Axes3D.set_zlabel('z')

plt.show()

由于这个角度不是观察的最佳视角

我们可以通过Axes3D.view_init()方法来调整视角

x=np.linspace(start=-6,stop=6,num=30)
y=np.linspace(start=-6,stop=6,num=30)

x_grid,y_grid=np.meshgrid(x,y)
z_grid=np.sin(np.sqrt(x_grid**2+y_grid**2))

Fig=plt.figure()
Axes3D=plt.axes(projection='3d')

Axes3D.contour3D(x_grid,y_grid,z_grid,50,cmap='RdGy')
Axes3D.set_xlabel('x')
Axes3D.set_ylabel('y')
Axes3D.set_zlabel('z')

Axes3D.view_init(60,35)				#调整俯视角为60度,方位角是35度

plt.show()

Axes3D.plot_wireframe()绘制线框图

线框图和下面将要讲解的曲面图都是二维平面直接映射到三维空间中的曲面.我们可以使用plot_wireframe()

x=np.linspace(start=-6,stop=6,num=30)
y=np.linspace(start=-6,stop=6,num=30)

x_grid,y_grid=np.meshgrid(x,y)
z_grid=np.sin(np.sqrt(x_grid**2+y_grid**2))

Fig=plt.figure()
Axes3D=plt.axes(projection='3d')

Axes3D.plot_wireframe(x_grid,y_grid,z_grid,color='black')
Axes3D.set_title('wireframe')

plt.show()

Axes3D.plot_surface()绘制曲面图

和前面的绘制线框图类似,我们为plot_surface传入三个网格数据即可绘制

x=np.linspace(start=-6,stop=6,num=30)
y=np.linspace(start=-6,stop=6,num=30)

x_grid,y_grid=np.meshgrid(x,y)
z_grid=np.sin(np.sqrt(x_grid**2+y_grid**2))

Fig=plt.figure()
Axes3D=plt.axes(projection='3d')

Axes3D.plot_surface(x_grid,y_grid,z_grid,rstride=1,cstride=1,cmap='viridis')
Axes3D.set_title('surface')

plt.show()

此外,我们实际上也可以使用柱坐标来绘图,此时只需要将x和y转变为r和theta

r=np.linspace(start=0,stop=6,num=20)
theta=np.linspace(start=0,stop=1.7*np.pi,num=50)

r_grid,theta_grid=np.meshgrid(r,theta)
x_grid=r_grid*np.sin(theta_grid)
y_grid=r_grid*np.cos(theta_grid)
z_grid=np.sin(np.sqrt(x_grid**2+y_grid**2))

Fig=plt.figure()
Axes3D=plt.axes(projection='3d')

Axes3D.plot_surface(x_grid,y_grid,z_grid,rstride=1,cstride=1,cmap='RdGy',edgecolor='none')
plt.show()

我们通过定义域的转变的方式来实现了对曲面的剖分

Axes3D.plot_trisurf()拟合曲面

上各种绘制曲面的方法,要么需要笛卡尔坐标,要么使用极坐标,并且需要对每个数据点生成网格,而后才能后进行绘图曲面

但在进行正常工作的时候,我们往往可能遇到的情况是由多个数据点而非严格的,均匀的网格数据

为此,我们想要对这些数据点所在的曲面进行拟合,这就需要Axes3D对象的plot_trisurf()方法来绘制

我们首先创建一些随机点,并观察这些点的分布来对数据的分布具有一个大概的认识

# 随机生成点
theta=2*np.pi*np.random.random(1000)
r=6*np.random.random(1000)
x=np.ravel(r*np.sin(theta))					#使用np.ravel来确保x是一维的数组
y=np.ravel(r*np.cos(theta))					#同上,确保y是一维的数组
z=np.sin(np.sqrt(x**2+y**2))

Fig=plt.figure()
Axes3D=plt.axes(projection='3d')

Axes3D.scatter3D(x,y,z)

plt.show()

接下来我们调用plot_trisurf()来进行曲面拟合

# 随机生成点
theta=2*np.pi*np.random.random(1000)
r=6*np.random.random(1000)
x=np.ravel(r*np.sin(theta))
y=np.ravel(r*np.cos(theta))
z=np.sin(np.sqrt(x**2+y**2))

Fig=plt.figure()
Axes3D=plt.axes(projection='3d')

Axes3D.plot_trisurf(x,y,z,cmap='viridis',edgecolor='none')

plt.show()

案例:莫比乌斯环

下面我们将使用上面讲解的内容来绘制一个莫比乌斯环

绘制莫比乌斯环的关键就是确定每个点的三个坐标信息.

莫比乌斯环的纸带我们可以理解有两种旋转关系,一种是围绕z轴进行旋转,另外一种是围绕自己的坐标轴旋转

借此,我们可以确定莫比乌斯环上每个点的坐标

r = 1 + w * np.cos(phi)

x = np.ravel(r * np.cos(theta))
y = np.ravel(r * np.sin(theta))
z = np.ravel(w * np.sin(phi))

from matplotlib.tri import Triangulation
tri = Triangulation(np.ravel(w), np.ravel(theta))

ax = plt.axes(projection='3d')
ax.plot_trisurf(x, y, z, triangles=tri.triangles,
                cmap='viridis', linewidths=0.2)

ax.set_xlim(-1, 1)
ax.set_ylim(-1, 1)
ax.set_zlim(-1, 1)

plt.show()

Matplotlib使用Basemap绘制地图图像

Basemap是Matplotlib的一个第三方开发的工具箱.

Basemap用起来会比较笨重,但是对于Python用户来来说还是很方便的

下面我们首先需要安装Basemap,然后在此基础上进行绘制

安装Basemap

Basemap实际不仅仅基于Matplotlib,还基于pyproj等地图管理库

因此想要正确的安装并能运行Basemap,安装的各个依赖库的版本必须要对,否则就要用pip重装

下面首先指出各个依赖库的版本

库名	适用版本
six	1.15.0
matplotlib	2.2.0
geos	0.2.2
pyproj	1.9.6

注意.pyproj需要使用1.9.6,否则会报错(我一开始使用的2.6.1的pyproj),matplotlib需要使用2.x版本,我一开始使用的是matplotlib3.x,这样就会报错

geos和basemap的安装步骤参考csdn博文即可

Basemap的基础概念

蓝色弹珠

下面我们将创建一个地球的蓝色弹珠图像，来讲解Basemap中的基础概念

from mpl_toolkits.basemap import Basemap

Fig=plt.figure(figsize=(8,8))
m=Basemap(projection='ortho',resolution=None,lat_0=50,lon_0=-100)

m.bluemarble(scale=0.5)
print(type(m))

plt.show()
>>>
<class 'mpl_toolkits.basemap.Basemap'>

这里我们首先从basemap库中导入了Basemap对象，Basemap对象是Basemap库中的核心对象，就像Pandas中的Series和DataFrame对象，Numpy中的Ndarry对象一样

Basemap库内置了全球的地理数据，包括经度(Longitude)，纬度（Latitude)以及对应的图像，因此我们如果想要显示某个地区的图像，那么首先就需要指定显示的地点的经纬度以及图像的范围，所以我们实例化一个Basemap对象其实就是确定其地理数据范围，注意Basemap仅仅作为地理数据范围，而非图像

显示出具体的图像需要我们调用Basemap的不同方法来绘制图像，例如我们这里调用bluemarble方法来将这些数据可视化

此外，由于地球是一个三维球体，因此其表面的地理图像自然分布在球面上，所以这个时候就有一个至关重要的问题：如何将三维图像(地球上的地理图像)表现在二位图像上

这个问题最早出现在绘图领域，当时的人们面临如何绘制地图

解决这个问题的关键就是投影，我们可以通过各种不同的方式将三维球面投影到二维平面，从而使得二维图像能够表示出地理图像

不同的投影方式的优点不同，因此在我们关注于不同的重点时，往往会选择不同的投影方式来绘制地图

因此，我们实例化一个Basemap对象的时候，不仅需要指定绘制地点的经度和纬度，还需要指定投影方式，至于图形的大小，可以直接使用缺失值也可以指定

这里我们首先创建了一个空白的Figure对象，接下来我们指定了投影的方式为正射投影(ortho)，这样便于我们从高空观察地球整体的图像

接下来我们指定显示的地点的经纬度为-100和50

此外，Basemap对象还具有许多方法，这里我们调用了bluemarble方法来绘制地球的蓝色弹珠图像

球面坐标

下面我们将绘制西安的位置来进一步掌握Basemap的基础概念

from mpl_toolkits.basemap import Basemap


Fig=plt.figure(figsize=(8,8))
m=Basemap(projection='lcc',resolution=None,\
    lat_0=33.8,lon_0=108.4,\
    width=8E6,height=8E6)

Longitude_XiAn,Latitude_XiAn=m(108.4,33.8)

m.etopo(scale=0.5,alpha=0.5)
plt.plot(Latitude_XiAn,Longitude_XiAn,'or',markersize=5)
plt.text(Latitude_XiAn+8E4,Longitude_XiAn+8E4,"Xi'an")

plt.show())

这里我们首先实例化一个Basemap对象，确定其数据范围为以西安为中心(其实就是代码中的经纬度)，投影方式为兰伯特等角圆锥投影(Lambert conforml conic projection, lcc),这种投影方式用来显示局部的地区效果会很好，具体将在后面讲解

而实际上，我们通过各种方法可视化的地理图像(这里我们使用Etopo方法，这是一种能够陆地与海底地形特征的绘图方式)其实是Matplotlib中的Axes对象下的一个子类AxesImage对象

但是作为子类，AxesImage对象并不能像Axes对象一样调用所有方法,例如Axes.plot()方法，所以我们不能够使用AxesImage.plot()绘制散点图来,进而添加西安这个点

因此实际上,我们调用Basemap的可视化方法绘制的其实是背景,即地图背景

还需要说明的是，Basemap可视化得到的AxesImage对象是基于球面坐标的，因此我们想要为AxesImage对象添加点的时候，必须先要将点的经纬度转换为球面坐标，这样才能够正确的绘制，即代码中的Longitude_XiAn,Latitude_XiAn=m(108.4,33.8)语句

下面我们来检验一下上面的内容

from mpl_toolkits.basemap import Basemap


Fig=plt.figure(figsize=(8,8))
m=Basemap(projection='lcc',resolution=None,\
    lat_0=33.8,lon_0=108.4,\
    width=8E6,height=8E6)

Longitude_XiAn,Latitude_XiAn=m(108.4,33.8)

AxesImage_XiAn=m.etopo(scale=0.5,alpha=0.5)
plt.plot(Latitude_XiAn,Longitude_XiAn,'or',markersize=5)
plt.text(Latitude_XiAn+8E4,Longitude_XiAn+8E4,"Xi'an")


print(type(AxesImage_XiAn))
print(type(Longitude_XiAn),Longitude_XiAn)
print(type(Latitude_XiAn),Latitude_XiAn)
>>>
<class 'matplotlib.image.AxesImage'>
<type 'float'>　 3999999.9999999986
<type 'float'>　 3999999.999999997

Basemap的投影方式

Basemap中具有三十多种投影方式,下面我们将讲解Basemap中的常见的投影方式

我们首先定义一个可以绘制经纬线的函数,在下面每种投影风格的演示中,我们将使用这个函数来查看效果\

def DrawMap(m,scale=0.2):

    from itertools import chain
    
    #绘制地貌晕染图
    m.shadedrelief(scale=scale)

    #绘制经纬度
    Latitude=m.drawparallels(np.linspace(start=-90,stop=90,num=13))
    Longitude=m.drawmeridians(np.linspace(start=-180,stop=180,num=13))

    #生成经纬线数据以便于设置经纬线样式
    LatitudeLines=chain(*(tup[1][0] for tup in Latitude.items()))
    LongitudeLines=chain(*(tup[1][0] for tup in Longitude.items()))
    AllLines=chain(LatitudeLines,LongitudeLines)

    #设置经纬线样式
    for line in AllLines:
        line.set(linestyle='-',alpha=0.3,color='w')

    plt.show()

这里我们首先导入了itertools这个库,itertool的主要用途就是进行高效的迭代,主要用于后面我们绘制经纬线时使用

我们调用了Basemap的shadedrelief方法来绘制地貌晕染图.在前面讲过,我们如果需要在Basemap绘制的地形图上绘制点或线,那么首先就需要将点和线(其实就是以点为元素的列表)的直角坐标转换为球坐标.

这里我们使用的Basemap对象的drawparallels和drawmeridian方法来绘制纬度线和经度线,注意,通过这两种方法得到的返回值其实是一个字典,字典的键是plt.Line2D对象

接下来我们使用itertools库的chain方法,chain方法是将多个对象串联成一个新的可迭代对象的方法,这里我们调用字典对象的items()方法将所有字典的键值对(元组形式)的键取出,传给chain方法.但是由于取出的所有的键是一个不定长的数组,因此要加上*

最后,通过chain方法生成的新的可迭代对象中的所有的元素都是plt.Line2D对象,包含所有的经纬线,我们迭代设置样式即可

圆柱投影　Cylindrical Projection

圆柱投影是最简单的地图投影类型,我们将经度线和纬度线分别映射成水平线与竖直线,从而将球形地图映射为一个矩形图像

但是通过圆柱投影得到的地图在赤道部分的表现会很好,但是在两个极地部分的表现会因为产生了畸变而很差

此外,根据经纬线的排布方式,我们还有等距圆柱投影(cyl),等积圆柱投影(cea,cylindrical equal-area),墨卡托圆柱投影(merc, Mercator)

Fig=plt.figure(figsize=(8,6),edgecolor='w')
m=Basemap(projection='cyl',resolution=None,llcrnrlat=-90,urcrnrlat=90,llcrnrlon=-180,urcrnrlon=180)
DrawMap(m)

这里我们调用了Basemap对象的llcrnrlat,urcrnrlat,llcrnrlon,urcrnrlon参数,来分别指定左下角(llcrnr,lower left corner)和右上角(urcrnr, up right corner)的经度(lon)和纬度(lat)设置

可以看到格林兰岛的畸变非常严重

伪圆柱投影 Pseudo-Cylindrical Projection

伪圆柱投影是在圆柱投影的基础上进行了改进,经度线不必是竖直的,而可以是弯曲的,这样将使得南北极地区区域的地形更加真实

例如摩尔威德投影(Moollweide,moll),他所有的经度线都是椭圆弧线

其他的伪圆柱投影还有正弦投影(sinusoidal,sinu),罗宾森投影(Robinson,robin)

Fig=plt.figure(figsize=(8,6),edgecolor='w')
m=Basemap(projection='moll',resolution=None,lat_0=0,lon_0=0)
DrawMap(m)

这里我们指定lat_0和lon_0来指定地区中心的纬度和经度

我们可以看到格林兰岛的外观稍微好点了

透视投影 Perspective Projection

透视投影是从某一个头十点对地球进行透视获得的投影,就好像站在地球晶体空的某一点给地球照相一样

一个典型的透视投影就是我们前面使用的正射投影(orthographic, ortho),即从无限远处观察地球的一侧

此外还有球心投影(gnomonic,gnom)和球极平面投影(stereographic,stere)

Fig=plt.figure(figsize=(8,6),edgecolor='w')
m=Basemap(projection='ortho',resolution=None,lat_0=50,lon_0=0)
DrawMap(m)

圆锥投影 Conic Projection

圆锥投影是先将地图投影成一个圆锥体,然后在将其展开,这样将获得非常好的局部效果,但是距离圆锥顶点较远的区域会发生严重的畸变

典型的例子就是兰伯特等角圆锥投影(Lambert conforml conic projection, lcc),它指定两个标准的纬线作为构成圆锥(分别是lat_1与lat_2)

此外,常用的圆锥投影还有等距圆锥投影(equidistant,eqdc)和阿尔伯斯等积圆锥(Albers equal-area,aea)

Fig=plt.figure(figsize=(8,6),edgecolor='w')
m=Basemap(projection='lcc',resolution=None,lat_0=50,lon_0=0,lat_1=45,lat_2=55,width=1.6E7,height=1.2E7)

DrawMap(m)

其他投影类型

实际上在Basemap程序包的介绍文档里还有更多关于投影类型的介绍,自行查阅即可

Basemap绘制地图背景

前面介绍,我们可以调用Basemap对象的bluemarble()方法和shadedrelief()绘制地球的投影,用drawparallels()绘制纬线,drawmeridians()方法绘制经线

其实我们还有许多地理图像可视化方法,下面将介绍一些方法

• 物理边界与水域
  • drawcoastlines(): 绘制大陆海岸线
  • drawlsmask():为陆地和海洋设置填充色,从而可以在陆地或海洋投影其他图像
  • drawmapboundary(): 绘制地图边界,包括为海洋填充颜色
  • drawrivers(): 绘制河流
  • fillcontinents(): 使用一种颜色填充大陆,可选参数用另外一种参数填充湖泊
• 政治边界
  • drawcountries(): 绘制国界线
  • drawstates(): 绘制美国州界线
  • drawcounties(): 绘制美国县界线
• 地图特性
  • drawgreatcircle(): 在两点间绘制圆
  • drawparallels(): 绘制纬线
  • drawmeridians(): 绘制经线
  • drawmapscale(): 绘制现行比例尺
• 地球图像
  • bluemarble():绘制NASA蓝色弹珠地球投影
  • shadedrelief(): 在地图上绘制地埋晕染图
  • etopo(): 绘制地形晕染图
  • warpimage(): 将用户提供的图像投影到地图上

此外,如果我们使用边界类图像绘制方法,就必须指定分辨率.我们可以指定resolution参数来设置分辨率,分别是

值	说明
c	原始分辨率
l	低分辨率
i	中分辨率
h	高分辨率
f	全画质分辨率

我们绘制一个海岸线为例

for i,res in enumerate(['l','h']):

    m=Basemap(projection='ortho',lat_0=57.3,lon_0=-6.2,width=90000,height=12000,resolution=res,ax=Axes[i])
    m.fillcontinents(color='#FFDDCC',lake_color='#DDEEFF')
    m.drawmapboundary(fill_color='#DDEEFF',)
    m.drawcoastlines()
    Axes[i].set_title("Resolution='{0}'".format(res))

plt.show()

这里我们使用Python自带的enumerate()函数,来将遍历的对象的元素和索引结合起来,类似于Series

Basemap在地图上绘制数据

最后,Basemap中最实用的功能就是以地图为背景绘制各种数据

我们可以使用任意的plt函数来绘制图形和文字,当然,我们需要将直角坐标转化为球面坐标

实际上Basemap对象中有许多方法也可以绘制地图,只不过多了latlon参数,如果我们将其设置为True,就表示使用原来的经纬度坐标,而非直角坐标

Basemap对现代格部分绘图方法如下:

• contour() / contourf() : 绘制等高线 / 填充等高线
• imshow() : 绘制图像
• pcolor() / pcolormesh() : 绘制带规则 / 不规则网格的伪彩图
• plot() :绘制线条和散点
• scatter() : 绘制带标签的点
• quiver() : 绘制箭头
• barbs() : 绘制风羽
• drawgreatcircle() : 绘制大圆圈

下面我们用前面使用过的美国加州城市数据为例,绘制散点图

cities=pd.read_csv('california_cities.csv')

cities_lat=cities['latd'].values
cities_lon=cities['longd'].values
cities_population=cities['population_total'].values
cities_area=cities['area_total_km2'].values

Fig=plt.figure(figsize=(8,8))

m=Basemap(projection='lcc',resolution='h',lat_0=37.5,lon_0=-119,width=1E6,height=1.2E6)
m.shadedrelief()
m.drawcoastlines(color='gray')
m.drawcountries(color='gray')

m.scatter(cities_lon,cities_lat,latlon=True,c=np.log10(cities_population),s=cities_area,cmap='Reds',alpha=0.5)

plt.colorbar(label=r"$\log_{10}({\rm population})$")
plt.clim(3,7)
for a in [100,300,500]:
    plt.scatter([],[],c='red',alpha=0.5,s=a,label=str(a)+'km$^2$')
plt.legend(scatterpoints=1,frameon=False,labelspacing=1,loc='lower left')

plt.show()

我们可以发现,加州所有的城市都沿着平坦的中央谷地的高速公路延伸,几乎完全避开了加州的山区地带