前言

嗨,你好,我是来自点点GIS的南南

我与Cartopy的认识起源于“气象水文科研猫”的这个推文,那时候的我觉得,用代码画地图好酷,arcgis就感觉low爆了。但是一直没有时间学习。前段时间放暑假,磕磕绊绊装完包以后就不想动弹了,折腾环境折腾的我头皮发麻。gdal和cartopy真的是我学python以来装的最麻烦的包。

很多时候我想学习,却感觉无从下手,大佬们的博文一大堆,我看了以后只有一句“卧槽,牛逼”,然后没然后了。前几天某人欺负我,把我整抑郁了,大晚上睡不着觉,就想着填个坑,虽然我技术菜,不能“炫技”,但写个入门教程难不倒我吧。然后就有了这篇文章。

这篇文章参考了诸多大佬的博文,如气象学家,云台书使,气象学人,好奇心Log,等等公众号大佬。我作为初学者,写的文肯定错误不少,欢迎大家给我留言指正,要是有老板叫我去打工就更好了,嘿嘿嘿

Emil:2531788189@qq.com

Cartopy介绍

Cartopy 是一个开源免费的第三方 Python 扩展包,由英国气象办公室的科学家们开发,支持 Python 2.7 和 Python 3,致力于使用最简单直观的方式生成地图,并提供对 matplotlib 友好的协作接口。

在常用的python绘图库中,basemap,geopandas,pyecharts等,其中basemap已经停止维护了,在前文中我已经讲到过,pyecharts是用于数据可视化等专业图表的绘制,之前我的文章也介绍过,pyecharts在地学可视化中功能实在过于简陋;geopandas是基于pandas的,一般用于商业数据分析,所以我更推荐大家学习Cartopy,虽然Cartopy气象专业用得很多哈哈哈

Cartopy安装

https://mp.weixin.qq.com/s/GW6odDBGLPVRZTKx0YZLVg

Cartopy绘图入门

在Cartopy的使用过程中,我们常常需要搭配其他库一起使用,常用的有如Matplotlib ,numpy,pandas,gma等

Cartopy投影

Cartopy 是利用 Matplotlib 来画图的,因此首先要导入 pyplot 模块。在 Cartopy 中,每种投影都是一个类,被存放在 cartopy.crs 模块中,crs 即坐标参考系统(Coordinate Reference Systems)之意。所以接着要导入这个模块。

投影用法示例:

proj = ccrs.PlateCarree()   #proj = ccrs.投影()

这里简单说明一下常用的三种投影,Cartopy 大概有三四十种投影,如国想详细了解可以参考下面的文章

https://scitools.org.uk/cartopy/docs/v0.15/crs/projections.html
默认投影PlateCarree
墨卡托投影Mercator
兰勃脱投影Lambert

默认投影适合单独省份或者地级市的绘制,兰勃脱投影适合中纬度大范围绘制,比如绘制全中国大公鸡、东亚形势、西北太平洋等;墨卡托投影适合低纬度赤道附近的绘制,一般研究台风、纬向环流等

绘制地图示例

将 cartopy 集成到 matplotlib 的主要类是 GeoAxes,它是普通 matplotlib Axes 的子类。GeoAxes 类为特定于绘制地图的轴添加了额外的功能。创建一个 GeoAxes 对象的办法是,在创建 axes(或 subplot)时,通过参数 projection 指定一个 ccrs 中的投影。这里便利用这一方法生成了一个等距圆柱投影下的 ax。

最后调用 ax 的方法 coastlines 画出海岸线,默认以本初子午线为中心,比例尺为 1:110m(m 表示 million)。

因此用 Cartopy 画地图的基本流程并不复杂:

  • 创建画布。

  • 通过指定 projection 参数,创建 GeoAxes 对象。

  • 调用 GeoAxes 的方法画图。

  • GeoAxes 用法扩展(部分常用)
    • set_global:让地图的显示范围扩展至投影的最大范围。例如,对 PlateCarree 投影的 ax 使用后,地图会变成全球的。
    • set_extent:给出元组 (x0, x1, y0, y1) 以限制地图的显示范围。
    • set_xticks:设置 x 轴的刻度。
    • set_yticks:设置 y 轴的刻度。
    • gridlines:给地图添加网格线。
    • coastlines:在地图上绘制海岸线。
    • stock_img:给地图添加低分辨率的地形图背景。
    • add_feature:给地图添加特征(例如陆地或海洋的填充、河流等)。
