python绘制地图的利器Cartopy使用说明

Posted mishidemudong

tags:

篇首语:本文由小常识网(cha138.com)小编为大家整理,主要介绍了python绘制地图的利器Cartopy使用说明相关的知识,希望对你有一定的参考价值。

python绘制地图一般使用Basemap绘图包,但该包配置相对较繁琐,自定义性不强,这里介绍一个绘制地图的利器Cartopy,个人认为该工具方便、快捷,附上一些自己写的程序。
准备工作,工欲善其事,必先利其器
(1)先下载主角:Cartopy
a)下载地址:
linux平台直接去官网下载:http://scitools.org.uk/cartopy/download.html
windows平台下的Cartopy下载地址:
http://www.lfd.uci.edu/~gohlke/pythonlibs/#cartopy
还需要一些必须的软件包:Shapely、pyshp也都从上面的网址下载
官网自带的示例网址:
http://scitools.org.uk/cartopy/docs/latest/gallery.html

b)Cartopy下载安装完毕后,打开python,输入以下代码:
  1. import matplotlib.pyplot as plt
  2. import cartopy.crs as ccrs

  3. plt.figure(figsize=(6, 3))
  4. ax = plt.axes(projection=ccrs.PlateCarree(central_longitude=180))
  5. ax.coastlines(resolution='110m')
  6. ax.gridlines()
复制代码 运行后,如果顺利导入了cartopy包,不管有没有出现地图,都不要紧,第二步走起。
(2)再下载配角:地图数据
下载地址:http://www.naturalearthdata.com/downloads/
里面有三种分辨率的shape地图数据可选,方便起见,分别下载三种分辨率中的physical数据中的Coastline和Land数据,每种数据下载后都是一个压缩包,如下载了1:10分辨率的physical中的coastline数据压缩包:ne_10m_coastline.zip,解压缩后有6个文件,其中“ne_10m_coastline.README”和“ne_10m_coastline.VERSION”可直接删除,剩余4个,进行改名操作,扩展名前面的文件名,如“ne_10m_coastline”,修改为“10m_coastline”,即去掉“ne_”,4个文件分别这样更改。再下载1:50和1:110的文件分别进行此操作。所有地图文件下载、解压、更名完毕后,拷贝到一个文件夹下。我的文件夹列表如下图,把这些文件全选(注意进入文件夹目录,全选文件,不带文件夹),复制粘贴到D:\\Program Files\\WinPython-32bit-2.7.9.3\\settings\\.local\\share\\cartopy\\shapefiles\\natural_earth\\physical 目录下(该目录根据自己所装的python而定,运行(1)中的程序后,程序会自动创建physical文件夹,具体该文件夹在哪,搜索电脑文件找找看),我安装的是winpython2.7.9.3,physical目录就位于上面这个目录中,所以我把所有shape地图文件拷贝到了该physical目录下。

地图文件列表.jpg (40.09 KB, 下载次数: 3)

下载附件  保存到相册

地图书文件列表

2015-3-25 19:40 上传



