目录

上一个主题

8. 使用Basemap进行地图可视化

下一个主题

8.2. 绘制地图背景

关注公众号


常见问题

  1. Windows下的安装说明
  2. Jupyter免费在线实验环境
  3. 勘误与补充


>>> from helper import info; info()
页面更新时间: 2020-02-21 18:03:21
操作系统/OS: Linux-4.19.0-8-amd64-x86_64-with-debian-10.3 ;Python: 3.7.3

8.1. 简介

8.1.1. Basemap简介

Matplotlib是Python常用的数据绘制包。它基于NumPy的数组运算功能。 Matplotlib绘图功能强大,可以轻易的画出各种统计图形,比如散点图,条行图,饼图等。 Matplotlib常与Numpy和Scipy相配合,用于许多研究领域。他们是免费工具,但其功能足可以与科研界的大佬Matlab竞争。

Matplotlib中的Basemap它具有专业标准的地图绘制工具。它可以与matplotlib的一般绘图功能相结合,并在地图上绘制数据。 Matplotlib的Basemap工具包是一个用于在Python中的地图上绘制2D数据的库。它在功能上类似于MATLAB地图工具箱,IDL地图工具,GrADS或通用地图工具。 PyNGL和CDAT是在Python中提供类似功能的类库。 Basemap本身不会进行任何绘图,但提供了将坐标转换为25个不同地图投影之一的功能。 Matplotlib也可以用于绘制变换坐标中的轮廓,图像,向量,线或点。 提供了海岸线,河流和政治边界数据集(来自通用地图工具),以及绘制它们的方法。 在Baseap底层使用了GEOS库,用来将海岸线和边界特征剪切到所需的地图投影区域。 Basemap提供读取shapefile的功能。

Basemap适合地球科学家,特别是海洋学家和气象学家的需求。 最初编写Basemap是用来帮助和研究气候和天气预报的,当时CDAT是Python中唯一的用于绘制地图投影数据。 多年来,Basemap的功能随着各个学科(如生物学,地质学和地球物理学)的科学家的要求和贡献的新功能而演变。

Basemap是Matplotlib的一个子包,负责地图绘制。在数据可视化过程中,我们需要将数据在地图上画出来。 比如说我们在地图上画出城市人口,飞机航线,军事基地,矿藏分布等等。这样的地理绘图有助于读者理解空间相关的信息。

Basemap 在 2020 年前随着 Python 2.7 版本一直有更新维护的。2020 年以后 Python 2.7 将停止更新,Basemap 会按照官方计划也迁移到 Cartopy 模块。

8.1.2. 安装

Basemap是Python的软件包,但是目前由于其与操作系统底层的一些类库耦合比较紧密, 无法直接通过 PIP 工具进行安装。

在 Debain Stretch中,可以使用下面的命令进行安装:

# apt install python3-mpltoolkits.basemap

或者,自行下载源代码编译安装。源代码可以直接下载压缩包,或者使用 Git 使用源代码。 要注意的是,Basemap使用了GEOS库,在编译安装Basemap时, 要将环境变量 GEOS_DIR 指向 libgeos_cgeos_c.h 的位置 (如果 libgeos_c/usr/local/lib 中, geos_c.h/usr/local/include 中, 则将 GEOS_DIR 设置为 /usr/local )。

安装时,请按照下列步骤进行操作:

解开Basemap版本X.Y.Z源tar.gz文件,和cd到底图-XY.Z目录。

安装GEOS库。如果你已经安装上了,只需设置。然后转到下一步。如果没有,可以通过下面步骤从包含在底图中的源代码构建它:

cd geos-3.3.3
export GEOS_DIR=<where you want the libs and headers to go>
./configure --prefix=$GEOS_DIR
make; make install

cd回到顶级Basemap目录(basemap-X.Y.Z)并运行安装 python setup.py 。 通过在python提示符下运行mpl_toolkits.basemap import basemap来检查并安装。

要测试,cd到examples目录并运行 python simpletest.py 。 如果运行所有的示例(除了那些有额外的依赖或需要互联网连接),请执行 python run_all.py

8.1.3. Basemap基本用法

>>> %matplotlib inline
>>> import warnings
>>> warnings.filterwarnings('ignore')
>>> from mpl_toolkits.basemap import Basemap
>>> import matplotlib.pyplot as plt

导入库中包括Basemap库和matplotlib。两者都是必要的:

>>> map = Basemap()

