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目录下。
地图文件列表.jpg (40.09 KB, 下载次数: 3)
下载附件 保存到相册
地图书文件列表
2015-3-25 19:40 上传
准备工作完成后,进入正题: (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()
复制代码
地图1.png (49.05 KB, 下载次数: 5)
下载附件 保存到相册
2015-3-25 19:54 上传
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)
复制代码
地图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对极坐标投影没有自动标注经纬度功能,需要自己设置,调了好久标注位置,请大家切用且珍惜哈。
- 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)
复制代码
北极地图.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即调用到了自己的地图,不用担心边界不准确了。 不知道为什么最后多出来一张图片,麻烦版主帮忙删了吧。 |