python gdal 按照属性提取矢量，并进行影像裁切

＃ -  *  -  coding：utf-8  -  *  - 
“”“ 
/ ******************* ************ 
ClipRaster.py:1 
    。按照属性提取矢量（在上层进行操作）
    2。矢量转为同源数据的栅格掩膜（某些矢量存在交叉点，无法使用gdalwarp进行裁切）
    3.栅格掩膜与源数据进行运行，得到裁切后的影像@version <1.1> 2018-07-03 Wujd：创建。“” 来自osgeo import ogr 来自osgeo import gdal 来自osgeo import gdal_array作为ga import os，math，ogr，osr ＃解决SHAPE文件的属性值










gdal.SetConfigOption('SHAPE_ENCODING', '')
# #解决中文路径乱码问题
# gdal.SetConfigOption("GDAL_FILENAME_IS_UTF8", "NO")

"""
路径设置
"""
input_vector ='./data/lns_210000_city.shp'
input_raster ='./data/GF1_WFV4_E124.8_N43.3_20180524_L1A0003210738_cal_orth.tiff'
extent_vector = './result/clip/lns_210000_city_42000102_clip.shp'
extent_raster = './result/clip/image_extent.tiff'
clip_raster = './result/clip/image_42000102_clip.tiff'


def copyInLayer(input_vector,extent_vector,filterStr):

    ds = ogr.Open(input_vector)
    driver = ogr.GetDriverByName("ESRI Shapefile")

    if os.access(extent_vector, os.F_OK):
        driver.DeleteDataSource(extent_vector)
    layer = ds.GetLayer()

    #findLayer = ds.ExecuteSQL("select * from lns_210000_city where 行政代码 = 211400 ")
    #filterStr="{} = {}".format('ID', 266)
    layer.SetAttributeFilter(filterStr)
    newds = driver.CreateDataSource（extent_vector）
    newds.CopyLayer（layer，''）
    newds.Destroy（）def polygon2mask（extent_vector，input_raster，extent_raster）：    “”“将具有多个功能的ESRI形状文件刻录到掩码文件        中输入栅格文件的扩展和分辨率。        多个特征由id的值表示    ：param str extent_vector：输入ESRI形状文件    ：param srt input_raster：基于烧录形状文件的输入栅格文件    ：param srt extent_raster：out文件名，格式：csv     ：param str字段：用于区分不同功能                 的字段，字段值应为数字类型，默认为“id”










    :return: status, 0 for success, 1 for failure
    """
    '''Open the shape file'''
    shp_driver = ogr.GetDriverByName("ESRI shapefile")
    shp_ds = shp_driver.Open(extent_vector)

    layer = shp_ds.GetLayer(0)

    '''Open the raster file'''
    raster_ds = gdal.Open(input_raster)
    ncol = raster_ds.RasterXSize
    nrow = raster_ds.RasterYSize
    proj = raster_ds.GetProjection()
    gt = raster_ds.GetGeoTransform()

    driver = gdal.GetDriverByName("GTiff")
    out_ds = driver.Create(extent_raster, ncol, nrow, 1, gdal.GDT_Byte)
    out_ds.SetProjection(proj)
    out_ds.SetGeoTransform(gt)
    band = out_ds.GetRasterBand(1)
    band.Fill(0)
    band.SetNoDataValue(0)

    #attribute = 'ATTRIBUTE=' + field
    status = gdal.RasterizeLayer(out_ds,  # output to our new dataset
                                 [1],  # output to our new dataset's first band
                                 layer,  # rasterize this layer
                                 None, None,  # No projection transform
                                 [1],  # burn value 0
                                 ['All_TOUCHED=TRUE'])  # rasterize all pixels touched by polygons
    del out_ds, shp_ds, raster_ds
    return status

def clipRaster(extent_raster,input_raster,clip_raster):

    extentArray = ga.LoadFile(extent_raster)
    inputArray = ga.LoadFile(input_raster)

    outputArray = extentArray*inputArray

    print(outputArray)

    out = ga.SaveArray(outputArray,clip_raster,format="GTiff",prototype=input_raster)    #prototype 原型设置为输入影像的格式
    out = None

if __name__ =='__main__':

    filterStr = "ID = 266"
    copyInLayer(input_vector, extent_vector, filterStr)
    polygon2mask(extent_vector, input_raster, extent_raster)
    clipRaster(extent_raster, input_raster, clip_raster)

    #copyInLayer(inshp, outputfile, filterStr)

    #polygon2mask(outputfile,rasterfile,outname)

    #clipRaster()

结果：
原始影像和矢量

裁切后的影像和矢量

待优化的问题：

 filterStr = "行政代码 = 211200"
 layer.SetAttributeFilter(filterStr)

gdal执行sql表达式，含有中文字段名的无法运行
C:\ProgramData\Anaconda3\python.exe F:/20180703/ClipRaster_f.py
ERROR 1: SQL Expression Parsing Error: syntax error, unexpected end of string. Occurred around :
行政代码 = 211200
————----------
C:\ProgramData\Anaconda3\python.exe F:/20180703/ClipRaster_f.py
ERROR 1: “行政代码” not recognised as an available field.

解决遇到的问题：
1.gdal中无法读取到中文属性值,使用SetConfigOption

gdal.SetConfigOption(‘SHAPE_ENCODING’, ‘’)

2.按照属性提取矢量,在layer层次上进行拷贝

layer.SetAttributeFilter(filterStr)
newds = driver.CreateDataSource（extent _ vector）> newds.CopyLayer（layer，’’）3。矢量存在交叉点无法裁切，使用矢量转栅格掩膜的方法> out _ ds = driver.Create（extent _ raster，ncol，nrow，1，gdal.GDT _ Byte）> status = gdal.RasterizeLayer（out _ ds，[1]，layer，None，None，[1]，[‘All _ TOUCHED = TRUE’]）4 。使用源数据信息输出裁切影像> ga.SaveArray（outputArray，clip _ raster，format =“GTiff”，prototype = input _ raster）