使用具有许多选项的Basemap类创建地图。 在没有传递任何选项的情况下,地图具有以经度=0和纬度= 0为中心的Plate Carrée投影。

在Mercator投影发明之前,航海中使用的是最简单的Plate Carrée投影(现在天地图还在用), 最早由Ptolemy发明,投影公式简单到不能再简单了:x=lon, y=lat。 但这个投影既不等角也不等积,特别在高纬度地区,与实际相差很大,所以并不实用。

>>> map.drawcoastlines()
>>> plt.show()
_images/basemap-basic_7_0.png

设置地图后,我们可以画出我们想要的。在这种情况下,海岸线图层,已经与库,使用方法drawcoastlines()

最后,地图必须显示已保存。使用mathplotlib的方法。在这个例子中,plt.show() 表示打开一个窗口来查看结果。 plt.savefig('file_name') 表示会将地图保存到图像文件中。这点以后不再赘述,书中还是以一般情况下还是使用 show() 方法来说明。

>>> plt.savefig('xx_test.png')
<Figure size 432x288 with 0 Axes>

在上面的使用中,我们没有声明地图的投影参数。 The default value is cyl, or Cylindrical Equidistant projection, also known as Equirectangular projection or Plate Carrée

更改投影很容易,只需将投影参数和lat_0和lon_0添加到“地图”中。 下面我们换另外一种投影 ortho ,看一下效果。 关于 Basemap 中投影的设置与使用,在后面会详细说明。 即使有新的投影,地图仍然有点空,所以让我们用一些颜色填充海洋和大陆,方法 fillcontinents()drawmapboundary() 将做到:

>>> map = Basemap(projection='ortho',lat_0=0, lon_0=105)
>>> map.drawmapboundary(fill_color='aqua')
>>> map.fillcontinents(color= 'coral',lake_color='aqua')
>>> map.drawcoastlines()
>>> map.drawcountries()
>>> plt.show()
_images/basemap-basic_13_0.png

用蓝色填充海洋,用陆地颜色填充大陆。整合在一起后,会发现海洋与大陆有着明显的区别。

>>> map = Basemap(projection='ortho',lat_0=0, lon_0=-15)
>>> map.drawmapboundary(fill_color='aqua')
>>> map.fillcontinents(color= 'coral',lake_color='aqua')
>>> map.drawcoastlines()
>>> map.drawcountries()
>>> plt.show()
_images/basemap-basic_15_0.png

8.1.4. 制作一个简单的地图

现在开始做一个简单的世界地图。运行下面的代码,你会得到一个漂亮的地图的地球,且有良好干净的海岸线:

>>> %matplotlib inline
>>> import warnings
>>> warnings.filterwarnings('ignore')
>>> import os
>>> from mpl_toolkits.basemap import Basemap
>>> import matplotlib.pyplot as plt
>>> import numpy as np

注意要确保分辨率的值是小写L,表示“低”,不是数字1。

>>> my_map = Basemap(resolution='i', projection='merc', llcrnrlat=10.0, urcrnrlat=55.0, llcrnrlon=60., urcrnrlon=140.0)

使用具有许多选项的Basemap类创建映射,并且传递选项数据。

>>> my_map.drawcoastlines()
>>> plt.show()
_images/basemap-basic_22_0.png

8.1.5. 添加详细信息

从国家边界开始,我们向这张地图添加一些细节。在 map.drawcoastlines() 后面添加以下行:

>>> my_map.drawcountries()
>>> my_map.fillcontinents(color='coral', lake_color='aqua')
>>> my_map.drawmapboundary(color='k', linewidth=1.0, fill_color='b', zorder=None, ax=None)
>>> plt.show()
_images/basemap-basic_24_0.png

数据中没有海洋, 则海洋作为背景处理。

在上图中你应该看到大陆被充满了。现在让我们来清理地球的边缘:

>>> my_map = Basemap(resolution='i', projection='merc', llcrnrlat=10.0, urcrnrlat=55.0, llcrnrlon=60., urcrnrlon=140.0)
>>> my_map.drawcoastlines()
>>> my_map.drawcountries()
>>> my_map.fillcontinents(color='coral')
>>> my_map.drawmapboundary()
>>>
>>> my_map.drawmeridians(np.arange(0, 360, 8))
>>> my_map.drawparallels(np.arange(-90, 90, 6))
>>>
>>> my_map.drawmapboundary(color='k', linewidth=1.0, fill_color='b', zorder=None, ax=None)
>>> plt.show()
_images/basemap-basic_27_0.png

