学习笔记5 | 地图绘制——Cartopy

目录

5.1 Cartopy简介

5.2 多种投影介绍

5.3 地图信息添加

5.3.1 地理特征添加

5.3.2 网格和刻度线设置

5.4 填色图结合地图绘制

---

5.1 Cartopy简介

Cartopy是开源免费的Python扩展包，致力于用最简单直观的方式生成地图。

主要功能：

1. 设置地图投影信息及投影转换
2. 添加地图地理信息
3. 绘制shp文件

Cartopy除了跟Matplotlib有良好的接口外，跟Xarray也可以完美兼容，这样我们在使用Xarray的函数绘制图像的时候也可以非常方便地添加地理信息。

#Cartopy和Matplotlib的连接
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
fig = plt.figure()
ax = fig.add_subplot() #ax类型是axesSubplot

proj = ccrs.PlateCarree()
fig = plt.figure()
ax = fig.add_subplot(projection=proj) #ax类型是GeoaxesSubplot
ax.coastlines() #画海岸线
ax.gridlines() #画网格线
ax.stock_img() #叠加地形背景

#Cartopy和Xarray
import xarray as xr
ds = xr.open_dataset('filename.nc')
t = ds['t2m']-273.15
t.isel(time=0).plot()

proj = ccrs.PlateCarree()
fig = plt.figure()
ax = fig.add_subplot(projection=proj)
ax.coastlines()
t.isel(time=0).plot(ax=ax)

5.2 多种投影介绍

1. 简易圆柱投影（Plate Carree）
   1. 适于索引地图以及用于绘制随经度变化的现象
2. 兰勃特投影（Lambert）
   1. 最适合用于对中纬度东西方向分布的大陆板块进行等角制图
   2. 通常将标准纬线置于在要制图区域顶部下方和底部上方六分之一纬度范围处
3. 墨卡托投影（Mercator）
   1. 等角圆柱地图投影
   2. 适用于绘制赤道附近地区（例如印度尼西亚和太平洋部分海域）的大比例地图

crs即坐标参考系统（Coordinate Reference Systems）之意

https://scitools.org.uk/cartopy/docs/v0.15/crs/projections.hmtl

import cartopy.crs as ccrs
import matplotlib.pyplot as plt
projections = [
    ccrs.PlateCarree(central_longitude=130.0),
    ccrs.Mercator(central_longitude=130.0),
    ccrs.LambertConformal(central_longitude=130.0),
    ccrs.Orthographic(central_longitude=130.0)  ]
for proj in projections:
    fig = plt.figure()
    ax = fig.add_subplot(projection=proj)
    ax.coastlines()
    ax.stock_img()
    ax.set_title(f'{type(proj)}')

5.3 地图信息添加

ATT！ Cartopy自带的中国地图不符合我国的地图标准

5.3.1 地理特征添加

import cartopy.crs as ccrs
import matplotlib.pyplot as plt
import cartopy.feature as cfeat
fig = plt.figure()
proj = ccrs.PlateCarree()
ax = fig.add_subplot(projection=proj)
ax.coastlines()
ax.add_feature(cfeat.RIVERS.with_scale('110m'))
ax.add_feature(cfeat.LAKES.with_scale('110m'))
ax.add_feature(cfeat.LAND.with_scale('110m'))
ax.add_feature(cfeat.OCEAN.with_scale('110m'))
#黄河流域没显示
#三种分辨率 10m 50m 110m

for res in ['110m','50m','10m']: #50m效果比较好
    proj = ccrs.PlateCarree()
    ax = fig.add_subplot(projection=proj)
    ax.coastlines()
    ax.add_feature(cfeat.RIVERS.with_scale(res))
    ax.add_feature(cfeat.LAKES.with_scale(res))
    ax.add_feature(cfeat.OCEAN.with_scale(res))
    ax.add_feature(cfeat.LAND.with_scale(res))
    ax.set_title(f'{res}')
    ax.set_extent([90,180,10,60])

5.3.2 网格和刻度线设置