import matplotlib.pyplot as plt###引入库包
import cartopy.crs as ccrs
proj = ccrs.PlateCarree()
fig = plt.figure(figsize=(4, 4), dpi=200)  # 创建画布
ax = plt.axes(projection=ccrs.PlateCarree())# 创建子图
ax.coastlines()# 调用ax的方法画海岸线
'''
fig语法:
figure(num=None, figsize=None, dpi=None, facecolor=None, edgecolor=None, frameon=True)
num:图像编号或名称,数字为编号 ,字符串为名称
figsize:指定figure的宽和高,单位为英寸;
dpi参数指定绘图对象的分辨率,即每英寸多少个像素,缺省值为80      1英寸等于2.5cm,A4纸是 21*30cm的纸张 
facecolor:背景颜色
edgecolor:边框颜色
frameon:是否显示边框
'''

image-20220718224019347

如需修改中央经线,可以在创建子图的时候指定一些参数

ax = plt.axes(projection=ccrs.PlateCarree(central_longitude=180))#全球图像的中央位于太平洋的 180 度经线处

image-20220718224340658

也可以在绘制线条的时候指定线宽

ax.coastlines(lw=0.3)#线条宽度0.3

image-20220718224700994

在地图上添加特征

上面的地图是一个十分粗糙的地图,仅有海岸线,尝试一下添加更多地理信息

添加预定义要素

首先需要导入一个cartopy.feature 常量,为了简化一些非常常见的情况,如大陆海洋国界等cartopy都已经进行了预定义,但需要注意的是直接导入的中国国界线并不是标准的,在使用时需要注意

这些信息带有三种分辨率,分别为110m,50m和10m。

image-20220718225823227

import matplotlib.pyplot as plt###引入库包
import cartopy.crs as ccrs
import cartopy.feature as cfeature#预定义常量

proj = ccrs.PlateCarree()
fig = plt.figure(figsize=(4, 4), dpi=200)  # 创建画布
ax = plt.axes(projection=ccrs.PlateCarree(central_longitude=180))#全球图像的中央位于太平洋的 180 度经线处

ax.add_feature(cfeature.LAND)#添加陆地
ax.add_feature(cfeature.COASTLINE,lw = 0.3)#添加海岸线
ax.add_feature(cfeature.RIVERS,lw = 0.25)#添加河流
#ax.add_feature(cfeat.RIVERS.with_scale('50m'),lw = 0.25)  # 加载分辨率为50的河流
ax.add_feature(cfeature.LAKES)#添加湖泊
ax.add_feature(cfeature.BORDERS, linestyle = '-',lw = 0.25)#不推荐,因为该默认参数会使得我国部分领土丢失
ax.add_feature(cfeature.OCEAN)#添加海洋

ax.coastlines(lw=0.3)

image-20220718230518430

在cartopy中,所有的线条都可以改变他的粗细。只用改变其参数lw = *即可

ax.add_feature(cfeature.RIVERS,lw = 2)#改变河流粗细

image-20220718231338174

同样的,所有的线,面也可以改变它的颜色

ax.add_feature(cfeature.LAKES,color='r')#指定湖泊颜色为红色

image-20220718232653107

添加经纬度标签

这里还是要用到前文说过的GeoAxes 用法

  • set_xticks:设置 x 轴的刻度。
  • set_yticks:设置 y 轴的刻度。

例:x轴为(-180,-120,-60,0,60,120,180),y轴为(-90,-60, -30, 0, 30, 60,90)

x_extent = [ 0, 60, 120, 180, 240, 300, 360]#经纬度范围,#直接使用(-180,-120,-60,0,60,120,180)会异常,需要写成(0, 60, 120, 180, 240, 300, 360)的形式
y_extent = [ -90,-60, -30, 0, 30, 60,90]