绘制经线,0度开始360度结束,8个步长。纬线-90开始,90结束,6个步长。且海洋作为背景处理。

下面是另外一个投影的例子。与上个例子类似。步长为30.

>>> my_map = Basemap(projection='ortho', lat_0=0, lon_0=-100, resolution='l', area_thresh=1000.0)
>>> my_map.drawcoastlines()
>>> my_map.drawcountries()
>>> my_map.fillcontinents(color='coral')
>>> my_map.drawmapboundary()
>>> my_map.drawmeridians(np.arange(0, 360, 30))
>>> my_map.drawparallels(np.arange(-90, 90, 30))
>>> my_map.drawmapboundary(color='k', linewidth=1.0, fill_color='b', zorder=None, ax=None)
>>> plt.show()
_images/basemap-basic_30_0.png

8.1.6. 设置地图投影

为了在二维地图上表示地球的曲面,则需要进行地图投影。地图投影的方法有许多种,每种方法都有自己的优点和缺点。 Basemap提供了24种地图投影方法。 有些是全球性的,有些只能代表区域。 在创建Basemap类实例时,必须指定所需的地图投影和地图投影将描述的地球表面部分的信息。 有两种基本的方法。一个是提供矩形映射投影区域的四个角的每一个的纬度和经度值。 另一个是提供地图投影区域中心的 lat/lon 值以及地图投影坐标中的区域的宽度和高度。 类变量 supported_projections 是一个字典(非 Python 数据结构中的字典,自然语言中的字典之意,实际上是 Python 的字符串), 包含有关Basemap支持的所有投影的信息。

>>> %matplotlib inline
>>> import warnings
>>> warnings.filterwarnings('ignore')
>>>
>>> import mpl_toolkits.basemap
>>>
>>> print(mpl_toolkits.basemap.supported_projections)
cyl              Cylindrical Equidistant
merc             Mercator
tmerc            Transverse Mercator
omerc            Oblique Mercator
mill             Miller Cylindrical
gall             Gall Stereographic Cylindrical
cea              Cylindrical Equal Area
lcc              Lambert Conformal
laea             Lambert Azimuthal Equal Area
nplaea           North-Polar Lambert Azimuthal
splaea           South-Polar Lambert Azimuthal
eqdc             Equidistant Conic
aeqd             Azimuthal Equidistant
npaeqd           North-Polar Azimuthal Equidistant
spaeqd           South-Polar Azimuthal Equidistant
aea              Albers Equal Area
stere            Stereographic
npstere          North-Polar Stereographic
spstere          South-Polar Stereographic
cass             Cassini-Soldner
poly             Polyconic
ortho            Orthographic
geos             Geostationary
nsper            Near-Sided Perspective
sinu             Sinusoidal
moll             Mollweide
hammer           Hammer
robin            Robinson
kav7             Kavrayskiy VII
eck4             Eckert IV
vandg            van der Grinten
mbtfpq           McBryde-Thomas Flat-Polar Quartic
gnom             Gnomonic
rotpole          Rotated Pole

键是短名称(与在创建Basemap类实例时用于定义投影的projection关键字一起使用),对应的后面是长的、更具描述性的名称。 类变量projection_params是一个字典,提供可用于定义每个投影的属性的参数列表。 以下是说明如何设置每个支持的投影的示例。 注意,许多地图投影具有两个期望的属性之一 - 它们可以是等面积(保留特征的面积)或保形(保留特征的形状)。 因为没有地图投影可以同时具有两者,所以在使用的时候选择哪种投影要在两者之间的做许多妥协。

所有地图必须进行地图投影。投影及其特征都是在创建对象 Basemap 时要明确的。 这样做与其他库(如 GDAL)有很大的不同,理解这一点对于使用Basemap是非常重要的。

>>> from mpl_toolkits.basemap import Basemap
>>>
>>> import matplotlib.pyplot as plt
>>> import numpy as np
>>> map = Basemap(projection='cyl')
>>>
>>> map.drawcoastlines()
>>> map.drawmeridians(np.arange(0, 360, 30))
>>> map.drawparallels(np.arange(-90, 90, 30))
>>> plt.show()
_images/basemap-basic_34_0.png

projection 参数是设置地图投影时要使用的。 默认值为cyl或圆柱等距投影,等距投影也称为等方矩形投影或板卡雷。

许多投影需要额外的参数:

>>> map = Basemap(projection='aeqd', lon_0=180, lat_0=50)
>>> map.drawmapboundary(fill_color='aqua')
>>> map.fillcontinents(color='coral', lake_color='aqua')
>>> map.drawcoastlines()
>>> plt.show()
_images/basemap-basic_36_0.png

