Python-Cartopy制图学习01-中国区域SPEI空间制图

多做事，少说话，多运动，忌久坐，早点睡，少熬夜。

---

前言

1. 概述

• 此系列博文记录Python的Cartopy库的学习笔记
• 此篇博文绘制中国区域2010年5月SPEI03的空间分布
• Cartopy库的安装可以参考Cartopy 简介与安装。

2. 版本

2.1 2021年6月1日 Version1

---

一、绘图数据

• 2010年5月份的3个月尺度的SPEI
• 中国边界
• 中国省界
• 九段线

二、制图程序

# -*- coding: utf-8 -*-
"""
1. 程序目的
   (1) 绘制2010年5月份中国SPEI03的空间分布图
   
2. 山东青岛 2021年06月1日

3. 数据

4. 参考资料
   'https://scitools.org.uk/cartopy/docs/v0.13/matplotlib/advanced_plotting.html'

"""

# 0. 相关包的导入
import numpy as np
import cartopy.crs as ccrs

import matplotlib.colors as mcolors
from cartopy.io.shapereader import Reader

import matplotlib.pyplot as plt
import matplotlib.ticker as mticker
import cartopy.feature as cfeature
from matplotlib import rcParams
from osgeo import gdal

# 设置matplotlib的字体
config = {"font.family":'SimHei', "font.size": 16, "mathtext.fontset":'stix'}
rcParams.update(config)
rcParams['font.sans-serif'] = ['SimHei']  # 用来正常显示中文标签
rcParams['axes.unicode_minus'] = False  # 用来正常显示负号

# 1. 路径处理和基本变量定义
   
rootdir = r'C:\Users\Desktop\ChinaSPEI03\\'
shpname = rootdir + 'china\\china_country'
shpname2 = rootdir + 'china\\china_nine_dotted_line'
chinaprovince = rootdir + 'china\\china'

# 2. tif数据读取

tif_file = rootdir + 'China_201005SPEI03_Mask.tif'

in_ds = gdal.Open(tif_file)

# proj_wkt = in_ds.GetProjection()
# print(proj_wkt)

# proj = osr.SpatialReference()
# proj.ImportFromWkt(proj_wkt)
# print(proj)

 # 获取tif数据的基本信息，用于制图设置
geo_transform = in_ds.GetGeoTransform()
#print(geo_transform)

origin_x = geo_transform[0]
origin_y = geo_transform[3]
pixel_width = geo_transform[1]
pixel_height = geo_transform[5]

 # 行数和列数
n_rows = in_ds.RasterYSize
n_cols = in_ds.RasterXSize

#print(n_rows,n_cols)
     
in_band = in_ds.GetRasterBand(1) #打开波段1
in_data = in_band.ReadAsArray()

 # 数据范围
lon_range_new = np.linspace(origin_x,(origin_x+pixel_width*n_cols-pixel_width),n_cols)
lat_range_new = np.linspace(origin_y,(origin_y + pixel_height*n_rows-pixel_height),n_rows)
 
 # 数组无效值掩膜
mask = (in_data < -999)
in_data_mask = np.ma.array(in_data,mask=mask)
 
# 3. 绘图
fig = plt.figure(figsize=(6,4),dpi=600)  #创建figure对象

 # 创建画图空间, 使用兰伯特投影
   # 投影设置
proj = ccrs.LambertConformal(central_longitude=107.5, \
                                 central_latitude=36.0, standard_parallels=(25,47))
ax = fig.subplots(1, 1, subplot_kw = {'projection': proj})  #主图
   # 图名
ax.set_title('中国区域2010年5月SPEI03',
             fontsize = 12,
             loc='center')
  # 设置颜色显示范围
norm = mcolors.Normalize(vmin=-2.5,vmax=2.5)
 
 # 设置网格点属性，以下用于兰伯特投影
region = [80, 130, 15.5, 52.5]
ax.set_extent(region)
lb = ax.gridlines(draw_labels=None,x_inline=False, y_inline=False,\
     linewidth=0.1, color='gray', alpha=0.8, linestyle='--' )
lb.xlocator = mticker.FixedLocator(range(0,180,10))
lb.ylocator = mticker.FixedLocator(range(0,90,10))
lb = ax.gridlines(draw_labels=True,x_inline=False, y_inline=False, \
     linewidth=0.1, color='gray', alpha=0.9, linestyle='--' )