fig = plt.figure(figsize=(4, 4), dpi=200)  
ax = plt.axes(projection=ccrs.PlateCarree(central_longitude=180))

ax.add_feature(cfeature.LAND)#添加陆地
ax.add_feature(cfeature.COASTLINE,lw = 0.3)#添加海岸线
ax.add_feature(cfeature.RIVERS,lw = 0.25)#添加河流
ax.add_feature(cfeature.LAKES)#指定湖泊颜色为红色#添加湖泊
#ax.add_feature(cfeature.BORDERS, linestyle = '-',lw = 0.25)#不推荐,因为该默认参数会使得我国部分领土丢失
ax.add_feature(cfeature.OCEAN)#添加海洋

ax.set_xticks(x_extent, crs=ccrs.PlateCarree())#添加经纬度
ax.set_yticks(y_extent, crs=ccrs.PlateCarree())

image-20220718233555539

将经纬度标签转换为具有单位的形式

from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter #导入Cartopy专门提供的经纬度的Formatter

x_extent = [ 0, 60, 120, 180, 240, 300, 360]#经纬度范围,#直接使用(-180,-120,-60,0,60,120,180)会异常,需要写成(0, 60, 120, 180, 240, 300, 360)的形式
y_extent = [ -90,-60, -30, 0, 30, 60,90]

fig = plt.figure(figsize=(4, 4), dpi=200)  
ax = plt.axes(projection=ccrs.PlateCarree(central_longitude=180))

ax.add_feature(cfeature.LAND)#添加陆地
ax.add_feature(cfeature.COASTLINE,lw = 0.3)#添加海岸线
ax.add_feature(cfeature.RIVERS,lw = 0.25)#添加河流
ax.add_feature(cfeature.LAKES)#指定湖泊颜色为红色#添加湖泊
#ax.add_feature(cfeature.BORDERS, linestyle = '-',lw = 0.25)#不推荐,因为该默认参数会使得我国部分领土丢失
ax.add_feature(cfeature.OCEAN)#添加海洋

ax.set_xticks(x_extent, crs=ccrs.PlateCarree())#添加经纬度
ax.set_yticks(y_extent, crs=ccrs.PlateCarree())

# 利用Formatter格式化刻度标签
lon_formatter = LongitudeFormatter(zero_direction_label=False)
lat_formatter = LatitudeFormatter()
ax.xaxis.set_major_formatter(lon_formatter)
ax.yaxis.set_major_formatter(lat_formatter)

image-20220718235401412

经纬网

在前文中GeoAxes 用法提到过添加经纬网的函数gridlines,只用在上述代码末尾中添加ax.gridlines(),即

ax.gridlines()#添加网格

image-20220718235745298

这实在是太丑了,给它换个线样式

ax.gridlines(linestyle='--')#添加网格,线样式为'--'

image-20220718235924923

配上网格总感觉刻度朝外怪怪的

ax.tick_params(color = 'blue',direction='in')#更改刻度指向为朝内,颜色设置为蓝色

image-20220719000410778

更多GeoAxes 用法可以参考以下文章

https://scitools.org.uk/cartopy/docs/v0.14/matplotlib/feature_interface.html

局部地图

很多时候我们绘图不会用到世界地图,所以可以尝试更改你的范围。cartopy提供更改的方式是通过定义经纬度范围来进行更改的。在上述代码末尾添加下方样例代码试试

set_extent:给出元组 (x0, x1, y0, y1) 以限制地图的显示范围。

ax.set_extent([0,180,0,35],crs = ccrs.PlateCarree()) #选取经度为0°E-180°E,纬度为0°N-90°N的区域

image-20220719000849112

添加标题

ax.set_title('Cartopy') #添加标题Cartopy

image-20220719002030473

Cartopy绘图进阶

在前文中提到过,Cartopy的中国地图边界是有问题的,那么在日常使用中,我们该如何避免这些问题呢?

Cartopy中国地图绘制方法

方法一

用 Cartopy 的 shapereader 读取后即可绘制。 Cartopy 的 shapereader用法示例可参考以下文章

https://scitools.org.uk/cartopy/docs/v0.15/tutorials/using_the_shapereader.html