该地图在欧洲区域已经具有以经度= 10和纬度= 50为中心的等距投影。一些投影需要更多参数,在手册的每个投影页中都有描述。 Basemap对象具有字段proj4string,该字段具有要与proj4一起使用的字符串,用于计算投影参数。

使用EPSG设置投影

EPSG代码是使用数字代码来命名投影的标准的。 Basemap允许使用此表示法来创建地图,但仅在某些情况下。 若要使用它,需要将epsg参数传递给带有代码的Basemap构造函数。 Basemap支持的epsg代码在文件 /mpl_toolkit /basemap/data/epsg 中。 即使所需的epsg出现在文件中,但有时库并不能使用投影,例如:

ValueError:23031不是受支持的EPSG代码

这里并没有很好地支持名称为“utm”的投影(即23031或15831),但可以使用名为tmerc的投影。首先要做的是打开文件,然后寻找一个合适的选项。

>>> map = Basemap(projection='mbtfpq', lon_0=105)

mbtfpq 投影。

>>> map.drawcoastlines()
>>> map.drawmeridians(np.arange(0, 360, 30))
>>> map.drawparallels(np.arange(-90, 90, 30))
>>> plt.show()
_images/basemap-basic_40_0.png

范围

>>> my_map = Basemap(projection='ortho', lat_0=0, lon_0=-100,
>>>                  resolution='l', area_thresh=1000.0)
>>> my_map.drawcoastlines()
>>> my_map.drawcountries()
>>> my_map.fillcontinents(color='coral')
>>> my_map.drawmapboundary()
>>> my_map.drawmeridians(np.arange(0, 360, 30))
>>> my_map.drawparallels(np.arange(-90, 90, 30))
>>> plt.show()
_images/basemap-basic_42_0.png

根据所有的例子,直到现在采取整个地球。仅绘制区域可以通过边界框或地图的中心和地图大小来完成。官方文档说,这两种方法都可以在大多数时候使用,但有很多例外。

注:使用cyl,merc,mill,cea和gall投影,如果没有设置,角默认假设为-180,-90,180,90(全球)。

传递边界框

>>> map = Basemap(projection='cyl')
>>> map.drawcoastlines()
>>> map.drawmeridians(np.arange(0, 360, 30))
>>> map.drawparallels(np.arange(-90, 90, 30))
>>> plt.show()
_images/basemap-basic_44_0.png

左下角和右上角以经度 - 纬度为单位作为参数而非地图单位原来参数。这就是为什么一些投影失败的原因,因为经度 - 纬度的正方形可能不会在投影单位中给出良好的边界框。 在该示例中,使用的是UTM(横向墨卡托投影)。在这种情况下,边界框方法更容易,因为从地图中心计算UTM单位的宽度要困难得多。 注:使用sinu,moll,hammer,npstere,spstere,nplaea,splaea,npaeqd,spaeqd,robin,eck4,kav7或mbtfpq投影时,则此方法不能使用。

在地图坐标中传递边界框

>>> map = Basemap(projection='mbtfpq', lon_0=105)
>>> map.drawcoastlines()
>>> map.drawmeridians(np.arange(0, 360, 30))
>>> map.drawparallels(np.arange(-90, 90, 30))
>>> plt.show()
_images/basemap-basic_47_0.png

有些投影(看起来像卫星图像的投影)接受使用地图坐标设置扩展。首先必须设置投影参数(中心点),然后要显示的区域只能是地球的一部分。

注:只有ortho,geos和nsper投影可以使用此方法设置地图扩展。

通过中心,宽度和高度

>>> map = Basemap(projection='sinu', lon_0=105, lat_0=39)
>>> map.drawcoastlines()
>>> map.drawmeridians(np.arange(0, 360, 30))
>>> map.drawparallels(np.arange(-90, 90, 30))
>>> plt.show()
_images/basemap-basic_49_0.png

在这种情况下,投影的中心,投影的宽度和高度作为参数传递。找到投影中心很容易,只是通过它在经度 - 纬度。但是大小有点棘手: 单位是以米为单位的投影单位,点(0,0)是左下角,点(宽度,高度)是右上角。因此,位置的原点不是由GDAL中的投影定义的。投影只是定义所使用的单位的大小,而不是原点。 该示例使用绘图函数来显示几个点的位置,以显示坐标从0到宽度和高度的范围。