lb.top_labels = False
lb.right_labels = None
lb.xlocator = mticker.FixedLocator(range(90, 129, 10))
lb.ylocator = mticker.FixedLocator(range(20, 50, 10))
lb.ylabel_style = {'size': 10, 'color': 'k'}
lb.xlabel_style = {'size': 10, 'color': 'k' }
lb.rotate_labels = False

 
  # 画pcolormesh图
cs = ax.pcolormesh(lon_range_new, lat_range_new, in_data_mask,\
     cmap = 'Spectral', norm=norm,transform = ccrs.PlateCarree())   
 
  # 添加海岸线
ax.add_feature(cfeature.COASTLINE.with_scale('50m'),linewidth = 0.5)  
    
  # 绘制中国省界
ax.add_geometries(Reader(chinaprovince + '.shp').geometries(), ccrs.PlateCarree(), \
     facecolor = 'none', edgecolor = 'gray', linewidth = 0.8)   
    
  # 绘制中国边界
ax.add_geometries(Reader(shpname + '.shp').geometries(), ccrs.PlateCarree(), \
     facecolor = 'none', edgecolor = 'k', linewidth = 1)     
       
  # 绘制九段线
ax.add_geometries(Reader(shpname2 + '.shp').geometries(), ccrs.PlateCarree(), \
    facecolor = 'none', edgecolor = 'r', linewidth = 1)
    
  # 添加文字
fontdict = {'size':10,'color':'k','family':'Times New Roman'}
ax.text(67,47.5,'2010-05-SPEI03',
transform=ccrs.PlateCarree(),
fontdict=fontdict,
weight='normal')
 
  # 单独设置图例                
l = 0.21
b = 0.06 # 距离主图的位置:上下
h = 0.02
w = 1 - 2*0.2    

rect = [l,b,w,h] 
cbar_ax = fig.add_axes(rect)  

cb = plt.colorbar(cs, cax=cbar_ax,
                  orientation = 'horizontal',
                  label = 'drought intensity',
                  extend='both'
                  )

cb.ax.tick_params(labelsize=11)

x1_label = cb.ax.get_xticklabels() 
[x1_label_temp.set_fontname('Times New Roman') for x1_label_temp in x1_label]

font = {'family' : 'Times New Roman',
    'color'  : 'black',
    'weight' : 'normal',
    'size'   : 11.5,
    }
cb.set_label('Drought Degree' ,fontdict=font) #设置colorbar的标签字体及其大小


#############添加南海，实际上就是新建一个子图覆盖在之前子图的右下角##########
# 设置兰伯特投影
proj2 = ccrs.LambertConformal(central_longitude=115, \
                              central_latitude=12.5, standard_parallels=(3,20))
    

# 定位南海子图的位置，四个数字需要调整
# 四个数字分别是(left, bottom, width, height)：
# 子图左下角水平方向相对位置，左下角垂直方向相对位置，子图宽，子图高
ax2 = fig.add_axes([0.73, 0.11, 0.1, 0.25], projection = proj2)

ax2.set_extent([105, 125, 0, 26]) 

# 设置网格点
lb=ax2.gridlines(draw_labels=None,x_inline=False, y_inline=False,\
    linewidth=0.1, color='gray', alpha=0.8, linestyle='--' )
lb.xlocator = mticker.FixedLocator(range(90, 135, 5))
lb.ylocator = mticker.FixedLocator(range(0, 90, 5))

# 添加海岸线
ax2.add_feature(cfeature.COASTLINE.with_scale('50m'),linewidth = 0.5)

  #绘制中国省界
ax2.add_geometries(Reader(chinaprovince + '.shp').geometries(), ccrs.PlateCarree(), \
     facecolor = 'none', edgecolor = 'gray', linewidth = 0.8)   

# 添加中国边界
china = Reader(shpname + '.dbf').geometries()
ax2.add_geometries(china, ccrs.PlateCarree(), facecolor = 'none', \
                    edgecolor = 'black', zorder = 1)
    
# 添加数据
cs2 = ax2.pcolormesh(lon_range_new, lat_range_new, in_data_mask,\
     cmap = 'Spectral', norm=norm,transform = ccrs.PlateCarree())
    

# 九段线
ax2.add_geometries(Reader(shpname2 + '.shp').geometries(), ccrs.PlateCarree(), \
    facecolor = 'none', edgecolor = 'r', linewidth = 1)
    

# 存储制图结果
plt.savefig(rootdir+'ChinaDrought.jpg',bbox_inches='tight',dpi=300)    
       
print('Finished.')

---

三、制图结果