关于正确的行政边界获取可以参考以下博文

https://mp.weixin.qq.com/s/aUCjTRrV6Cz-k7mtXOEiGw

使用cartopy读取shapefile绘制共有两种方法,分别是add_featureadd_geometries,但无论你使用哪种方法,你都需要先读取文件才能够添加。Cartopy提供了一个基于pyshp的接口以实现对地理文件的简单读取和操作:

首先导入相关包

import numpy as np
import cartopy.crs as ccrs
import matplotlib.pyplot as plt
import cartopy.feature as cfeat
import cartopy.io.shapereader as shapereader

你找他总得给他来一个地址吧

shp_path = "D:/data2/china/china/province_9south.shp"#切记路径一定要英文,血的教训!!!

绘图三部曲,画布,投影,子图

fig  = plt.figure(figsize=(6,6)) #创建画布
proj = ccrs.PlateCarree()#创建投影,选择cartopy的默认投影
ax   = fig.subplots(1, 1, subplot_kw={'projection': proj})  #子图
extent = [70,140,2,60]#限定绘图范围
add_feature加载自定义地图
china_map = cfeat.ShapelyFeature(shapereader.Reader(china_map).geometries(), proj, edgecolor='k', facecolor='none')
ax.add_feature(china_map, linewidth=1)
ax.set_extent(extent, crs=proj)

image-20220719223944436

add_geometries加载自定义地图
ax.add_geometries(shapereader.Reader(china_map).geometries(), crs=proj, facecolor='none',edgecolor='k',linewidth=1)
ax.set_extent(extent, crs=proj)

image-20220719223944436

完整代码如下
import numpy as np
import cartopy.crs as ccrs
import matplotlib.pyplot as plt
import cartopy.feature as cfeat
import cartopy.io.shapereader as shapereader

shp_path = "D:/data2/china/china/province_9south.shp"#切记路径一定要英文,血的教训!!!

fig  = plt.figure(figsize=(6,6)) #创建画布
proj = ccrs.PlateCarree()#创建投影,选择cartopy的默认投影
ax   = fig.subplots(1, 1, subplot_kw={'projection': proj})  #子图
extent = [70,140,2,60]#限定绘图范围

china_map = cfeat.ShapelyFeature(shapereader.Reader(shp_path).geometries(), proj, edgecolor='k', facecolor='none')
ax.add_feature(china_map, linewidth=1)#添加读取的数据,设置线宽
ax.set_extent(extent, crs=proj)

#ax.add_geometries(shapereader.Reader(shp_path).geometries(), crs=proj, #facecolor='none',edgecolor='k',linewidth=1)
#ax.set_extent(extent, crs=proj)
方法二

使用GMT 中文社区上提供的省界的经纬度数据进行绘图。GMT 中文社区为用户提供了一套名为CN-border的基本准确的数字化中国国界,、省界数据。在此向 GMT 中文社区的维护和贡献者表示感谢!

但需要注意的是,即便正确绘制了国界、省界,所绘地图如果要在境内公开展示,依然需要审图。个人没有提交审图申请的资格,需要付费,委托给有资质的企事业单位代为提交审图申请。下为GMT 中文社区网址,大家可自行下载

https://docs.gmt-china.org/latest/dataset-CN/CN-border/

CN-border 数据提供了三个数据文件:

  • CN-border-La.gmt:中国国界、省界、十段线以及南海诸岛数据
  • CN-border-L1.gmt:中国国界、十段线以及南海诸岛数据,不含省界数据
  • ten-dash-line.gmt:仅十段线数据
import numpy as np
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cfeature

#加载边界数据,CN-border-La.dat下载
#https://github.com/gmt-china/china-geospatial-data/releases
with open('C:/Users/啊啊啊啊/Desktop/china-geospatial-data-GB2312/CN-border-La.gmt') as src:
    context = src.read()
    blocks = [cnt for cnt in context.split('>') if len(cnt) > 0]
    borders = [np.fromstring(block, dtype=float, sep=' ') for block in blocks]

