[原创]python绘制地图的利器Cartopy使用说明

阿阿飞飞 · 发表于 2015-3-25 20:22:03

登录后查看更多精彩内容~

您需要登录才可以下载或查看，没有帐号？立即注册

x

本帖最后由阿阿飞飞于 2015-3-25 21:08 编辑

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，输入以下代码：

import matplotlib.pyplot as plt
import cartopy.crs as ccrs
plt.figure(figsize=(6, 3))
ax = plt.axes(projection=ccrs.PlateCarree(central_longitude=180))
ax.coastlines(resolution='110m')
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目录下。

地图书文件列表

准备工作完成后，进入正题：
（3）绘制地图
前面两步虽有些繁琐，但会一劳永逸，下面是示例程序，scale参数用于调整使用哪种分辨率的地图，全球建议用1：110的，小区域可以用1:50的或1:10的。
a）绘制全球地图程序：

#===================================================
#使用cartopy绘制地图
#需要从http://www.naturalearthdata.com/downloads/下载shape文件
#下载后，解压缩，文件名统一去掉"ne_"开头，拷贝至D:\Program Files\
#WinPython-32bit-2.7.9.3\settings\.local\share\cartopy\shapefiles\natural_earth\physical\
#路径下面，coastline文件对应ax.coastlines命令，land文件对应land命令
#===================================================
scale='110m'
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cfeature
from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter
fig=plt.figure(figsize=(8, 10))
ax=plt.axes(projection=ccrs.PlateCarree(central_longitude=180))
ax.set_global()
#===================================================
#需要填充陆地颜色时使用
#ax.add_feature(cfeature.LAND, facecolor='0.75') #默认为110m，其它分辨率需用下面命令
land = cfeature.NaturalEarthFeature('physical', 'land', scale,edgecolor='face',
facecolor=cfeature.COLORS['land'])
ax.add_feature(land, facecolor='0.75')
#===================================================
#改变ax.add_feature和ax.coastlines的先后使用顺序可实现边界线的显示或完全填充覆盖
ax.coastlines(scale)
#===================================================
#标注坐标轴
ax.set_xticks([0, 60, 120, 180, 240, 300, 360], crs=ccrs.PlateCarree())
ax.set_yticks([-90, -60, -30, 0, 30, 60, 90], crs=ccrs.PlateCarree())
#zero_direction_label用来设置经度的0度加不加E和W
lon_formatter = LongitudeFormatter(zero_direction_label=False)
lat_formatter = LatitudeFormatter()
ax.xaxis.set_major_formatter(lon_formatter)
ax.yaxis.set_major_formatter(lat_formatter)
#添加网格线
gl = ax.gridlines()

复制代码

b）绘制全球地图（函数形式）

#===================================================
#函数形式，调用cartopy，绘制全球地图
#===================================================
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cfeature
from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter
def make_map(scale):
fig=plt.figure(figsize=(8, 10))
ax=plt.axes(projection=ccrs.PlateCarree(central_longitude=180))
ax.set_global()
land = cfeature.NaturalEarthFeature('physical', 'land', scale,edgecolor='face',
facecolor=cfeature.COLORS['land'])
ax.add_feature(land, facecolor='0.75')
ax.coastlines(scale)
#标注坐标轴
ax.set_xticks([0, 60, 120, 180, 240, 300, 360], crs=ccrs.PlateCarree())
ax.set_yticks([-90, -60, -30, 0, 30, 60, 90], crs=ccrs.PlateCarree())
#zero_direction_label用来设置经度的0度加不加E和W
lon_formatter = LongitudeFormatter(zero_direction_label=False)
lat_formatter = LatitudeFormatter()
ax.xaxis.set_major_formatter(lon_formatter)
ax.yaxis.set_major_formatter(lat_formatter)
#添加网格线
#gl = ax.gridlines()
ax.grid()
return fig,ax
fig,ax=make_map(scale='110m')

复制代码

此程序结果与上面一样。
c）绘制区域地图（函数形式）

#===================================================
#函数形式，调用cartopy，绘制区域地图
#===================================================
import numpy as np
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cfeature
from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter
def make_map(scale,box,xstep,ystep):
fig=plt.figure(figsize=(8, 10))
ax=plt.axes(projection=ccrs.PlateCarree())
#set_extent需要配置相应的crs，否则出来的地图范围不准确
ax.set_extent(box,crs=ccrs.PlateCarree())
land = cfeature.NaturalEarthFeature('physical', 'land', scale,edgecolor='face',
facecolor=cfeature.COLORS['land'])
ax.add_feature(land, facecolor='0.75')
ax.coastlines(scale)
#===================================================
#图像地址D:\Program Files\WinPython-32bit-2.7.9.3\python-2.7.9\Lib\site-packages\
#cartopy\data\raster\natural_earth\50-natural-earth-1-downsampled.png
#如果有其它高精度图像文件，改名替换即可
ax.stock_img()
#===================================================
#标注坐标轴
ax.set_xticks(np.arange(box[0],box[1]+xstep,xstep), crs=ccrs.PlateCarree())
ax.set_yticks(np.arange(box[2],box[3]+ystep,ystep), crs=ccrs.PlateCarree())
#zero_direction_label用来设置经度的0度加不加E和W
lon_formatter = LongitudeFormatter(zero_direction_label=False)
lat_formatter = LatitudeFormatter()
ax.xaxis.set_major_formatter(lon_formatter)
ax.yaxis.set_major_formatter(lat_formatter)
#添加网格线
ax.grid()
return fig,ax
box=[100,150,0,50]
fig,ax=make_map(scale='50m',box=box,xstep=10,ystep=10)