准备工作完成后,进入正题:
(3)绘制地图
前面两步虽有些繁琐,但会一劳永逸,下面是示例程序,scale参数用于调整使用哪种分辨率的地图,全球建议用1:110的,小区域可以用1:50的或1:10的。
a)绘制全球地图程序:
  1. #===================================================
  2. #使用cartopy绘制地图
  3. #需要从http://www.naturalearthdata.com/downloads/下载shape文件
  4. #下载后,解压缩,文件名统一去掉"ne_"开头,拷贝至D:\\Program Files\\
  5. #WinPython-32bit-2.7.9.3\\settings\\.local\\share\\cartopy\\shapefiles\\natural_earth\\physical\\
  6. #路径下面,coastline文件对应ax.coastlines命令,land文件对应land命令
  7. #===================================================
  8. scale='110m'
  9. import matplotlib.pyplot as plt
  10. import cartopy.crs as ccrs
  11. import cartopy.feature as cfeature
  12. from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter
  13. fig=plt.figure(figsize=(8, 10))
  14. ax=plt.axes(projection=ccrs.PlateCarree(central_longitude=180))
  15. ax.set_global()
  16. #===================================================
  17. #需要填充陆地颜色时使用
  18. #ax.add_feature(cfeature.LAND, facecolor='0.75') #默认为110m,其它分辨率需用下面命令
  19. land = cfeature.NaturalEarthFeature('physical', 'land', scale,edgecolor='face',
  20.                                                               facecolor=cfeature.COLORS['land'])
  21. ax.add_feature(land, facecolor='0.75')
  22. #===================================================
  23. #改变ax.add_feature和ax.coastlines的先后使用顺序可实现边界线的显示或完全填充覆盖
  24. ax.coastlines(scale)
  25. #===================================================
  26. #标注坐标轴
  27. ax.set_xticks([0, 60, 120, 180, 240, 300, 360], crs=ccrs.PlateCarree())
  28. ax.set_yticks([-90, -60, -30, 0, 30, 60, 90], crs=ccrs.PlateCarree())
  29. #zero_direction_label用来设置经度的0度加不加E和W
  30. lon_formatter = LongitudeFormatter(zero_direction_label=False)
  31. lat_formatter = LatitudeFormatter()
  32. ax.xaxis.set_major_formatter(lon_formatter)
  33. ax.yaxis.set_major_formatter(lat_formatter)
  34. #添加网格线
  35. gl = ax.gridlines()
复制代码

地图1.png (49.05 KB, 下载次数: 5)

下载附件  保存到相册

2015-3-25 19:54 上传



b)绘制全球地图(函数形式)
  1. #===================================================
  2. #函数形式,调用cartopy,绘制全球地图
  3. #===================================================
  4. import matplotlib.pyplot as plt
  5. import cartopy.crs as ccrs
  6. import cartopy.feature as cfeature
  7. from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter
  8. def make_map(scale):
  9.     fig=plt.figure(figsize=(8, 10))
  10.     ax=plt.axes(projection=ccrs.PlateCarree(central_longitude=180))
  11.     ax.set_global()
  12.     land = cfeature.NaturalEarthFeature('physical', 'land', scale,edgecolor='face',
  13.                                                               facecolor=cfeature.COLORS['land'])
  14.     ax.add_feature(land, facecolor='0.75')
  15.     ax.coastlines(scale)
  16.     #标注坐标轴
  17.     ax.set_xticks([0, 60, 120, 180, 240, 300, 360], crs=ccrs.PlateCarree())
  18.     ax.set_yticks([-90, -60, -30, 0, 30, 60, 90], crs=ccrs.PlateCarree())
  19.     #zero_direction_label用来设置经度的0度加不加E和W
  20.     lon_formatter = LongitudeFormatter(zero_direction_label=False)
  21.     lat_formatter = LatitudeFormatter()
  22.     ax.xaxis.set_major_formatter(lon_formatter)
  23.     ax.yaxis.set_major_formatter(lat_formatter)
  24.     #添加网格线
  25.     #gl = ax.gridlines()
  26.     ax.grid()
  27.     return fig,ax
  28. fig,ax=make_map(scale='110m')