# 画布
fig = plt.figure(figsize=[10, 8])
# 设置投影并绘制主图形
ax = plt.axes(projection=ccrs.LambertConformal(central_latitude=90,
                                               central_longitude=105))
# 绘制线条
for line in borders:
    ax.plot(line[0::2], line[1::2], '-', lw=1, color='k',
            transform=ccrs.Geodetic())
# 绘制边界线
for line in borders:
    sub_ax.plot(line[0::2], line[1::2], '-', lw=1, color='k',
                transform=ccrs.Geodetic())
# 设置图形范围
sub_ax.set_extent([105, 125, 0, 25])
# 显示图形
plt.show()

image-20220720144752772

方法三

cnmaps绘图,cnmaps是Cartopy的一个扩展包,内置了中国地图数据源,该数据源自高德地图API,大家可以自行去该项目进行学习

https://github.com/Clarmy/cnmaps

重点区域突出显示

该示例是在前文方法一的基础上进行绘制

Cartopy 默认对所有地物使用相同的外观,如果需要突出显示某些地物,就必须进行筛选。这里从province_9south.shp这份全国省级行政区划数据中选中山西(通过属性表SHENG_ID字段),然后使用 ax.add_geometries() 方法将它们加入到原有地图元素中。

for state, geos in zip(shapereader.Reader(shp_path).records(), shapereader.Reader(shp_path).geometries()) :
    if state.attributes['SHENG_ID'] in [14 ]:#如果'14'在'SHENG_ID'这个字段中
        ax.add_geometries([geos], crs=proj,facecolor='None', edgecolor='blue')  # 重点区域上色

[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-f6fujWKR-1658307845790)(../AppData/Roaming/Typora/typora-user-images/image-20220720114405809.png)]

关于属性表,你可以通过gdal来打印出该shpfile文件所具有的属性表,也可以通过arcgis,qgis等GIS软件打开查询

仅绘制山西省

注意在前文中我加粗的那个词“加入”这意味着上文中重点区域突出显示的山西省图层是经过筛选而添加到地图中的,而不是对山西省这个范围进行符号化,如果我们要只显示山西省的话,就将上述代码中显示全国范围的代码注释掉即可,参见下方示例代码

fig = plt.figure(figsize=(6,6))  
ax = fig.subplots(1, 1, subplot_kw={'projection': proj})  
#ax.add_geometries(shapereader.Reader(shp_path).geometries(), crs=proj,facecolor='none',edgecolor='k',linewidth=1)
ax.set_extent(extent, crs=proj)

for state, geos in zip(shapereader.Reader(shp_path).records(), shapereader.Reader(shp_path).geometries()) :
    if state.attributes['SHENG_ID'] in [14]:
        ax.add_geometries([geos], crs=proj,facecolor='None', edgecolor='red')  # 重点区域上色

image-20220720120049159

添加南海小地图

在之前的的学习中我们知道,cartopy绘制的地图称为子图,在绘制中国地图时候,有时候由于地图大小的限制,我们无法展示部分地区如南海,常规的方法是绘制两幅地图,比如一张为全国地图,一张为局部地图,也就是常说的南海小地图。

常见的subplot和subplot2grid函数一般来说绘制的地图大小是一样的,不容易展示比例大小,所以我们选择add_axes()命令来绘制两个大小不一样的子图。需要注意的是在绘制时我们需要定义两个extent参数,即分别为总的地图和南海小地图分别定义画布大小绘图范围

extent = [105, 133, 15, 45]#限定绘图范围

fig = plt.figure(figsize=(6,8))  #注意画布大小
ax = fig.subplots(1, 1, subplot_kw={'projection': proj})  
ax.add_geometries(shapereader.Reader(shp_path).geometries(), crs=proj,facecolor='none',edgecolor='k',linewidth=1)
ax.set_extent(extent, crs=proj)

left, bottom, width, height = 0.67, 0.16, 0.23, 0.27#南海小地图大小
ax2 = fig.add_axes([left, bottom, width, height], projection=proj)
ax2.add_geometries(shapereader.Reader(shp_path).geometries(), crs=proj, facecolor='none',edgecolor='k',linewidth=1)
ax2.set_extent([105, 125, 0, 25])#注意南海小地图的范围

