赞
踩
绘图效果如图所示:
- import cartopy.crs as ccrs
- from cartopy.io.shapereader import Reader
- import matplotlib.pyplot as plt
- import rioxarray as rxr
- import numpy as np
- import shapefile
- from matplotlib.path import Path
- from matplotlib.patches import PathPatch
- from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter
-
- '''
- 绘图数据准备:使用ArcGIS查看栅格数据和shp文件的坐标系信息,使之保持一致,才能保证数据的位置信息和shp重合,正确绘图
- 绘图思路:先绘制栅格数据的填色图,再利用shp2clip进行白化处理,最后叠加shp边界
- created by sunny_xx.
- '''
-
-
- def shp2clip(originfig, ax, shpfile, region): # 用于白化选定区域以外的数据
- sf = shapefile.Reader(shpfile) # 若使用的shp中有中文报错,可以添加 ,encoding='gbk' 在shpfile之后
- vertices = []
- codes = []
- '''使用maskout中的函数shp2clip(originfig, ax, shpfile, region)必须使得region与shape_rec.record中的标识符对应
- 否则应该修改region中的参数或者修改if语句中shape_rec.record[]的索引使二者匹配'''
- for shape_rec in sf.shapeRecords():
- # print(shape_rec.record[0:]) # 查看shape_rec.record中的内容
- if shape_rec.record[1] in region: # 这里需要找到和region匹配的唯一标识符,record[]中必有一项是对应的。
- pts = shape_rec.shape.points
- prt = list(shape_rec.shape.parts) + [len(pts)]
- for i in range(len(prt) - 1):
- for j in range(prt[i], prt[i + 1]):
- vertices.append((pts[j][0], pts[j][1]))
- codes += [Path.MOVETO]
- codes += [Path.LINETO] * (prt[i + 1] - prt[i] - 2)
- codes += [Path.CLOSEPOLY]
- clip = Path(vertices, codes)
- clip = PathPatch(clip, transform=ax.transData)
- for contour in originfig.collections:
- contour.set_clip_path(clip)
- return clip
-
-
- def location_index(shp_records_generator, loc): # 找到指定绘图区域需要额外添加的市或县的边界
- n = 0
- loc_index = []
- for e in shp_records_generator:
- for f in loc:
- if f in e.attributes.get('市') or f in e.attributes.get('省'):
- loc_index.append(n)
- # print(e.attributes)
- n = n + 1
- return loc_index
-
-
- # 字体设置
- plt.rcParams['axes.unicode_minus'] = False # 显示负号
- config = {
- "font.family": 'serif', # 英文显示为Times New Roman字体
- "font.size": 12, # 字体大小
- "mathtext.fontset": 'stix', # 公式字体为stix
- "font.serif": ['SimSun'], # 中文显示为宋体
- }
- plt.rcParams.update(config)
- # 读取栅格数据和省市边界的shp文件
- tiff_file = r'***路径\**.tif' # 栅格数据存储路径
- shp_file = r'***路径\省.shp'
- shp_city = r"***路径\市.shp"
- shpFile = Reader(shp_file)
- shpCity = Reader(shp_city)
- TiffData = rxr.open_rasterio(tiff_file)
-
- # 指定需要显示的区域,注意名称的准确性,可通过函数shp2clip中语句 print(shape_rec.record[0:])查看列表city中的名称是否正确
- city = ['北京市', '天津市', '河北省', '山西省']
- recCity = shpCity.records()
- City_Index = location_index(recCity, city) # 获取city中所有城市的索引
- shpFile_list = list(shpFile.geometries())
- shpCity_list = list(shpCity.geometries())
- # 读取栅格数据的经纬度网格
- lon, lat = np.meshgrid(TiffData['x'], TiffData['y'])
- # 读取栅格数据的浓度值
- data = TiffData[0]
- '''
- 若需要计算浓度差,则可以按相同方法读取另一栅格数据,直接做差即可
- data0 = ds0[0]
- data_diff = data - data0
- '''
- # 设置绘图经纬度范围,可自行修改绘图区域
- lon1 = min(TiffData['x'])
- lon2 = 125
- lat1 = min(TiffData['y'])
- lat2 = max(TiffData['y'])
- extent = [lon1, lon2, lat1, lat2]
- # 绘图
- proj = ccrs.PlateCarree()
- fig, ax = plt.subplots(1, 1, dpi=100, subplot_kw={'projection': proj})
- ax.set_extent(extent, crs=proj)
- ax.xaxis.set_major_formatter(LongitudeFormatter()) # 将横坐标转换为经度格式
- ax.yaxis.set_major_formatter(LatitudeFormatter()) # 将纵坐标转换为纬度格式
- ax.set_xticks(np.arange(int(lon1), int(lon2), 4)) # 设置需要显示的经度刻度
- ax.set_yticks(np.arange(int(lat1), int(lat2), 3))
- # 绘制填色图(cmap=plt.cm.GnBu为颜色设置,GnBu代表从绿色到蓝色的渐变色;levels=np.linspace(0, 400, 401)指bar的最小值最大值)
- cf = ax.contourf(lon, lat, data, transform=proj, cmap=plt.cm.GnBu, levels=np.linspace(0, 400, 401))
- # 根据指定的范围进行裁剪(city中包含的省市名称)
- clip = shp2clip(cf, ax, shp_file, city)
- # 设置colorbar
- cb = fig.colorbar(cf, ax=ax, shrink=0.9, extendfrac='auto', extendrect=True, location='right', fraction=0.05, pad=0.05)
- lvs = np.arange(0, 400, 40)
- cb.set_ticks(lvs) # 设置色标刻度
- cb.set_ticklabels(lvs) # 设置色标刻度,如果想设置不同于xy轴labels的字体大小,使用, fontsize=8修改即可
- cb.set_label('PM$_{2.5}$ ($\mu$g/m$^3$)')
- # 叠加绘制shp边界
- Map = ax.add_geometries(shpFile_list, proj, edgecolor='grey', facecolor='none', lw=0.6) # 叠加省界
- for m in City_Index:
- Map = ax.add_geometries(shpCity_list[m:m + 1], proj, edgecolor='k', facecolor='none', lw=0.6) # 叠加指定城市
- plt.savefig(r'**路径\test_PM2.5.png', dpi=400) # *********图片保存路径************
- plt.show()

数据准备阶段要使shp文件和tiff数据的坐标系一致,这部分在python中做了尝试,发现还是在arcgis中操作更简单快捷。可自行查找arcgis坐标系转换教程。
包括arcgis在内的全部详细步骤以及示例tiff数据和shp文件见附件:
Copyright © 2003-2013 www.wpsshop.cn 版权所有,并保留所有权利。