根据点矢量删除栅格中的对应像元（主要操作为矢量转栅格）

之前写过一篇文档，是利用ArcGIS完成的本文的操作，但过程十分繁琐，其中主要的一个过程在于点矢量转栅格的过程，因此受文章python实现使用GDAL实现矢量转栅格的启发，利用gdal.RasterizeLayer进行矢量转栅格的操作后，实现本文的要求。
代码如下所示：

from osgeo import gdal, gdalconst
from osgeo import ogr
import gdalconst
import os
from os.path import join
import numpy as np

raster_root = r'D:\file\novel coronavirus\Infrared remote sensing\Annual output of steel\TangShan\Temperature Cut'  # 原影像
shp_root = r'D:\file\novel coronavirus\Infrared remote sensing\Annual output of steel\TangShan\Dian to shp 1'  # 裁剪矩形
save_root = r'D:\file\novel coronavirus\Infrared remote sensing\Annual output of steel\TangShan\Dian to shp result'

raster_file = []
shp_file = []
raster_names = sorted(os.listdir(raster_root))
# 得到所需的shp文件及tif栅格图像
for raster_name in raster_names:
    raster_nameend = os.path.splitext(raster_name)[-1]
    if (raster_nameend == '.tif'):
        raster_file.append(raster_name)
shp_names = sorted(os.listdir(shp_root))
for shp_name in shp_names:
    shp_nameend = os.path.splitext(shp_name)[-1]
    if (shp_nameend == '.shp'):
        shp_file.append(shp_name)

# 理论上一个矢量对应一张栅格，因此用栅格的数量来做循环的判断
for i in range(len(raster_file)):
    raster_for = os.path.splitext(raster_file[i])[0]
    shp_for = os.path.splitext(shp_file[i])[0]
    save_name = shp_for + 'Dian.tif'
    save_path = join(save_root, save_name)

    # 得到栅格文件的最终路径
    rasterFile = join(raster_root, raster_file[i])
    shpFile = join(shp_root, shp_file[i])
    dataset = gdal.Open(rasterFile, gdalconst.GA_ReadOnly)
    raster_data = dataset.GetRasterBand(1).ReadAsArray().astype(np.float32)

    # GetGeoTransform：得到左上角横坐标、纵坐标，水平空间及垂直空间分辨率等内容
    geo_transform = dataset.GetGeoTransform()
    cols = dataset.RasterXSize  # 列数
    rows = dataset.RasterYSize  # 行数

    x_min = geo_transform[0]
    y_min = geo_transform[3]
    pixel_width = geo_transform[1]

    # 打开矢量文件
    shp = ogr.Open(shpFile, 0)
    # GetLayerByIndex：获取图层个数，shp文件一般为一层
    m_layer = shp.GetLayerByIndex(0)
    # 实现矢量转换为栅格
    target_ds = gdal.GetDriverByName('GTiff').Create(save_path, xsize=cols, ysize=rows, bands=1,
                                                     eType=gdal.GDT_Float32)
    target_ds.SetGeoTransform(geo_transform)
    target_ds.SetProjection(dataset.GetProjection())

    band = target_ds.GetRasterBand(1)
    band.SetNoDataValue(0)
    band.FlushCache()
    
    # gdal.RasterizeLayer：根据shp或geojson矢量创建栅格
    # 跟shp字段给栅格像元赋值（本次试验中是RASTERVALU字段的值）
    gdal.RasterizeLayer(target_ds, [1], m_layer, options=["ATTRIBUTE=RASTERVALU"])
    # gdal.RasterizeLayer(target_ds, [1], m_layer) # 多边形内像元值的全是255
    del dataset
    del target_ds
    shp.Release()

    # 重新打开矢量转栅格的结果
    ShpToDian_file = gdal.Open(save_path, gdalconst.GA_ReadOnly)
    ShpToDian_data = ShpToDian_file.GetRasterBand(1).ReadAsArray().astype(np.float32)

    # 用原始栅格图像减去矢量转栅格的结果，以此来删除目标像元
    differance = raster_data - ShpToDian_data
    differance = differance.astype(np.float32)
    differance[differance == 0] = np.nan

    # 删除结果保存的文件名及路径设置
    differance_savename = 'R' + shp_for.split('H')[0] + 'T.tif'
    differance_saveroot = r'D:\file\novel coronavirus\Infrared remote sensing\Annual output of steel\TangShan\Temperature Result'
    differance_savepath = join(differance_saveroot, differance_savename)

    # 输出数据
    # 设置坐标系等
    gtiff_driver = gdal.GetDriverByName("GTiff")
    out_ds = gtiff_driver.Create(differance_savepath,
                                 differance.shape[1], differance.shape[0], 1, gdal.GDT_Float32)
    # 设置输出数据坐标投影为原始影像的坐标投影
    out_ds.SetProjection(ShpToDian_file.GetProjection())
    out_ds.SetGeoTransform(ShpToDian_file.GetGeoTransform())
    out_band = out_ds.GetRasterBand(1)
    out_band.WriteArray(differance)
    out_band.FlushCache()