image-20220720122112863

添加地形图背景

一个简单的尝试

在看完了上述内容后还是觉得有些不足,让我们结合来Cartopy绘图入门中所讲述的内容来给他添加一个背景吧

import numpy as np
import cartopy.crs as ccrs
import matplotlib.pyplot as plt
import cartopy.feature as cfeature
import cartopy.io.shapereader as shapereader

shp_path = "D:/data2/china/china/province_9south.shp"#切记路径一定要英文,血的教训!!!

fig = plt.figure(figsize=(6,8))  #注意画布大小
proj = ccrs.PlateCarree()#创建投影,选择cartopy的默认投影
ax   = fig.subplots(1, 1, subplot_kw={'projection': proj})  #子图
extent = [105, 133, 15, 45]#限定绘图范围

ax = fig.subplots(1, 1, subplot_kw={'projection': proj})  
ax.add_geometries(shapereader.Reader(shp_path).geometries(), crs=proj,facecolor='none',edgecolor='k',linewidth=1)
ax.set_extent(extent, crs=proj)
#添加其他要素
ax.add_feature(cfeature.OCEAN.with_scale('50m'))
ax.add_feature(cfeature.LAND.with_scale('50m'))
ax.add_feature(cfeature.RIVERS.with_scale('50m'))
ax.add_feature(cfeature.LAKES.with_scale('50m'))

left, bottom, width, height = 0.67, 0.16, 0.23, 0.27#南海小地图大小
ax2 = fig.add_axes([left, bottom, width, height], projection=proj)
ax2.add_geometries(shapereader.Reader(shp_path).geometries(), crs=proj, facecolor='none',edgecolor='k',linewidth=1)
ax2.set_extent([105, 125, 0, 25])#注意南海小地图的范围

image-20220720145541535

这时候我们会疑惑,为什么南海小地图还是这样。在上文中提到过,南海小地图是另一幅地图,他与原来的中国地图你可以理解为是两幅地图。所以也需要在南海小地图中添加要素,也就是ax2中

import numpy as np
import cartopy.crs as ccrs
import matplotlib.pyplot as plt
import cartopy.feature as cfeature
import cartopy.io.shapereader as shapereader

shp_path = "D:/data2/china/china/province_9south.shp"#切记路径一定要英文,血的教训!!!

fig = plt.figure(figsize=(6,8))  #注意画布大小
proj = ccrs.PlateCarree()#创建投影,选择cartopy的默认投影
ax   = fig.subplots(1, 1, subplot_kw={'projection': proj})  #子图
extent = [105, 133, 15, 45]#限定绘图范围

ax = fig.subplots(1, 1, subplot_kw={'projection': proj})  
ax.add_geometries(shapereader.Reader(shp_path).geometries(), crs=proj,facecolor='none',edgecolor='k',linewidth=1)
ax.set_extent(extent, crs=proj)
#添加其他要素_中国地图
ax.add_feature(cfeature.OCEAN.with_scale('50m'))
ax.add_feature(cfeature.LAND.with_scale('50m'))
ax.add_feature(cfeature.RIVERS.with_scale('50m'))
ax.add_feature(cfeature.LAKES.with_scale('50m'))

left, bottom, width, height = 0.67, 0.16, 0.23, 0.27#南海小地图大小
ax2 = fig.add_axes([left, bottom, width, height], projection=proj)
ax2.add_geometries(shapereader.Reader(shp_path).geometries(), crs=proj, facecolor='none',edgecolor='k',linewidth=1)
ax2.set_extent([105, 125, 0, 25])#注意南海小地图的范围
#添加其他要素_南海小地图
ax2.add_feature(cfeature.OCEAN.with_scale('50m'))
ax2.add_feature(cfeature.LAND.with_scale('50m'))
ax2.add_feature(cfeature.RIVERS.with_scale('50m'))
ax2.add_feature(cfeature.LAKES.with_scale('50m'))

image-20220720150307392

添加图片作为背景
方法一

使用cartopy的add_image方法,这里以添加osm在线地图的图片为例子

imagery = OSM()
ax.add_image(imagery, 4)#四为影像级别

