tiff（栅格）数据绘图及白化（python详细注释版）

绘图效果如图所示：

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文件见附件：

(9条消息) tiff数据python绘图叠加shp边界（包括代码、shp和arcgis坐标转换说明文档）资源-CSDN文库https://download.csdn.net/download/weixin_43304836/87685304