简介

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 模块。

安装

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

Basemap基本用法

In [60]:
%matplotlib inline 
import warnings
warnings.filterwarnings('ignore')    
from mpl_toolkits.basemap import Basemap
import matplotlib.pyplot as plt

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

In [61]:
map = Basemap()

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

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

In [62]:
map.drawcoastlines()
plt.show()

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

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

In [63]:
plt.savefig('xx_test.png')
<matplotlib.figure.Figure at 0x7f022fb331d0>

在上面的使用中,我们没有声明地图的投影参数。 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() 将做到:

In [64]:
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()

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

In [65]:
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()

制作一个简单的地图

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

In [66]:
%matplotlib inline
import warnings
warnings.filterwarnings('ignore')
In [67]:
import os
from mpl_toolkits.basemap import Basemap
import matplotlib.pyplot as plt
import numpy as np

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

In [68]:
my_map = Basemap(resolution='i', projection='merc', llcrnrlat=10.0, urcrnrlat=55.0, llcrnrlon=60., urcrnrlon=140.0)

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

In [69]:
my_map.drawcoastlines()
plt.show()

添加详细信息

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

In [70]:
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()

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

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

In [71]:
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()

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

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

In [72]:
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()

设置地图投影

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

In [73]:
%matplotlib inline
import warnings
warnings.filterwarnings('ignore')

import mpl_toolkits.basemap

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

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

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

In [74]:
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()

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

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

In [75]:
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()

该地图在欧洲区域已经具有以经度= 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的投影。首先要做的是打开文件,然后寻找一个合适的选项。

In [76]:
map = Basemap(projection='mbtfpq', lon_0=105)

用‘mbtfpq’投影。

In [77]:
map.drawcoastlines()
map.drawmeridians(np.arange(0, 360, 30))
map.drawparallels(np.arange(-90, 90, 30))
plt.show()

范围

In [78]:
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()

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

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

传递边界框

In [79]:
map = Basemap(projection='cyl')
map.drawcoastlines()
map.drawmeridians(np.arange(0, 360, 30))
map.drawparallels(np.arange(-90, 90, 30))
plt.show()

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

在地图坐标中传递边界框

In [80]:
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()

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

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

通过中心,宽度和高度

In [81]:
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()

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

绘制地图背景

Basemap包括GSSH海岸线数据集,以及来自GMT的河流、州和国家边界的数据集。这些数据集可用于在地图上以几种不同的分辨率绘制海岸线,河流和政治边界。相关的底图方法是:

  • drawcoastlines():绘制海岸线。
  • fillcontinents():为大陆内部着色(通过填充海岸线多边形)。不幸的是,fillcontinents方法并不总是做正确的。因为Matplotlib总是尝试填充多边形的内部。在某些情况下,海岸线多边形的内部可以是模糊的,并且有时需要填充外部而不是内部。在这些情况下,建议的解决方法是使用drawlsmask()方法覆盖具有为陆地和水域指定的不同颜色的图像。
  • drawcountries():绘制国家边界。
  • drawstates():绘制北美的州界。
  • drawrivers():绘制河流。

绘制地图并不只包括绘制海岸线和政治边界,也包括地图背景的绘制。Basemap提供了绘制地图背景的几个方法:

  • drawlsmask():绘制高分辨率海陆掩码作为图像,指定土地和海洋颜色。陆地海面掩模源自GSHHS海岸线数据,并且有多个海岸线选项和像素大小可供选择。
  • bluemarble():绘制一张NASA蓝色大理石图像作为地图背景。
  • shadedrelief():绘制一个阴影浮雕图像作为地图背景。
  • etopo():绘制etopo浮雕图像作为地图背景。
  • warpimage():使用abitrary图像作为地图背景。图像必须是全球的,从国际数据线向东,南极向北,以纬度/经度坐标覆盖世界。

1、使用arcgis REST API服务下载并绘制图像。

In [82]:
%matplotlib inline
import warnings
warnings.filterwarnings('ignore')

from mpl_toolkits.basemap import Basemap
import matplotlib.pyplot as plt

map = Basemap(llcrnrlon=3.75,llcrnrlat=39.75,urcrnrlon=4.35,urcrnrlat=40.15, epsg=5520)
#http://server.arcgisonline.com/arcgis/rest/services

map.arcgisimage(service='ESRI_Imagery_World_2D', xpixels = 1500, verbose= True)
plt.show()
http://server.arcgisonline.com/ArcGIS/rest/services/ESRI_Imagery_World_2D/MapServer/export?bbox=1564270.9620172717,4401598.726056214,1615017.0560380095,4446612.184939825&bboxSR=5520&imageSR=5520&size=1500,1330&dpi=96&format=png32&f=image

下面使用了另外一种底图:

In [83]:
map.arcgisimage(service='World_Shaded_Relief', xpixels = 1500, verbose= True)
plt.show()
http://server.arcgisonline.com/ArcGIS/rest/services/World_Shaded_Relief/MapServer/export?bbox=1564270.9620172717,4401598.726056214,1615017.0560380095,4446612.184939825&bboxSR=5520&imageSR=5520&size=1500,1330&dpi=96&format=png32&f=image
arcgisimage(server='http://server.arcgisonline.com/ArcGIS', service='ESRI_Imagery_World_2D',
    xpixels=400, ypixels=None, dpi=96, verbose=False, **kwargs)
  • 服务器可以使用相同的REST API连接到另一个服务器。
  • 服务器显示的是绘制的层。要获取可用图层的完整列表,请检查API网。
  • xpixels实际上是用来设置图像的缩放。数字越大将会要求图像更大,所以图像将有更多的细节。因此,当缩放更大时,xsize必须更大以保持分辨率 。
  • y像素可以用于强制图像在y方向上具有与以宽高比定义的像素不同数量的像素。默认情况下,宽高比保持不变,这似乎是一个很好的价值。
  • dpi是输出的图像分辨率。更改其值将更改像素数,但不会更改缩放级别。
  • verbose打印用于获取远程图像的URL。
  • 使用这种方法的一个重要的点是,投影必须使用epsg参数设置,除非使用底图表示法中的4326或cyl。

2、bluemarble

在地图上绘制bluemarble图像。

bluemarble(ax=None, scale=None, **kwargs)

  • 该比例对于降低原始图像分辨率以加速处理是有用的。将图像的大小除以4值为0.5。
  • 图像的最终投影被扭曲,所以所有投影方法均正常工作。
In [84]:
map = Basemap(llcrnrlon=-10.5,llcrnrlat=33,urcrnrlon=10.,urcrnrlat=46.,
             resolution='i', projection='cass', lat_0 = 39.5, lon_0 = 0.)

map.bluemarble()

map.drawcoastlines()

plt.show()

3、drawcoastlines

绘制海岸线。

drawcoastlines(linewidth=1.0, linestyle=’solid’, color=’k’, antialiased=1, ax=None, zorder=None)

  • 线宽集。线宽以像素为单位。
  • linestyle设置线型。默认情况下是固体,但可以是虚线,或任何matplotlib选项。
  • 颜色默认为k(black)。遵循matplotlib约定。
  • 抗锯齿默认为true
  • zorder设置图层位置。默认情况下,顺序由Basemap设置。
In [85]:
map = Basemap()

map.drawcoastlines()

plt.show()

4、drawcounties

从库包括的图层绘制国家边界。

drawcounties(linewidth=0.1, linestyle=’solid’, color=’k’, antialiased=1, facecolor=’none’, ax=None, zorder=None, drawbounds=False)

  • 线宽集。线宽以像素为单位。
  • linestyle设置线型。默认情况下是固体,但可以是虚线,或任何matplotlib其他选项。
  • 颜色默认为k(black)。遵循matplotlib约定。
  • 抗锯齿默认为true。
  • zorder设置图层位置。默认情况下,顺序由Basemap设置。

注:facecolor参数,它应该是县的颜色,但在一些Basemap版本不能运行。

注意:

  • 分辨率是固定的,不依赖于传递给类构造函数的分辨率。
  • 海岸线是另一个功能,而乡镇海岸线不被认为是海岸线,因此有必要把这个方法和其他的方法结合起来以获得好的地图。
In [86]:
# map = Basemap(llcrnrlon=-93.,llcrnrlat=40.,urcrnrlon=-75.,urcrnrlat=50.,
#              resolution='i', projection='tmerc', lat_0 = 40., lon_0 = -80)

# map.drawmapboundary(fill_color='aqua')
# map.fillcontinents(color='#cc9955', lake_color='aqua')

# map.drawcounties()

# plt.show()

5、drawcountries

从库中包含的图层绘制国家/地区边界。

drawcountries(linewidth=1.0, linestyle=’solid’, color=’k’, antialiased=1, ax=None, zorder=None)

  • 线宽集。线宽以像素为单位。
  • linestyle设置线型。默认情况下是固体,但可以是虚线,或任何matplotlib选项。
  • 颜色默认为k(black)。遵循matplotlib约定。
  • 抗锯齿默认为true zorder设置图层位置。默认情况下,顺序由Basemap设置。

注意:

  • 创建Basemap实例时会产生具有更好或粗糙分辨率的图层。
  • 海岸线是另一个功能,乡镇海岸线不是通常的海岸线,因此要结合其它函数来绘制地图。
In [87]:
map = Basemap(projection='ortho', 
              lat_0=0, lon_0=0)

map.drawmapboundary(fill_color='aqua')
map.fillcontinents(color='coral',lake_color='aqua')

map.drawcountries()

plt.show()