安装basemap模块: https://zhuanlan.zhihu.com/p/3450984
1 前言
地理信息可视化是我们测绘、气象、地球科学方向的学子常会面对的任务。一副好图胜过大段的文字,可以直观揭示事物的特点,帮助我们探索问题的本质。高质量的成图也可以提高我们科研论文和工作报告的表现力,提高文章的档次。
表1中对笔者和周围朋友常用的地球科学绘图软件做了简单的汇总和对比。一般来说,Excel成图美观度较差,难以通过脚本完成自动化操作,主要用于快速的数据分析。因此,用于地理信息可视化的软件主要是GMT、Matplotlib的basemap模块,和Matlab的M_map模块。
GMT和M_map笔者早先使用过,目前笔者使用的是Python的basemap模块。除了功能强大,语法简单、清晰,轻量级(安装包占用空间不大)等优势之外,另一个主要原因在于Python强大的生态和活跃的社区。无论是数据处理、分析,还是科学计算,日常工作自动化等方面,Python都有丰富且热门的模块。笔者使用Python作为主要开发语言之一,因此就选择了basemap模块,将其集成在数据处理和分析任务中。
笔者周围不少小伙伴对basemap画图很感兴趣。本文根据小伙伴们对basemap模块关心的一些问题,收集资料,结合笔者的一些初步经验,对basemap画图做了初步总结,供大家参考。
文章全文约6000字,包括15个代码(片段),21幅样例图。受篇幅所限,部分较长的代码并未在文中全部贴出。如果您有兴趣,可在后台留言"Python-basemap"和您的邮箱地址,笔者会将本期PDF文稿和相应数据、脚本发送给您。
2 简介
在数据可视化过程中,我们需要将数据在地图上画出来。 比如说我们在地图上画出城市人口,飞机航线,军事基地,矿藏分布等等。这样的地理绘图有助于读者理解空间相关的信息。basemap 是Python的一个强大的负责实现地理信息可视化的库,是Matplotlib的一个附加工具包,通过结合 matplotlib 可以绘制出很多漂亮的地图。
basemap包括GSSH海岸线数据集,以及来自GMT的河流、州和国家边界的数据集。这些数据集可用于在地图上以几种不同的分辨率绘制海岸线,河流和政治边界。basemap底层使用了Geometry Engine-Open Source(GEOS)库,用来将海岸线和边界特征剪切到所需的地图投影区域。此外,basemap还提供读取shapefile的功能。
basemap面向地球科学家的需求,特别是海洋学家和气象学家。起初,杰夫·惠特克(Jeff Whitaker)写“底图”(basemap)用来帮助他进行气候和天气预报的研究。到现在,basemap 不断被开发扩展已经具备了很多功能。多年来,basemap的功能随着各个学科(如生物学,地质学和地球物理学)的科学家的要求和贡献的新功能而演变。
3 关于安装
问1:basemap模块的大小和安装的难易程度?
答:basemap模块约为120M左右,其依赖库pyproj约为20M左右。
由于basemap还没有pip检索,因此传统的(pip install basemap)安装方法不适用。
Windows用户需要到https://www.lfd.uci.edu/~gohlke/pythonlibs/下载对应的wheel文件到本地,然后控制台进入其所在目录,使用pip install xxxx.whl文件安装。其依赖于pyproj库。具体安装过程参考https://zhuanlan.zhihu.com/p/34509847 。Linux用户安装步骤可参考 https://matplotlib.org/basemap/users/installing.html ,包括下载源码、安装GEOS库、设置GEOS_DIR环境变量等步骤,安装完成后可进入到examples文件夹,运行测试例子。 此外还可以通过Anaconda的Python发行版安装。
建议同时安装Pillow库, 可用于以更多不同格式保存图片文件,以及调用bluemarble()
, etopo()
, shadedrelief()
and warpimage()
函数。
问2:basemap使用是否方便?有没有相关介绍或例子?
答:basemap基于GEOS的地图二维数据,其底图数据库与GMT相同,封装了大量常用的地图投影、坐标转换功能,利用简洁的Python语法支持绘出多种多样的地理地图。
如下所示,basemap的功能有且不限于绘制地图、等高线、降水量图、海平面气压天气图、晨昏线等,具体功能及对应实例代码可参考https://matplotlib.org/basemap/users/examples.html。本文后续内容也给出了多个样例。
4 基本操作
问3:能否在地图上画线、画箭头之类的标注线?
答:可以。比如需要连接地图上已知经纬度(lat1,lon1;lat2,lon2)的两个点,可使用plot来进行连线,basemap中暂没有画箭头的函数,可用plt.arrow绘制箭头。一些更高级的用法,比如画球面距离,可使用drawgreatcircle函数。下图实例为绘制伦敦到纽约的直线以及球面距离。
问4:如何导出矢量图?
答:使用函数savefig(),保存矢量图只需要三个参数,即文件名fname,DPI和文件格式format。比如保存文件名为“vectorgraph”,dpi为600,eps格式的图形,可使用如下代码。目前savefig支持的格式包括:eps,jpeg,jpg,pdf,pgf,png,ps,raw,rgba,svg,svgz,tif,tiff。值得注意的是,如果文件名含后缀,则format参数可省略。
问5:经纬度怎么更简洁的表示,全部显示有点多和乱,国际上有惯例吗?
答:经纬度的表示没有国际惯例,需要用户自己设置。drawparallels和drawmeridians函数可通过设置经纬度间隔来显示图上的经纬度。且能分别设置图片左右边框和上下边框是否显示经纬度。
问6:地图怎么画比例?
答:画比例尺的函数为drawmapscale。下图给出了两种比例尺示例。
问7:能否调节地图分辨率?
答:可以在创建 basemap 图像时设置所需的分辨率。basemap
类的resolution
参数设置分辨率级别,他们是'c'
(原始),'l'
(低),'i'
(中),'h'
(高),'f'
(完整)或None
(如果没有使用边界)。这个选项很重要,在全局地图上设置高分辨率边界可能非常慢。
问8:如何设置投影方式?
答:basemap的投影方式是由属性projection控制,提供了包括:Robinson 罗宾逊投影、Orthographic 正射投影、Mercator 魔卡托投影、Sinusoidal 正弦投影、Lambert 兰勃特投影在内的近30种投影方式。默认投影为‘cyl’(圆柱投影)。具体的投影方式对应的设置方法可参考https://matplotlib.org/basemap
问9:地图能否设置长宽比?
答:basemap地图的长宽比是固定的。例如,robin投影为椭圆,不能通过设置进行拉伸;cyl投影为矩形,长宽为固定比例,也不能改变。
问10:basemap能否显示省级信息?
答:全球的行政区划数据都可以在GADM下载(https://gadm.org/download_country_v3.html)。
如果想在basemap中显示中国区域的省级信息,需到上述网址下载中国区域的shapefile文件,解压之后的文件名中会有0,1,2,3后缀,分别对应于第零、一、二、三级行政区划界限,基本相当于中国大陆地区的边界、省界、市界和区县界。可通过readshapefile函数读取相关信息。(PS:数据中有领土纠纷的部分不代表笔者立场)
问11:假如数据是离散的,是怎么样的插值让图看起来不是单个网格?
答:如果有固定的数学模型,可根据数学模型加密数据,使图形更平滑。如果没有确定的数学模型,python中提供了一种样条插值办法,可使用scipy.interpolate模块下的interpld函数实现插值,加密数据。
问12:basemap画图有办法用滚轮进行缩放吗?
答:可以。缩放等鼠标交互绘制操作不属于basemap库,可使用fig.canvas.mpl_connect()函数来绑定相关fig的交互事件,此类交互事件包括:鼠标按下与释放,按键按下与释放,滚轮滚动等,其中滚轮操作属性为“button”,滚动属性为“up”和“down”。具体可参考https://www.jianshu.com/p/cf205a759
问13:能否绘制多个子图在一个图中,比如卫星弧段图,SNR图等?
答:Python可通过subplot函数画多个子图。调用形式如:subplot(nrows,ncols,index),图表的整个绘图区域被分成nrows行和ncols列,按照从左到右,从上到下的顺序对每个子区域进行编号,左上的子区域编号为1。index参数指定创建的Axes对象所在的区域。
matplotlib官网上的一个2*2多子图例子如下所示:
比如要将卫星弧段图和SNR图两个图按照上下顺序放置于同一张图时,可采用如下方式:
5 高级操作
问14.1:basemap可以实现画全球和局部测站分布图吗?
问14.2:测站分布图可以按照可观测卫星系统类型进行站点分类吗?
答:可以实现。对于全球或局部测站分布图,可通过控制经纬度显示范围得到。对于想表示的要素,比如卫星系统或者DOP值,需要将测站能观测到的卫星系统或者DOP值作为参数读取到代码中,用不同颜色对参数的值进行表示即可。以下将以按照可观测卫星系统的值给出全球测站分布图。标注测站ID,需要在代码中进行微调,将重叠部分的ID上下左右移动。
问15:basemap画全球气温分布图,该怎么在底图上叠加数据?
答:在绘制气温分布图时,需要用到读写netcdf文件的第三方库:netCDF4。对于标准的气温数据,可在NCEP/NCAR上下载.nc文件,调用netCDF4库读取此文件绘制气温分布图。下图用月平均气温为例,绘制了全球范围内某月月平均气温图。
问16.1:basemap能不能画全球范围可见卫星数图(用渐变色表示)?
问16.2:怎么画电离层热力图?
答:basemap可以画可见卫星数图和电离层热力图,两者是同类型。以前者为例,需要首先自己划分格网点,然后分别计算各格网点的单天平均可见卫星数,形成格网文件,然后画图。下图以COD20342.sp3文件中给出的GPS、Galileo、GLONASS、BDS和QZSS的卫星轨道参数为例,绘制当天全球范围内单天平均可见卫星数图。为了简化代码量,有的函数只注释了功能,未给出具体代码。
问17:能否把测站属性用渐变颜色表述出来?比如希望以DOP值为颜色变量,快速筛选站星分布较好的测站。
答:可以实现。首先根据测站计算当天平均DOP值,然后读取该文件画图。由于全球测站较多,现以北美地区为例,用不同颜色区分DOP大小,方便用户根据DOP值筛选测站。
问18:怎么画卫星星下点轨迹图?
答:可首先根据SP3文件计算当天各时刻各卫星的位置,转换为经纬度后,调用basemap的scatter函数进行绘制。下图将以某一天的SP3文件为例,绘制当天C06,G01和J01三种卫星的星下点轨迹图。由于代码量较多,只给出绘图部分的源码。
问19:能否绘制卫星分布图(skyplot)?
答:可以。绘制skyplot时不用调用basemap模块。首先计算出各历元各卫星高度角和方位角,然后当卫星可见时,使用极坐标即可绘制卫星分布图(雷达图)。下图给出在2019年1月1日测站abpo上空C13卫星的运行轨迹图。
问20:能否像GMT一样在测站分布图外围加黑白框?
答:这是笔者目前用basemap唯一无法达到和GMT一样效果的地方。想要加上可能需要额外写函数实现。
问21.1:背景图能否显示出像谷歌地图或谷歌卫星云图效果?
问21.2:能否在地图上显示出绿化/沙漠信息或者高程信息等三维信息?
答:在basemap模块中,可通过不同函数,实现不同地图背景地图的更换。如,使用arcgisimage可显示类似谷歌卫星云图的效果,使用etopo函数可显示类似浮雕图的高程信息等。下图以类似谷歌地图的显示效果为例给出示例图:
etopo函数绘制由NOAA提供的浮雕地形图。
bluemarble函数绘制蓝大理石背景地形图。
shadedrelief函数绘制来源于http://www-shadedrelief.com网站的阴影浮雕背景图。
更多用法可参考:https://basemaptutorial.readthedocs.io/en/latest/backgrounds.html#
问22:basemap能否画动图?比如很多测站不同时间的数据,随时间变化闪烁的动图?
答:可以。最简单的方式为每次画之前清除画布,然后将之前的数据叠加之后再一起画出,即可画出动图,下面代码以之前的skyplot为例,以动图的方式画出C13卫星一天的运行轨迹,感兴趣的读者可自行运行。
6 与其他工具对比
问23:basemap画图相比Matlab有哪些优点?学习起来哪个更容易?
答:Python和Matlab的语法都很简单易懂,是很容易学习的编程语言。
与M_map相比,basemap和Python的优势在于安装包占用空间小,免费开源。相比动不动占用几十G空间,还需购买正版license的Matlab相比,前者安装下来占用空间不到1G。
问24:basemap相比于GMT优势在哪?
答:优势在于Python的语法更简单易懂。
画海岸线就是drawcoastlines()
函数,填充大陆颜色就是fillcontinents()
函数。很容易理解。
basemap画图时后面的函数功能直接加到图上。GMT4和5版本绘制一张稍复杂的图时,经常需要在原有的代码中增添、删除或修改已有命令的顺序,尤其需要注意 -K
、 -O
以及重定向符号的使用。用户不注意很容易画图失败。这也是basemap较之GMT更方便的地方。
问25:对那些有固定GMT代码每次只需要修改数据的图(比如测站分布图),basemap有什么优点?
答:basemap在相同场景下同样能固定脚本,只需要修改数据绘制新图。
因为Python语言的简单,basemap的一个优点在于一段时间后再来阅读和修改代码都相当容易,不会出现看不懂的情况,可维护性高。
问26.1:如果标注测站ID的话,ID会覆盖站点,密集区域会重叠,需要怎么调整?
问26.2:GMT画全球测站分布图时,为了不让测站名重叠,需要人工微调,basemap有没有更好的解决办法?
答:据笔者所知,为使测站名不重叠,两个工具都需要人工微调测站名的位置,还没有智能到能自己调节的地步。
7 其他
问27:basemap主要参考网站和资料有哪些?
答:首先推荐basemap官网:
[1] https://basemaptutorial.readthedocs.io/en/latest/index.html
[2] https://matplotlib.org/basemap/index.html
其次就是本文小结:-)
还有basemap源码。
附:basemap基本函数
basemap 包包含一系列有用的函数,用于绘制物理特征的边界,例如大陆,海洋,湖泊和河流等,以及政治边界,例如国家地区。以下是一些可用绘图功能:
- 物理边界和水体
-
drawcoastlines()
:绘制大陆海岸线 -
drawlsmask()
:绘制高分辨率海陆掩码作为图像,指定土地和海洋颜色。陆地海面掩模源自GSHHS海岸线数据,并且有多个海岸线选项和像素大小可供选择 -
drawmapboundary()
:绘制地图边界,包括海洋的填充颜色 -
drawrivers()
:在地图上绘制河流 -
fillcontinents()
:用给定的颜色填充大陆;可选择用另一种颜色填充湖泊 - 政治边界
-
drawcountries()
:绘制国界 -
drawstates()
:绘制北美的州界 -
drawcounties()
:绘制美国县界 - 地图功能
-
drawgreatcircle()
:在两点之间绘制大圆圈 -
drawparallels()
:绘制指定纬度的线条 -
drawmeridians()
:绘制指定经度的线条 -
drawmapscale()
:在地图上绘制比例尺刻度 - 全球图像
-
bluemarble()
:将 NASA 的蓝色大理石图像作为地图背景 -
shadedrelief()
:将阴影浮雕图像作为地图背景 -
etopo()
:绘制etopo浮雕图像作为地图背景 -
warpimage()
:将用户提供的图像投影到地图上
附:basemap更多绘图例子