复制代码

（4）绘制极地投影地图
绘制该地图需要使用最新版的cartopy-0.12版本，windows下暂无此版本，可以安装好0.11后，删除cartopy安装目录下的mpl文件夹下的全部文件，然后拷贝0.12目录cartopy-0.12.0rc1\lib\cartopy\mpl文件夹下的全部文件进行替换，即可使用新版的功能。此程序长了点，因为Cartopy对极坐标投影没有自动标注经纬度功能，需要自己设置，调了好久标注位置，请大家切用且珍惜哈。

import matplotlib.path as mpath
import matplotlib.pyplot as plt
import numpy as np
import matplotlib.ticker as mticker
import cartopy.crs as ccrs
import cartopy.feature as cfeature
fig=plt.figure(figsize=(6, 6))
ax=plt.axes(projection=ccrs.NorthPolarStereo())
box=[-180, 180, 55, 90]
xstep,ystep=30,15
# Limit the map to -60 degrees latitude and below.
ax.set_extent(box, crs=ccrs.PlateCarree())
scale='50m'
land = cfeature.NaturalEarthFeature('physical', 'land', scale,edgecolor='face',
facecolor=cfeature.COLORS['land'])
ocean = cfeature.NaturalEarthFeature('physical', 'ocean', scale,edgecolor='face',
facecolor=cfeature.COLORS['water'])
ax.add_feature(land,facecolor='0.75')
ax.add_feature(ocean,facecolor='blue')
ax.coastlines(scale,linewidth=0.9)
#标注坐标轴
line=ax.gridlines(draw_labels=False)
line.ylocator=mticker.FixedLocator(np.arange(40,90,20))#手动设置x轴刻度
line.xlocator=mticker.FixedLocator(np.arange(-180,210,30))#手动设置x轴刻度
# Compute a circle in axes coordinates, which we can use as a boundary
# for the map. We can pan/zoom as much as we like - the boundary will be
# permanently circular.
theta = np.linspace(0, 2*np.pi, 100)
center, radius = [0.5, 0.5], 0.5
verts = np.vstack([np.sin(theta), np.cos(theta)]).T
circle = mpath.Path(verts * radius + center)
ax.set_boundary(circle, transform=ax.transAxes)
#创建要标注的labels字符串
ticks=np.arange(0,210,30)
etick=['0']+['%d$^\circ$E'%tick for tick in ticks if (tick !=0) & (tick!=180)]+['180']
wtick=['%d$^\circ$W'%tick for tick in ticks if (tick !=0) & (tick!=180)]
labels=etick+wtick
#创建与labels对应的经纬度标注位置
#xticks=[i for i in np.arange(0,210,30)]+[i for i in np.arange(-32,-180,-30)]
xticks=[-0.8,28,58,89.1,120,151,182.9,-36,-63,-89,-114,-140]
yticks=[53]+[53]+[54]+[55]*2+[54.5]+[54]+[50]+[49]*3+[50.6]
#标注经纬度
#ax.text(0.01,0.23,'60$^\circ$W',transform=ax.transAxes,rotation=25)
#ax.text(-63,50,'60$^\circ$W',transform=ccrs.Geodetic(),rotation=25)
for xtick,ytick,label in zip(xticks,yticks,labels):
ax.text(xtick,ytick,label,transform=ccrs.Geodetic())
x=[180, 180, 0, 0]
y=[50, 90, 90, 50]
ax.plot([-180,0],[80,80],':',transform=ccrs.Geodetic(),color='k',linewidth=0.4)
ax.plot([-90,90],[80,80],':',transform=ccrs.Geodetic(),color='k',linewidth=0.5)
#ax.plot([90,0],[50,50],'-.',transform=ccrs.Geodetic(),color='r',linewidth=6)
ax.text(11.9333,78.9166,r'$\bigstar,transform=ccrs.Geodetic(),size=15,color='r')
fig.savefig(u'c:\\北极.png',dpi=300)

复制代码

北极投影地图

介绍完毕，写一个帖子还真挺耗时间，不多说了，更多的功能，请去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即调用到了自己的地图，不用担心边界不准确了。
不知道为什么最后多出来一张图片，麻烦版主帮忙删了吧。

chongzika · 发表于 2017-4-13 15:57:32

chongzika 发表于 2017-4-13 15:41
linux系统下的怎么设置呢

它会默认下载到.local/share/cartopy/shapefiles/

yunpiao · 发表于 2015-3-25 21:11:25

yunpiao · 发表于 2015-3-25 21:11:37

平流层的萝卜 · 发表于 2015-3-25 21:57:21

非常好，感谢分享！

风之精灵 · 发表于 2015-3-26 09:06:07

不错哎~学习啦，谢谢楼主分享

乄belongU℡ · 发表于 2015-3-26 10:26:15

请教一下楼主，Grads能画出这个图吗？ forum.php.png

阿阿飞飞 · 发表于 2015-3-26 12:06:29

乄belongU℡ 发表于 2015-3-26 10:26
请教一下楼主，Grads能画出这个图吗？

应该可以的，grads也可以直接读取绘制shape地图，具体请去grads区看看

乄belongU℡ · 发表于 2015-3-27 08:57:24

好的。

taxueyueming · 发表于 2015-3-27 12:25:32

很好，感谢分享

Aplomb-天元 · 发表于 2015-3-28 14:42:05

python 和ncl 可以连用me。哪个用起来处理nc文件可能更好呢~

		自动登录	找回密码
密码			立即注册

[源代码] [原创]python绘制地图的利器Cartopy使用说明

登录后查看更多精彩内容~

评分