>>> from env_helper import info; info()
页面更新时间: 2022-05-27 12:17:06
运行环境:
    Linux发行版本: Debian GNU/Linux bookworm/sid
    操作系统内核: Linux-5.16.18-200.fc35.x86_64-x86_64-with-glibc2.33
    Python版本: 3.10.4

17.2. 开始使用 Basemap

现在开始在 Python 中使用 Basemap 。

17.2.1. 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/sec2_basic_8_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/sec2_basic_14_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/sec2_basic_16_0.png

17.2.2. 制作一个简单的地图

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

>>> %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/sec2_basic_23_0.png

17.2.3. 添加详细信息

从国家边界开始,我们向这张地图添加一些细节。在 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/sec2_basic_25_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/sec2_basic_28_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/sec2_basic_31_0.png

17.2.4. 设置地图投影

为了在二维地图上表示地球的曲面,则需要进行地图投影。地图投影的方法有许多种,每种方法都有自己的优点和缺点。 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/sec2_basic_35_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/sec2_basic_37_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/sec2_basic_41_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/sec2_basic_43_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/sec2_basic_45_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/sec2_basic_48_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/sec2_basic_50_0.png

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