fig = plt.figure()
proj = ccrs.PlateCarree()
ax = fig.add_subplot(projection=proj)
ax.coastlines()
ax.add_feature(cfeat.RIVERS.with_scale('110m'))
ax.add_feature(cfeat.LAKES.with_scale('110m'))
ax.add_feature(cfeat.LAND.with_scale('110m'))
ax.add_feature(cfeat.OCEAN.with_scale('110m'))

gl = ax.gridlines(
    linewidth=0.9,
    color='k',
    linestyle='--',
    alpha=0.5,
    draw_labels=False )
#gl.top_labels = gl.right_labels = False #刻度显示不完整

import numpy as np
ax.set_xticks(np.arange(-180,181,30))
ax.set_yticks(np.arange(-90,91,30))   #刻度完整但是不标注NSWE,而且网格不是30°间隔

from cartopy.mpl.ticker import LongitudeFormatter,LatitudeFormatter
ax.xaxis.set_major_formatter(LongitudeFormatter)
ax.yaxis.set_major_formatter(LatitudeFormatter)

ax.xaxis.set_minor_locator(plt.MultipleLocator(10))
ax.yaxis.set_minor_locator(plt.MultipleLocator(10))

5.4 填色图结合地图绘制

import xarray as xr
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cfeat

ds = xr.open_dataset('filename.nc')
t = ds['t2m']-273.15
fig = plt.figure(figsize = (12,12))
proj = ccrs.PlateCarree() #这个是底图投影
ax = fig.add_subplot(projection=proj)

# --- 设置特征
country = cfeat.ShapelyFeature(
    gpd.read_file('filename.shp')['geometry'],
    proj,
    edgecolor = 'k',
    facecoloe = None
)
ax.add_feature(country)
province = cfeat.ShapelyFeature(
    gpd.read_file('filename.shp')['geometry'],
    proj,
    edgecolor = 'k',
    facecoloe = None
)
ax.add_feature(province)

ax.add_feature(cfeat.RIVERS)
ax.add_feature(cfeat.OCEAN)
ax.add_feature(cfeat.LAKES)

# --- 设置经纬网格
gl = ax.gridlines(
    linewidth=0.9,
    color='k',
    linestyle='--',
    alpha=0.5,
    draw_labels=False )
#gl.top_labels = gl.right_labels = False #刻度显示不完整

import numpy as np
ax.set_xticks(np.arange(-180,181,5))
ax.set_yticks(np.arange(-90,91,5))   #刻度完整但是不标注NSWE,而且网格不是30°间隔

from cartopy.mpl.ticker import LongitudeFormatter,LatitudeFormatter
ax.xaxis.set_major_formatter(LongitudeFormatter)
ax.yaxis.set_major_formatter(LatitudeFormatter)

ax.xaxis.set_minor_locator(plt.MultipleLocator(1))
	ax.yaxis.set_minor_locator(plt.MultipleLocator(1))
# --- 设置区域
ax.set_extent([100,140,15,50],crs=proj)
# --- 可视化
cf = t.isel(time=0).plot(
    ax=ax,
    transform=ccrs.PlateCarree(), #transform是数据的坐标系统
    add_colorbar=False )
# --- 设置色标
fig.colorbar(cf,ax=ax,orientation='horizontal',pad=0.06,aspect=50)

import geopandas as gpd
china = gdp.read_file('filename.shp')
#修改在最上边

# --- 设置南海子图
ax1 = fig.add_axes(
    [0.657,0.284,0.2,0.24],
    projection = proj )

country = cfeat.ShapelyFeature(
    gpd.read_file('filename.shp')['geometry'],
    proj,
    edgecolor = 'k',
    facecoloe = None
)
ax1.add_feature(country)
province = cfeat.ShapelyFeature(
    gpd.read_file('filename.shp')['geometry'],
    proj,
    edgecolor = 'k',
    facecoloe = None
)
ax1.add_feature(province)
ax1.set_extent([105,125,0,25],crs=proj)
cf = t.isel(time=0).plot(
    ax=ax1,
    transform=ccrs.PlateCarree(), #transform是数据的坐标系统
    add_colorbar=False )
ax1.set_title('')