gdal矢量裁剪栅格（遥感影像）

裁剪（矢量裁剪栅格）

GDAL OGR

#矢量批量裁剪栅格影像
#@liuxy


from osgeo import gdal
import os

input_file = r'C:\Users\BNU\Desktop\temp\HNDK\GY_2020\TIFF' #栅格数据（可以是多个栅格）
output_file = r'C:\Users\BNU\Desktop\temp\HNDK\GY_2020\clip/' #裁剪后存放的文件夹
input_shape = r'C:\Users\BNU\Desktop\temp\HNDK\GY_2020\TIFF\TEST1.shp'  #矢量数据
prefix = '.tif'

file_all = os.listdir(input_file) 
# print(file_all)
for file_i in file_all:
    if file_i.endswith(prefix):
        file_name = input_file + "\\"+file_i

        data = gdal.Open(file_name)
        ds = gdal.Warp(output_file + os.path.splitext(file_i)[0] + '_Clip1.tif',
              data,
              format = 'GTiff',
              cutlineDSName = input_shape,      
            #   cutlineWhere= "TBLXMC = '农作物'",#特定字段筛选
              cropToCutline=True, # 保证裁剪后影像大小跟矢量文件的图框大小一致
              outputType=gdal.GDT_Int16,
            #   creationOptions=["TILED=YES", "COMPRESS=LZW"],
              dstNodata = 0)

效果

裁剪前（栅格+矢量）

裁剪后

遇到问题，解决问题

当进行耕地矢量进行裁剪时，如果矢量面存在拓扑错误，相交等出现如下报错

Warning1: RingSelf-intersectionatornearpoint112 .4866642030000334 .830899357000078 

ERROR1: Cutlinepolygonisinvalid.

参考 https://blog.csdn.net/qq_45723511 。这里给出一种使用python+ogr函数库的解决方法。

#矢量批量裁剪栅格影像
#@liuxy

from osgeo import gdal,ogr
import os

# # # 注册所有的驱动
ogr.RegisterAll()
gdal.SetConfigOption("GDAL_FILENAME_IS_UTF8", "YES")
gdal.SetConfigOption("SHAPE_ENCODING", "GBK")

#------------------------修改拓扑错误（相交）---------------------#
def s_shp_si(shpFile):
    out_shp = shpFile[:-4] + '_ssi' + '.shp'
    # print(out_shp)
    # 打开数据
    ds = ogr.Open(shpFile, 0)
    if ds is None:
        print("打开文件 %s 失败！" % shpFile)
        return
    print("打开文件%s成功" % shpFile)
    # 获取该数据源中的图层个数，一般shp数据图层只有一个，如果是mdb、dxf等图层就会有多个
    m_layer_count = ds.GetLayerCount()
    m_layer = ds.GetLayerByIndex(0)
    if m_layer is None:
        print("获取第%d个图层失败\n", 0)
        return
    # 创建输出文件
    driver = ogr.GetDriverByName('ESRI Shapefile')
    if os.path.exists(out_shp):
        driver.DeleteDataSource(out_shp)
    outds = driver.CreateDataSource(out_shp)
    outlayer = outds.CreateLayer(out_shp[:-4], geom_type=ogr.wkbPolygon, options=["ENCODING=GBK"])
    # 获取输出层的要素定义
    outLayerDefn = outlayer.GetLayerDefn()
    # 对图层进行初始化，如果对图层进行了过滤操作，执行这句后，之前的过滤全部清空
    m_layer.ResetReading()
    # 获取投影
    prosrs = m_layer.GetSpatialRef()
    # 添加字段
    inLayerDefn = m_layer.GetLayerDefn()
    for i in range(0, inLayerDefn.GetFieldCount()):
        fieldDefn = inLayerDefn.GetFieldDefn(i)
        outlayer.CreateField(fieldDefn)

    # loop through the input features
    m_feature = m_layer.GetNextFeature()
    while m_feature:
        # print(m_feature)
        o_geometry = m_feature.GetGeometryRef()
        # 关键，合并几何
        o_geometry = o_geometry.Union(o_geometry)
        outfeature = ogr.Feature(outLayerDefn)
        outfeature.SetGeometry(o_geometry)
        # 遍历每个要素的字段，并设置字段属性
        for i in range(0, outLayerDefn.GetFieldCount()):
            print(outLayerDefn.GetFieldDefn(i).GetNameRef())
            outfeature.SetField(outLayerDefn.GetFieldDefn(i).GetNameRef(), m_feature.GetField(i))
        outlayer.CreateFeature(outfeature)
        # dereference the features and get the next input feature
        outfeature = None
        m_feature = m_layer.GetNextFeature()

    ds.Destroy()
    outds.Destroy()
    # 保存投影文件

    prjfile = (out_shp[:-4] + '.prj')
    fn = open(prjfile, 'w')
    fn.write(prosrs.ExportToWkt())
    fn.close()

    return out_shp


#-----------------路径文件-------------
input_file = r'C:\Users\BNU\Desktop\temp\HNDK\GY_2020\TIFF' #栅格数据（可以是多个栅格）
output_file = r'C:\Users\BNU\Desktop\temp\HNDK\GY_2020\clip/' #裁剪后存放的文件夹
# input_shape = r'C:\Users\BNU\Desktop\temp\HNDK\GY_2020\TIFF\TEST1.shp'  #矢量数据（不存在拓扑错误）
input_shape = r'D:\矢量数据\shp\GY_SHP_CROP.shp'
prefix = '.tif'

## 矢量拓扑错误（相交）的改正
## 如果不存在拓扑错误（相交），可以将下面两行注销掉
input_shape = s_shp_si(input_shape) #经过了处理矢量会存在同级文件夹中
print(input_shape,'input_shape_correct')


file_all = os.listdir(input_file) 
# print(file_all)
for file_i in file_all:
    if file_i.endswith(prefix):
        file_name = input_file + "\\"+file_i

        data = gdal.Open(file_name)
        ds = gdal.Warp(output_file + os.path.splitext(file_i)[0] + '_Clip1.tif',
              data,
              format = 'GTiff',
              cutlineDSName = input_shape,      
            #   cutlineWhere= "TBLXMC = '农作物'",#特定字段筛选
              cropToCutline=True, # 保证裁剪后影像大小跟矢量文件的图框大小一致
              outputType=gdal.GDT_Int16,
            #   creationOptions=["TILED=YES", "COMPRESS=LZW"],
              dstNodata = 0)

参考链接

https://blog.csdn.net/qq_44935981/article/details/122141636
https://blog.csdn.net/weixin_42990464/article/details/117693553