ICode9

精准搜索请尝试: 精确搜索
首页 > 编程语言> 文章详细

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

2021-06-01 13:00:19  阅读:651  来源: 互联网

标签:01 lb geometries Python range ccrs ax linewidth 制图


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

文章目录


前言

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.')


三、制图结果

制图结果

标签:01,lb,geometries,Python,range,ccrs,ax,linewidth,制图
来源: https://blog.csdn.net/EWBA_GIS_RS_ER/article/details/117404469

本站声明: 1. iCode9 技术分享网(下文简称本站)提供的所有内容,仅供技术学习、探讨和分享;
2. 关于本站的所有留言、评论、转载及引用,纯属内容发起人的个人观点,与本站观点和立场无关;
3. 关于本站的所有言论和文字,纯属内容发起人的个人观点,与本站观点和立场无关;
4. 本站文章均是网友提供,不完全保证技术分享内容的完整性、准确性、时效性、风险性和版权归属;如您发现该文章侵犯了您的权益,可联系我们第一时间进行删除;
5. 本站为非盈利性的个人网站,所有内容不会用来进行牟利,也不会利用任何形式的广告来间接获益,纯粹是为了广大技术爱好者提供技术内容和技术思想的分享性交流网站。

专注分享技术,共同学习,共同进步。侵权联系[81616952@qq.com]

Copyright (C)ICode9.com, All Rights Reserved.

ICode9版权所有