复制代码 此程序结果与上面一样。
c)绘制区域地图(函数形式)
  1. #===================================================
  2. #函数形式,调用cartopy,绘制区域地图
  3. #===================================================
  4. import numpy as np
  5. import matplotlib.pyplot as plt
  6. import cartopy.crs as ccrs
  7. import cartopy.feature as cfeature
  8. from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter
  9. def make_map(scale,box,xstep,ystep):
  10.     fig=plt.figure(figsize=(8, 10))
  11.     ax=plt.axes(projection=ccrs.PlateCarree())
  12.     #set_extent需要配置相应的crs,否则出来的地图范围不准确
  13.     ax.set_extent(box,crs=ccrs.PlateCarree())
  14.     land = cfeature.NaturalEarthFeature('physical', 'land', scale,edgecolor='face',
  15.                                                               facecolor=cfeature.COLORS['land'])
  16.     ax.add_feature(land, facecolor='0.75')
  17.     ax.coastlines(scale)
  18.     #===================================================
  19.     #图像地址D:\\Program Files\\WinPython-32bit-2.7.9.3\\python-2.7.9\\Lib\\site-packages\\
  20.     #cartopy\\data\\raster\\natural_earth\\50-natural-earth-1-downsampled.png
  21.     #如果有其它高精度图像文件,改名替换即可
  22.     ax.stock_img()
  23.     #===================================================
  24.     #标注坐标轴
  25.     ax.set_xticks(np.arange(box[0],box[1]+xstep,xstep), crs=ccrs.PlateCarree())
  26.     ax.set_yticks(np.arange(box[2],box[3]+ystep,ystep), crs=ccrs.PlateCarree())
  27.     #zero_direction_label用来设置经度的0度加不加E和W
  28.     lon_formatter = LongitudeFormatter(zero_direction_label=False)
  29.     lat_formatter = LatitudeFormatter()
  30.     ax.xaxis.set_major_formatter(lon_formatter)
  31.     ax.yaxis.set_major_formatter(lat_formatter)
  32.     #添加网格线
  33.     ax.grid()
  34.     return fig,ax
  35. box=[100,150,0,50]
  36. fig,ax=make_map(scale='50m',box=box,xstep=10,ystep=10)
复制代码

地图2.png (149.94 KB, 下载次数: 2)

下载附件  保存到相册

2015-3-25 19:57 上传



