制作一个简单的地图

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

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

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

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

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

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

添加详细信息

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

In [13]:
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 [14]:
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个步长。且海洋作为背景处理。

In [ ]:
下面是另外一个投影的例子与上个例子类似步长为30.
In [15]:
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 [16]:
%matplotlib inline
import warnings
warnings.filterwarnings('ignore')

import mpl_toolkits.basemap

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

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

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

In [9]:
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 [10]:
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 [3]:
map = Basemap(projection='mbtfpq', lon_0=105)

用‘mbtfpq’投影。

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

范围

In [5]:
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 [6]:
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 [7]:
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 [8]:
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 [1]:
%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 [21]:
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 [22]:
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 [23]:
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 [24]:
# 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 [25]:
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()

没有绘制海岸线,结果有些奇怪:

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

6、Drawlsmask

Drawlsmask是绘制湖泊,陆地和海洋的方法。尽量避免使用fillcontinents和drawmapboundary方法。此外,Drawlsmask可以将数据源更改为自定义点阵列。 与其他方法的不同之处在于,在此方法中无法设置zorder。

drawlsmask(land_color=‘0.8’, ocean_color=’w’, lsmask=None, lsmask_lons=None, lsmask_lats=None, lakes=True, resolution=’l’, grid=5, **kwargs)

  • land_color设置给土地的颜色(默认为灰色)
  • ocean_color设置海洋的颜色(默认为白色)
  • lsmask具有备用数据的阵列。如果为“无”,则将采用库中的默认数据。该数组必须包含海洋的0和土地的1
  • lsmask_lons替换陆地海面掩码的经度
  • lsmask_lats替代陆地海面的纬度
  • 分辨率可以更改由Basemap实例定义的数据的分辨率
  • grid栅格数组栅格分辨率,单位为弧度。默认值为5分钟
In [27]:
# from mpl_toolkits.basemap import Basemap
# import matplotlib.pyplot as plt

# map = Basemap(projection='ortho', 
#               lat_0=0, lon_0=0)

# map.drawlsmask(land_color = "#ddaa66", 
#                ocean_color="#7777ff",
#                resolution = 'l')

# plt.show()

7、drawmapboundary

在地图上绘制地球边界,可选择填充。

drawmapboundary(color=’k’, linewidth=1.0,     fill_color=None, zorder=None, ax=None)

  • 线宽集。线宽以像素为单位
  • 颜色设置边缘颜色,默认为k(black)。遵循matplotlib约定
  • fill_color设置填充球体颜色,默认情况下为None。还遵循matplotlib约定
  • zorder设置图层位置。默认情况下,顺序由Basemap设置
In [28]:
plt.figure(121)
map = Basemap(projection='ortho',lon_0=0,lat_0=0,resolution='c')
map.drawmapboundary()

plt.figure(122)
map = Basemap(projection='sinu',lon_0=0,resolution='c')
map.drawmapboundary(fill_color='aqua')

plt.show()

8、Drawmeridians

在地图上绘制经线。

drawmeridians(meridians, color=’k’, linewidth=1.0, zorder=None, dashes=[1, 1], labels=[0, 0, 0, 0], labelstyle=None, fmt=’%g’, xoffset=None, yoffset=None, ax=None, latmax=None, **kwargs)

  • 经线是一个有经度的列表。如果值是整数,则可以使用range()创建。如果你需要浮点数,np.arange()是一个不错的选择
  • color设置线的颜色。
  • 线宽集。线宽以像素为单位
  • zorder可以设置线的位置,土地覆盖经线,或相反。
  • dashing可以设置dash样式。第一个元素是要绘制的像素数,第二个是要跳过的像素数
  • 更改绘制标签的位置。将值设置为1使得标签在地图的所选边上绘制。四个位置是[左,右,上,下]
In [29]:
map = Basemap(projection='aeqd', 
              lon_0=0.0, lat_0=0,             
              width=25000000, height=25000000)

map.drawmeridians(range(0, 360, 20))

plt.show()

9、parallels

在地图上绘制纬线。

drawparallels(circles, color=’k’, linewidth=1.0, zorder=None, dashes=[1, 1], labels=[0, 0, 0, 0], labelstyle=None, fmt=’%g’, xoffset=None, yoffset=None, ax=None, latmax=None, **kwargs)

  • rallels是要绘制的纬度的列表。如果值是整数,则可以使用range()创建。如果你需要浮点数,np.arange()是一个不错的选择。
  • color设置线的颜色。
  • 线宽集,当然,线宽以像素为单位
  • order可以设置线的位置,以便能够例如使得平面