image-20220720162919914

完整代码如下

import cartopy.crs as ccrs
import cartopy.feature as cfeat
from cartopy.io.shapereader import Reader
import matplotlib.pyplot as plt
from matplotlib.image import imread
from cartopy.io.img_tiles import OSM

def create_map():
    # --设置shp路径,数据集已公开
    shp_path = 'D:/data2/china/china/province_9south.shp'
    # --设置tif路径,数据集已公开
    tif_path = 'D:/data2/china/china/NE1_50M_SR_W.tif'
    # --创建画图空间
    proj = ccrs.PlateCarree()  # 创建坐标系
    fig = plt.figure(figsize=(6, 8))  # 创建页面
    ax = fig.subplots(1, 1, subplot_kw={'projection': proj})  

    # --设置地图属性
    provinces = cfeat.ShapelyFeature(Reader(shp_path).geometries(),proj, edgecolor='k',facecolor='none')
    # --加载省界线
    ax.add_feature(provinces, lw=0.6, zorder=2)
    ax.set_extent([105, 133, 15, 45])  #可根据需求自行定义
   #添加osm图片
    imagery = OSM()
    ax.add_image(imagery, 4)#四为影像级别
ax = create_map()
plt.show()
方法二

添加图片所使用的第二种方法是 ax.imshow,是基于matplotlib的,这里以加载地形图图片为例

ax.imshow(imread(tif_path),extent=[-180, 180, -90, 90])注意范围

image-20220720154440290

绘制的完整代码如下:

import cartopy.crs as ccrs
import cartopy.feature as cfeat
from cartopy.io.shapereader import Reader
import matplotlib.pyplot as plt
from matplotlib.image import imread

def create_map():
    # --设置shp路径,数据集已公开
    shp_path = 'D:/data2/china/china/province_9south.shp'
    # --设置tif路径,数据集已公开
    tif_path = 'D:/data2/china/china/NE1_50M_SR_W.tif'
    # --创建画图空间
    proj = ccrs.PlateCarree()  # 创建坐标系
    fig = plt.figure(figsize=(6, 8))  # 创建页面
    ax = fig.subplots(1, 1, subplot_kw={'projection': proj})  

    # --设置地图属性
    provinces = cfeat.ShapelyFeature(Reader(shp_path).geometries(),proj, edgecolor='k',facecolor='none')
    # --加载省界线
    ax.add_feature(provinces, lw=0.6, zorder=2)
    ax.set_extent([105, 133, 15, 45])  #可根据需求自行定义
    # ax.stock_img()  ##可直接使用低分辨率背景
    # --加载高分辨率地形
    ax.imshow(imread(tif_path),extent=[-180, 180, -90, 90])
ax = create_map()
plt.show()

你可以试试结合上述所讲的内容,做出一副这样的地图

image-20220720154648087

示例代码如下:

import cartopy.crs as ccrs
import cartopy.feature as cfeat
from cartopy.io.shapereader import Reader
import matplotlib.pyplot as plt
from matplotlib.image import imread

def create_map():
    # --设置shp路径,数据集已公开
    shp_path = 'D:/data2/china/china/province_9south.shp'
    # --设置tif路径,数据集已公开
    tif_path = 'D:/data2/china/china/NE1_50M_SR_W.tif'
    # --创建画图空间
    proj = ccrs.PlateCarree()  # 创建坐标系
    fig = plt.figure(figsize=(6, 8))  # 创建页面
    ax = fig.subplots(1, 1, subplot_kw={'projection': proj})  

    # --设置地图属性
    provinces = cfeat.ShapelyFeature(Reader(shp_path).geometries(),proj, edgecolor='k',facecolor='none')
    # --加载省界线
    ax.add_feature(provinces, lw=0.6, zorder=2)
    ax.set_extent([105, 133, 15, 45])  #可根据需求自行定义
    # ax.stock_img()  ##可直接使用低分辨率背景
    # --加载高分辨率地形
    ax.imshow(imread(tif_path),origin='upper', transform=proj,extent=[-180, 180, -90, 90])
    # --设置网格点属性
    gl = ax.gridlines(crs=ccrs.PlateCarree(),draw_labels=True,linewidth=1.2,color='k',alpha=0.5,linestyle='--')
    # --设置南海子图
    left, bottom, width, height = 0.67, 0.16, 0.23, 0.27
    ax2 = fig.add_axes([left, bottom, width, height], projection=proj)
    ax2.add_feature(provinces, linewidth=0.6, zorder=2)
    ax2.set_extent([105, 125, 0, 25])
    # ax2.stock_img()  ##可直接使用低分辨率背景
    # --加载高分辨率地形
    ax2.imshow(imread(tif_path),origin='upper',transform=proj,extent=[-180, 180, -90, 90])
    return ax