(4)绘制极地投影地图
绘制该地图需要使用最新版的cartopy-0.12版本,windows下暂无此版本,可以安装好0.11后,删除cartopy安装目录下的mpl文件夹下的全部文件,然后拷贝0.12目录cartopy-0.12.0rc1\\lib\\cartopy\\mpl文件夹下的全部文件进行替换,即可使用新版的功能。此程序长了点,因为Cartopy对极坐标投影没有自动标注经纬度功能,需要自己设置,调了好久标注位置,请大家切用且珍惜哈。
  1. import matplotlib.path as mpath
  2. import matplotlib.pyplot as plt
  3. import numpy as np
  4. import matplotlib.ticker as mticker
  5. import cartopy.crs as ccrs
  6. import cartopy.feature as cfeature

  7. fig=plt.figure(figsize=(6, 6))
  8. ax=plt.axes(projection=ccrs.NorthPolarStereo())
  9. box=[-180, 180, 55, 90]
  10. xstep,ystep=30,15
  11.     # Limit the map to -60 degrees latitude and below.
  12. ax.set_extent(box, crs=ccrs.PlateCarree())
  13. scale='50m'
  14. land = cfeature.NaturalEarthFeature('physical', 'land', scale,edgecolor='face',
  15.                                                               facecolor=cfeature.COLORS['land'])
  16. ocean = cfeature.NaturalEarthFeature('physical', 'ocean', scale,edgecolor='face',
  17.                                                              facecolor=cfeature.COLORS['water'])
  18. ax.add_feature(land,facecolor='0.75')
  19. ax.add_feature(ocean,facecolor='blue')
  20. ax.coastlines(scale,linewidth=0.9)
  21. #标注坐标轴
  22. line=ax.gridlines(draw_labels=False)

  23. line.ylocator=mticker.FixedLocator(np.arange(40,90,20))#手动设置x轴刻度
  24. line.xlocator=mticker.FixedLocator(np.arange(-180,210,30))#手动设置x轴刻度
  25.     # Compute a circle in axes coordinates, which we can use as a boundary
  26.     # for the map. We can pan/zoom as much as we like - the boundary will be
  27.     # permanently circular.
  28. theta = np.linspace(0, 2*np.pi, 100)
  29. center, radius = [0.5, 0.5], 0.5
  30. verts = np.vstack([np.sin(theta), np.cos(theta)]).T
  31. circle = mpath.Path(verts * radius + center)
  32. ax.set_boundary(circle, transform=ax.transAxes)

  33. #创建要标注的labels字符串
  34. ticks=np.arange(0,210,30)
  35. etick=['0']+['%d$^\\circ$E'%tick for tick in ticks if (tick !=0) & (tick!=180)]+['180']
  36. wtick=['%d$^\\circ$W'%tick for tick in ticks if (tick !=0) & (tick!=180)]
  37. labels=etick+wtick
  38. #创建与labels对应的经纬度标注位置
  39. #xticks=[i for i in np.arange(0,210,30)]+[i for i in np.arange(-32,-180,-30)]
  40. xticks=[-0.8,28,58,89.1,120,151,182.9,-36,-63,-89,-114,-140]
  41. yticks=[53]+[53]+[54]+[55]*2+[54.5]+[54]+[50]+[49]*3+[50.6]

  42. #标注经纬度   
  43. #ax.text(0.01,0.23,'60$^\\circ$W',transform=ax.transAxes,rotation=25)
  44. #ax.text(-63,50,'60$^\\circ$W',transform=ccrs.Geodetic(),rotation=25)
  45. for xtick,ytick,label in zip(xticks,yticks,labels):
  46.     ax.text(xtick,ytick,label,transform=ccrs.Geodetic())
  47. x=[180, 180, 0, 0]
  48. y=[50, 90, 90, 50]
  49. ax.plot([-180,0],[80,80],':',transform=ccrs.Geodetic(),color='k',linewidth=0.4)
  50. ax.plot([-90,90],[80,80],':',transform=ccrs.Geodetic(),color='k',linewidth=0.5)
  51. #ax.plot([90,0],[50,50],'-.',transform=ccrs.Geodetic(),color='r',linewidth=6)

  52. ax.text(11.9333,78.9166,r'$\\bigstar,transform=ccrs.Geodetic(),size=15,color='r')
  53. fig.savefig(u'c:\\\\北极.png',dpi=300)
复制代码

北极地图.png (64.49 KB, 下载次数: 3)

下载附件  保存到相册

北极投影地图

2015-3-25 20:24 上传



介绍完毕,写一个帖子还真挺耗时间,不多说了,更多的功能,请去cartopy官网http://scitools.org.uk/cartopy/docs/latest/gallery.html,看更多的例子。
另外一个notebook网址也很多示例可以借鉴:http://nbviewer.ipython.org/gist/pelson

注:如果想使用自己的地图,比如我有一个带我国法定边界的地图shape文件,名为:World.shp、World.shx、World.dbf、World.prj,重新命名为10m_coastline或10m_land文件(这里要根据该shape文件是line型还是polygon型),替换原目录下的shape文件,画图时使用scale为10m即调用到了自己的地图,不用担心边界不准确了。
不知道为什么最后多出来一张图片,麻烦版主帮忙删了吧。

北极地图.png (64.49 KB, 下载次数: 1)

下载附件  保存到相册

2015-3-25 20:11 上传


以上是关于python绘制地图的利器Cartopy使用说明的主要内容,如果未能解决你的问题,请参考以下文章

Python 3.4在生成一些 - 但不是全部 - 具有分段错误11的Cartopy地图时崩溃

使用 cartopy 绘制来自 netcdf 文件的 4 维变量的数据

matplotlib+cartopy+geopandas,实现专业地图可视化!

matplotlib+cartopy+geopandas,实现专业地图可视化!

matplotlib+cartopy+geopandas,实现专业地图可视化!

matplotlib+cartopy+geopandas,实现专业地图可视化!