ax = create_map()
plt.show()

注:本例教程参考至和鲸大佬【大自然搬运工】,针对大佬的代码做了一些简化方便大家理解;大佬绘制经纬网好像使用的是matplotlib库,与我之前所讲述的cartopy库绘制不一样。原文章如下:

https://www.heywhale.com/mw/project/5f3c95a3af3980002cbf560b

学习心得

本文作为基础入门教程就写到这里吧,这几天学习强度有点大了。在本文的学习过程中,我在和鲸社区找到了大量优质的学习博文,十分建议大家可以去看一看,同时以也才cartopy的文档里得到了非常多的帮助,里面还有很多地图的绘制方式,如果有机会,我希望我可以学习一下。

和鲸社区
https://www.heywhale.com/home
cartopyAPI参考
https://scitools.org.uk/cartopy/docs/latest/reference/generated/cartopy.mpl.geoaxes.GeoAxes.html#cartopy.mpl.geoaxes.GeoAxes.add_feature
cartopy更多地图
https://scitools.org.uk/cartopy/docs/latest/gallery/index.html

数据下载

行政区划_省界.shp

https://wwc.lanzouw.com/iQ4TT0843gpc

GMT 中文社区提供的数据

https://wwc.lanzouw.com/iYLrc0843gvi

地形图图片

https://wwc.lanzouw.com/ipPrM0843g8f

参考文档

https://zhajiman.github.io/post/cartopy_introduction/
https://mp.weixin.qq.com/s/FKNBOD2Yo4sawWTO76zZgA
https://mp.weixin.qq.com/s/m4Ao0--o1umSICEWO9UlKw
https://mp.weixin.qq.com/s/oe2ncVdqfdjaPSDZ2-YLQg
https://cloud.tencent.com/developer/article/1790590
https://mp.weixin.qq.com/s/wbyrkmVoaWgz6MdBI_6BZg
https://mp.weixin.qq.com/s/5RiuM2AQxNIvgNcLf5fCDw
https://mp.weixin.qq.com/s/FGCPIoiC5OUkGGf_IreWZQ
https://www.cnblogs.com/youxiaogang/p/14247184.html
https://www.heywhale.com/mw/project/629eefe52d8168866e54adc9
https://www.heywhale.com/mw/project/5f1a4df794d484002d2db14a
https://www.heywhale.com/mw/project/5f3c95a3af3980002cbf560b
https://blog.csdn.net/m0_37362454/article/details/81511427
https://www.osgeo.cn/pygis/cartopy-feature.html
https://vimsky.com/examples/detail/python-module-cartopy.mpl.ticker.html
https://scitools.org.uk/cartopy/docs/v0.14/matplotlib/feature_interface.html
https://mp.weixin.qq.com/s/jpZOpnFvMwi4ZTafflXIzw
https://gnss.help/2018/04/24/cartopy-gallery/index.html
https://blog.csdn.net/weixin_42372313/article/details/113572786
https://scitools.org.uk/cartopy/docs/latest/index.html
https://scitools.org.uk/cartopy/docs/v0.13/matplotlib/geoaxes.html
https://waterhackweek.github.io/visualization/04-geopandas-and-cartopy/
https://scitools.org.uk/cartopy/docs/v0.14/matplotlib/feature_interface.html
https://mp.weixin.qq.com/s/ASV-VI6gfbea90vtJbdpIw
Logo

为开发者提供学习成长、分享交流、生态实践、资源工具等服务,帮助开发者快速成长。

更多推荐