shp裁剪tif

# -*- coding: utf-8 -*-
import os
import numpy as np
from osgeo import gdal, gdalnumeric, ogr, osr, gdal_array
gdal.UseExceptions()

def world2Pixel(geoMatrix, x, y):
  """
  Uses a gdal geomatrix (gdal.GetGeoTransform()) to calculate
  the pixel location of a geospatial coordinate
  """
  ulX = geoMatrix[0]
  ulY = geoMatrix[3]
  xDist = geoMatrix[1]
  yDist = geoMatrix[5]
  rtnX = geoMatrix[2]
  rtnY = geoMatrix[4]
  pixel = int((x - ulX) / xDist)
  line = int((ulY - y) / xDist)
  return (pixel, line)

#
#  EDIT: this is basically an overloaded
#  version of the gdal_array.OpenArray passing in xoff, yoff explicitly
#  so we can pass these params off to CopyDatasetInfo
#
def OpenArray( array, prototype_ds = None, xoff=0, yoff=0 ):
    # ds = gdal.Open( gdalnumeric.GetArrayFilename(array))
    ds = gdal_array.OpenArray(array)

    if ds is not None and prototype_ds is not None:
        if type(prototype_ds).__name__ == 'str':
            prototype_ds = gdal.Open( prototype_ds )
        if prototype_ds is not None:
            gdalnumeric.CopyDatasetInfo( prototype_ds, ds, xoff=xoff, yoff=yoff )
    return ds


def write_img(filename,im_proj,im_geotrans,im_data):
    if 'int8' in im_data.dtype.name:
        datatype = gdal.GDT_Byte
    elif 'int16' in im_data.dtype.name:
        datatype = gdal.GDT_UInt16
    else:
        datatype = gdal.GDT_Float32

    if len(im_data.shape) == 3:
        im_bands, im_height, im_width = im_data.shape
    else:
        im_bands, (im_height, im_width) = 1,im_data.shape 

    driver = gdal.GetDriverByName("GTiff")
    dataset = driver.Create(filename, im_width, im_height, im_bands, datatype)

    dataset.SetGeoTransform(im_geotrans)
    dataset.SetProjection(im_proj)
    if im_bands == 1:
        dataset.GetRasterBand(1).WriteArray(im_data)
    else:
        for i in range(im_bands):
            dataset.GetRasterBand(i+1).WriteArray(im_data[i])

    del dataset

def shpClipRaster(shapefile_path, raster_path, save_path):
    # Load the source data as a gdalnumeric array
    # srcArray = gdalnumeric.LoadFile(raster_path)

    # Also load as a gdal image to get geotransform
    # (world file) info
    srcImage = gdal.Open(raster_path)
    geoTrans = srcImage.GetGeoTransform()
    geoProj = srcImage.GetProjection()

    # Create an OGR layer from a boundary shapefile
    shapef = ogr.Open(shapefile_path)
    lyr = shapef.GetLayer( os.path.split( os.path.splitext( shapefile_path )[0] )[1] )
    poly = lyr.GetNextFeature()

    # Convert the layer extent to image pixel coordinates
    minX, maxX, minY, maxY = lyr.GetExtent()
    ulX, ulY = world2Pixel(geoTrans, minX, maxY)
    lrX, lrY = world2Pixel(geoTrans, maxX, minY)

    # Calculate the pixel size of the new image
    pxWidth = int(lrX - ulX)
    pxHeight = int(lrY - ulY)

    # clip = srcArray[:, ulY:lrY, ulX:lrX]
    clip = srcImage.ReadAsArray(ulX,ulY,pxWidth,pxHeight)   #***只读要的那块***

    #
    # EDIT: create pixel offset to pass to new image Projection info
    #
    xoffset =  ulX
    yoffset =  ulY
    print ("Xoffset, Yoffset = ( %f, %f )" % ( xoffset, yoffset ))

    # Create a new geomatrix for the image
    geoTrans = list(geoTrans)
    geoTrans[0] = minX
    geoTrans[3] = maxY

    write_img(save_path, geoProj, geoTrans, clip)
    gdal.ErrorReset()

if __name__ == "__main__":
    shp = "dataset/E22_Bound.shp"
    img = "dataset/CGdomYRJ-114(CK0-17)_E_22.tif"
    out = "dataset/E22.tif"

    shpClipRaster(shp,img,out)
    print(img)

shp裁剪shp

import os
from osgeo import gdal, ogr

def ShapeClip(
		baseFilePath,
		maskFilePath,
		saveFolderPath):
	"""
	矢量裁剪
	:param baseFilePath: 要裁剪的矢量文件
	:param maskFilePath: 掩膜矢量文件
	:param saveFolderPath: 裁剪后的矢量文件保存目录
	:return:
	"""
	ogr.RegisterAll()
	gdal.SetConfigOption("GDAL_FILENAME_IS_UTF8", "YES")
	# 载入要裁剪的矢量文件

	baseData = ogr.Open(baseFilePath)
	print(os.path.split( os.path.splitext( baseFilePath )[0] )[1])
	baseLayer = baseData.GetLayer( os.path.split( os.path.splitext( baseFilePath )[0] )[1] )

	spatial = baseLayer.GetSpatialRef()
	geomType = baseLayer.GetGeomType()
	baseLayerName = baseLayer.GetName()
	# 载入掩膜矢量文件
	maskData = ogr.Open(maskFilePath)
	maskLayer = maskData.GetLayer()
	maskLayerName = maskLayer.GetName()
	# 生成裁剪后的矢量文件
	outLayerName = maskLayerName + "_Clip_" + baseLayerName
	outFilePath = saveFolderPath
	gdal.SetConfigOption("SHAPE_ENCODING", "GBK")
	driver = ogr.GetDriverByName("ESRI Shapefile")
	outData = driver.CreateDataSource(outFilePath)
	outLayer = outData.CreateLayer(outLayerName, spatial, geomType)
	baseLayer.Clip(maskLayer, outLayer)
	outData.Release()
	baseData.Release()
	maskData.Release()
	return outFilePath


if __name__ == "__main__":
	baseFilePath = 'dataset/veg_E_22.shp'
	maskFilePath = 'dataset/E22_Bound.shp'
	saveFolderPath = 'dataset/E22.shp'
	outFilePath=ShapeClip(baseFilePath,maskFilePath,saveFolderPath)
	print(outFilePath)

矢量转栅格

from osgeo import gdal, ogr, gdalconst
def shp2Raster(shp,templatePic,output,nodata):
    """
    shp:字符串，一个矢量，从0开始计数，整数
    templatePic:字符串，模板栅格，一个tif，地理变换信息从这里读，栅格大小与该栅格一致
    output:字符串，输出栅格，一个tif
    field:字符串，栅格值的字段
    nodata:整型或浮点型，矢量空白区转换后的值
    """
    ndsm = templatePic
    data = gdal.Open(ndsm, gdalconst.GA_ReadOnly)
    geo_transform = data.GetGeoTransform()
    proj=data.GetProjection()
    #source_layer = data.GetLayer()
    x_min = geo_transform[0]
    y_max = geo_transform[3]
    x_max = x_min + geo_transform[1] * data.RasterXSize
    y_min = y_max + geo_transform[5] * data.RasterYSize
    x_res = data.RasterXSize
    y_res = data.RasterYSize
    mb_v = ogr.Open(shp)
    mb_l = mb_v.GetLayer()
    pixel_width = geo_transform[1]
    #输出影像为24位整型
    target_ds = gdal.GetDriverByName('GTiff').Create(output, x_res, y_res, 1, gdal.GPI_RGB)

    target_ds.SetGeoTransform(geo_transform)
    target_ds.SetProjection(proj)
    band = target_ds.GetRasterBand(1)
    NoData_value = nodata
    band.SetNoDataValue(NoData_value)
    band.FlushCache()
    gdal.RasterizeLayer(target_ds, [1], mb_l, options=['ALL_TOUCHED=TRUE'])

    target_ds = None

if __name__ == "__main__":
    shp = "dataset/E22.shp"
    templatePic= "dataset/E22.tif"
    output = "dataset/E22_mask.tif"
    nodata=0
    shp2Raster(shp,templatePic,output,nodata)
    

多个shp裁剪shp

import os
import os
import numpy as np
from osgeo import gdal, gdalnumeric, ogr, osr, gdal_array
gdal.UseExceptions()

def world2Pixel(geoMatrix, x, y):
  """
  Uses a gdal geomatrix (gdal.GetGeoTransform()) to calculate
  the pixel location of a geospatial coordinate
  """
  ulX = geoMatrix[0]
  ulY = geoMatrix[3]
  xDist = geoMatrix[1]
  yDist = geoMatrix[5]
  rtnX = geoMatrix[2]
  rtnY = geoMatrix[4]
  pixel = int((x - ulX) / xDist)
  line = int((ulY - y) / xDist)
  return (pixel, line)

#
#  EDIT: this is basically an overloaded
#  version of the gdal_array.OpenArray passing in xoff, yoff explicitly
#  so we can pass these params off to CopyDatasetInfo
#
def OpenArray( array, prototype_ds = None, xoff=0, yoff=0 ):
    # ds = gdal.Open( gdalnumeric.GetArrayFilename(array))
    ds = gdal_array.OpenArray(array)

    if ds is not None and prototype_ds is not None:
        if type(prototype_ds).__name__ == 'str':
            prototype_ds = gdal.Open( prototype_ds )
        if prototype_ds is not None:
            gdalnumeric.CopyDatasetInfo( prototype_ds, ds, xoff=xoff, yoff=yoff )
    return ds


def write_img(filename,im_proj,im_geotrans,im_data):
    if 'int8' in im_data.dtype.name:
        datatype = gdal.GDT_Byte
    elif 'int16' in im_data.dtype.name:
        datatype = gdal.GDT_UInt16
    else:
        datatype = gdal.GDT_Float32

    if len(im_data.shape) == 3:
        im_bands, im_height, im_width = im_data.shape
    else:
        im_bands, (im_height, im_width) = 1,im_data.shape 

    driver = gdal.GetDriverByName("GTiff")
    dataset = driver.Create(filename, im_width, im_height, im_bands, datatype)

    dataset.SetGeoTransform(im_geotrans)
    dataset.SetProjection(im_proj)
    if im_bands == 1:
        dataset.GetRasterBand(1).WriteArray(im_data)
    else:
        for i in range(im_bands):
            dataset.GetRasterBand(i+1).WriteArray(im_data[i])

    del dataset

pre_path='dataset/pre/'
labellist = filter(lambda x: x.find('label')!=-1, os.listdir(pre_path))
list1 = list(map(lambda x: x[:], labellist))
label_name=pre_path +  list1[0]

boundarylist = filter(lambda x: x.find('shp')!=-1, os.listdir(pre_path+'boundary/'))
list2 = list(map(lambda x: x[:], boundarylist))


imagelist = filter(lambda x: x.find('tif')!=-1, os.listdir(pre_path))
list3 = list(map(lambda x: x[:], imagelist))
img_path=pre_path +  list3[0]


"""
矢量裁剪
:param label_name: 要裁剪的矢量文件
:param boundary_name: 掩膜矢量文件
img_path: 影像
:param saveFolderPath: 裁剪后的矢量文件保存目录
:return:
"""
ogr.RegisterAll()
gdal.SetConfigOption("GDAL_FILENAME_IS_UTF8", "YES")
# 载入要裁剪的矢量文件

labelData = ogr.Open(label_name)

labelLayer = labelData.GetLayer( os.path.split( os.path.splitext( label_name )[0] )[1] )

spatial = labelLayer.GetSpatialRef()
geomType = labelLayer.GetGeomType()


# 载入掩膜矢量文件

def new_func(outLayerName):
    return outLayerName

for i in list2:
    boundary_name=pre_path+'boundary/'+ i
    maskData = ogr.Open(boundary_name)
    maskLayer = maskData.GetLayer()
    #裁剪shp
    # 生成裁剪后的矢量文件
    save_shp_dir='./dataset/pre/shp/'
    if not os.path.exists(save_shp_dir):
        os.mkdir(save_shp_dir)
    outLayerName = (save_shp_dir+i)
    gdal.SetConfigOption("SHAPE_ENCODING", "GBK")
    driver = ogr.GetDriverByName("ESRI Shapefile")
    outData = driver.CreateDataSource(outLayerName)
    outLayer = outData.CreateLayer(new_func(outLayerName), spatial, geomType)
    labelLayer.Clip(maskLayer, outLayer)
    outData.Release()
    maskData.Release()

    #裁剪tif

    shp = "dataset/E22_Bound.shp"
    img = "dataset/CGdomYRJ-114(CK0-17)_E_22.tif"
    out = "dataset/E22.tif"
    # Load the source data as a gdalnumeric array
    # srcArray = gdalnumeric.LoadFile(raster_path)

    # Also load as a gdal image to get geotransform
    # (world file) info
    srcImage = gdal.Open(raster_path)
    geoTrans = srcImage.GetGeoTransform()
    geoProj = srcImage.GetProjection()

    # Create an OGR layer from a boundary shapefile
    shapef = ogr.Open(shapefile_path)
    lyr = shapef.GetLayer( os.path.split( os.path.splitext( shapefile_path )[0] )[1] )
    poly = lyr.GetNextFeature()

    # Convert the layer extent to image pixel coordinates
    minX, maxX, minY, maxY = lyr.GetExtent()
    ulX, ulY = world2Pixel(geoTrans, minX, maxY)
    lrX, lrY = world2Pixel(geoTrans, maxX, minY)

    # Calculate the pixel size of the new image
    pxWidth = int(lrX - ulX)
    pxHeight = int(lrY - ulY)

    # clip = srcArray[:, ulY:lrY, ulX:lrX]
    clip = srcImage.ReadAsArray(ulX,ulY,pxWidth,pxHeight)   #***只读要的那块***

    #
    # EDIT: create pixel offset to pass to new image Projection info
    #
    xoffset =  ulX
    yoffset =  ulY
    print ("Xoffset, Yoffset = ( %f, %f )" % ( xoffset, yoffset ))

    # Create a new geomatrix for the image
    geoTrans = list(geoTrans)
    geoTrans[0] = minX
    geoTrans[3] = maxY

    write_img(save_path, geoProj, geoTrans, clip)
    gdal.ErrorReset()
labelData.Release()

汇总：输入shp和影像，输出影像对应的标签

#影像裁剪shp,转为栅格，为该影像标签
#输入：存放影像文件夹dataset/sat_train，存放标签矢量文件夹dataset/mask_shp
#输出：标签（栅格），存放在dataset/mask_train

from osgeo import gdal, ogr, osr, gdalconst
import fnmatch
import os

def ShapeClip(
		baseFilePath,
		maskFilePath,
		saveFolderPath):
	"""
	矢量裁剪
	:param baseFilePath: 要裁剪的矢量文件
	:param maskFilePath: 掩膜矢量文件
	:param saveFolderPath: 裁剪后的矢量文件保存目录
	:return:
	"""
	ogr.RegisterAll()
	gdal.SetConfigOption("GDAL_FILENAME_IS_UTF8", "YES")
	# 载入要裁剪的矢量文件

	baseData = ogr.Open(baseFilePath)
	baseLayer = baseData.GetLayer( os.path.split( os.path.splitext( baseFilePath )[0] )[1] )

	spatial = baseLayer.GetSpatialRef()
	geomType = baseLayer.GetGeomType()
	baseLayerName = baseLayer.GetName()
	# 载入掩膜矢量文件
	maskData = ogr.Open(maskFilePath)
	maskLayer = maskData.GetLayer()
	maskLayerName = maskLayer.GetName()
	# 生成裁剪后的矢量文件
	outLayerName = maskLayerName + "_Clip_" + baseLayerName
	outFilePath = saveFolderPath
	gdal.SetConfigOption("SHAPE_ENCODING", "GBK")
	driver = ogr.GetDriverByName("ESRI Shapefile")
	outData = driver.CreateDataSource(outFilePath)
	outLayer = outData.CreateLayer(outLayerName, spatial, geomType)
	baseLayer.Clip(maskLayer, outLayer)
	outData.Release()
	baseData.Release()
	maskData.Release()
	return outFilePath

def shp2Raster(shp,templatePic,output,nodata):
    """
    shp:字符串，一个矢量，从0开始计数，整数
    templatePic:字符串，模板栅格，一个tif，地理变换信息从这里读，栅格大小与该栅格一致
    output:字符串，输出栅格，一个tif
    field:字符串，栅格值的字段
    nodata:整型或浮点型，矢量空白区转换后的值
    """
    ndsm = templatePic
    data = gdal.Open(ndsm, gdalconst.GA_ReadOnly)
    geo_transform = data.GetGeoTransform()
    proj=data.GetProjection()
    #source_layer = data.GetLayer()
    x_min = geo_transform[0]
    y_max = geo_transform[3]
    x_max = x_min + geo_transform[1] * data.RasterXSize
    y_min = y_max + geo_transform[5] * data.RasterYSize
    x_res = data.RasterXSize
    y_res = data.RasterYSize
    mb_v = ogr.Open(shp)
    mb_l = mb_v.GetLayer()
    pixel_width = geo_transform[1]
    #输出影像为24位整型
    target_ds = gdal.GetDriverByName('GTiff').Create(output, x_res, y_res, 1, gdal.GPI_RGB)

    target_ds.SetGeoTransform(geo_transform)
    target_ds.SetProjection(proj)
    band = target_ds.GetRasterBand(1)
    NoData_value = nodata
    band.SetNoDataValue(NoData_value)
    band.FlushCache()
    gdal.RasterizeLayer(target_ds, [1], mb_l, options=['ALL_TOUCHED=TRUE'])

    target_ds = None

print("开始制作标签")
ogr.RegisterAll()
img_path="dataset/sat_train/" #影像所在的文件夹
mask_shp_path="dataset/mask_shp/" #原始标签shp位置

shape_path="dataset/mask_boundary_shp/" #shape输出位置
mask_clip_path='dataset/mask_clip_train/'#裁剪后shp
mask_train_path='dataset/mask_train/'#最终输出标签文件夹
if not os.path.exists(shape_path):
    os.mkdir(shape_path)
if not os.path.exists(mask_clip_path):
    os.mkdir(mask_clip_path)

imagelist = filter(lambda x: x.find('shp')!=-1, os.listdir(mask_shp_path))
list = list(map(lambda x: x[:], imagelist))
mask_shp_name=mask_shp_path +  list[0]
img_list = fnmatch.filter(os.listdir(img_path), '*.tif')
for img in img_list:
    p_img=img_path+img
    outfilename = shape_path+img[:-4]+".shp"
    dataset = gdal.Open(p_img)
    oDriver = ogr.GetDriverByName('ESRI Shapefile')
    oDS = oDriver.CreateDataSource(outfilename)
    srs = osr.SpatialReference(wkt=dataset.GetProjection())
    geocd = dataset.GetGeoTransform()
    oLayer = oDS.CreateLayer("polygon", srs, ogr.wkbPolygon)
    oDefn = oLayer.GetLayerDefn()
    row = dataset.RasterXSize
    line = dataset.RasterYSize
    geoxmin = geocd[0]
    geoymin = geocd[3]
    geoxmax = geocd[0] + (row) * geocd[1] + (line) * geocd[2]
    geoymax = geocd[3] + (row) * geocd[4] + (line) * geocd[5]
    ring = ogr.Geometry(ogr.wkbLinearRing)
    ring.AddPoint(geoxmin, geoymin)
    ring.AddPoint(geoxmax, geoymin)
    ring.AddPoint(geoxmax, geoymax)
    ring.AddPoint(geoxmin, geoymax)
    ring.CloseRings()
    poly = ogr.Geometry(ogr.wkbPolygon)
    poly.AddGeometry(ring)
    outfeat = ogr.Feature(oDefn)
    outfeat.SetGeometry(poly)
    oLayer.CreateFeature(outfeat)
    outfeat = None
    oDS.Destroy()
    mask_train_name = mask_clip_path+img[:-4]+".shp"
    #裁剪
    outFilePath=ShapeClip(mask_shp_name,outfilename, mask_train_name)
    #矢量转栅格
    output = mask_train_path + img
    nodata=0
    shp2Raster(mask_train_name,p_img,output,nodata)
    print(output)
    
print('标签制作完成')

 做mask

'根据多个给定范围shp，对画好的标签进行裁剪并转栅格，做为标签样本，对影像进行裁剪，作为影像样本'
'输入：'
'输出'
import os
import os
import numpy as np
from osgeo import gdal, gdalnumeric, ogr, osr, gdal_array
gdal.UseExceptions()

def world2Pixel(geoMatrix, x, y):
  """
  Uses a gdal geomatrix (gdal.GetGeoTransform()) to calculate
  the pixel location of a geospatial coordinate
  """
  ulX = geoMatrix[0]
  ulY = geoMatrix[3]
  xDist = geoMatrix[1]
  yDist = geoMatrix[5]
  rtnX = geoMatrix[2]
  rtnY = geoMatrix[4]
  pixel = int((x - ulX) / xDist)
  line = int((ulY - y) / xDist)
  return (pixel, line)

#
#  EDIT: this is basically an overloaded
#  version of the gdal_array.OpenArray passing in xoff, yoff explicitly
#  so we can pass these params off to CopyDatasetInfo
#
def OpenArray( array, prototype_ds = None, xoff=0, yoff=0 ):
    # ds = gdal.Open( gdalnumeric.GetArrayFilename(array))
    ds = gdal_array.OpenArray(array)

    if ds is not None and prototype_ds is not None:
        if type(prototype_ds).__name__ == 'str':
            prototype_ds = gdal.Open( prototype_ds )
        if prototype_ds is not None:
            gdalnumeric.CopyDatasetInfo( prototype_ds, ds, xoff=xoff, yoff=yoff )
    return ds


def write_img(filename,im_proj,im_geotrans,im_data):
    if 'int8' in im_data.dtype.name:
        datatype = gdal.GDT_Byte
    elif 'int16' in im_data.dtype.name:
        datatype = gdal.GDT_UInt16
    else:
        datatype = gdal.GDT_Float32

    if len(im_data.shape) == 3:
        im_bands, im_height, im_width = im_data.shape
    else:
        im_bands, (im_height, im_width) = 1,im_data.shape 

    driver = gdal.GetDriverByName("GTiff")
    dataset = driver.Create(filename, im_width, im_height, im_bands, datatype)

    dataset.SetGeoTransform(im_geotrans)
    dataset.SetProjection(im_proj)
    if im_bands == 1:
        dataset.GetRasterBand(1).WriteArray(im_data)
    else:
        for i in range(im_bands):
            dataset.GetRasterBand(i+1).WriteArray(im_data[i])

    del dataset


def shp2Raster(shp,templatePic,output,nodata):
    """
    shp:字符串，一个矢量，从0开始计数，整数
    templatePic:字符串，模板栅格，一个tif，地理变换信息从这里读，栅格大小与该栅格一致
    output:字符串，输出栅格，一个tif
    field:字符串，栅格值的字段
    nodata:整型或浮点型，矢量空白区转换后的值
    """
    ndsm = templatePic
    data = gdal.Open(ndsm, gdalconst.GA_ReadOnly)
    geo_transform = data.GetGeoTransform()
    proj=data.GetProjection()
    #source_layer = data.GetLayer()
    x_min = geo_transform[0]
    y_max = geo_transform[3]
    x_max = x_min + geo_transform[1] * data.RasterXSize
    y_min = y_max + geo_transform[5] * data.RasterYSize
    x_res = data.RasterXSize
    y_res = data.RasterYSize
    mb_v = ogr.Open(shp)
    mb_l = mb_v.GetLayer()
    pixel_width = geo_transform[1]
    #输出影像为24位整型
    target_ds = gdal.GetDriverByName('GTiff').Create(output, x_res, y_res, 1, gdal.GPI_RGB)
 
    target_ds.SetGeoTransform(geo_transform)
    target_ds.SetProjection(proj)
    band = target_ds.GetRasterBand(1)
    NoData_value = nodata
    band.SetNoDataValue(NoData_value)
    band.FlushCache()
    gdal.RasterizeLayer(target_ds, [1], mb_l, options=['ALL_TOUCHED=TRUE'])
 
    target_ds = None


pre_path='dataset/pre/'

mask_train_path='dataset/mask_train/'#最终输出标签文件夹

labellist = filter(lambda x: x.find('label')!=-1, os.listdir(pre_path))
list1 = list(map(lambda x: x[:], labellist))
label_name=pre_path +  list1[0]

boundarylist = filter(lambda x: x.find('.shp')!=-1, os.listdir(pre_path+'boundary/'))
list2 = list(map(lambda x: x[:], boundarylist))


imagelist = filter(lambda x: x.find('tif')!=-1, os.listdir(pre_path+'img'))
list3 = list(map(lambda x: x[:], imagelist))
img_path=pre_path +  list3[0]


"""
矢量裁剪
:param label_name: 要裁剪的矢量文件
:param boundary_name: 掩膜矢量文件
img_path: 影像
:param saveFolderPath: 裁剪后的矢量文件保存目录
:return:
"""
print('开始用矢量范围裁剪影像')
ogr.RegisterAll()
gdal.SetConfigOption("GDAL_FILENAME_IS_UTF8", "YES")
# 载入要裁剪的矢量文件

labelData = ogr.Open(label_name)

labelLayer = labelData.GetLayer( os.path.split( os.path.splitext( label_name )[0] )[1] )

spatial = labelLayer.GetSpatialRef()
geomType = labelLayer.GetGeomType()


# 载入掩膜矢量文件

def new_func(outLayerName):
    return outLayerName

for i in list2:
    boundary_name=pre_path+'boundary/'+ i
    maskData = ogr.Open(boundary_name)
    maskLayer = maskData.GetLayer()
    #裁剪shp
    # 生成裁剪后的矢量文件
    save_shp_dir='./dataset/pre/shp/'
    if not os.path.exists(save_shp_dir):
        os.mkdir(save_shp_dir)
    outLayerName = (save_shp_dir+i)
    gdal.SetConfigOption("SHAPE_ENCODING", "GBK")
    driver = ogr.GetDriverByName("ESRI Shapefile")
    outData = driver.CreateDataSource(outLayerName)
    outLayer = outData.CreateLayer(new_func(outLayerName), spatial, geomType)
    labelLayer.Clip(maskLayer, outLayer)

    lyr = maskData.GetLayer( os.path.split( os.path.splitext( boundary_name )[0] )[1] )
    shpminX, shpmaxX, shpminY, shpmaxY = lyr.GetExtent()



    #裁剪tif
    flag=0
    for j in list3:
        raster_path = pre_path+'img/'+j
        srcImage = gdal.Open(raster_path)
        geocd = srcImage.GetGeoTransform()
        geoProj = srcImage.GetProjection()
        row = srcImage.RasterXSize
        line = srcImage.RasterYSize
        tifxmin = geocd[0]
        tifymin = geocd[3]
        tifxmax = geocd[0] + (row) * geocd[1] + (line) * geocd[2]
        tifymax = geocd[3] + (row) * geocd[4] + (line) * geocd[5]
        if shpminX>=tifxmin and shpmaxX<=tifxmax and shpminY<=tifymin and shpmaxY>=tifymax:
            ulX, ulY = world2Pixel(geocd, shpminX, shpmaxY)
            lrX, lrY = world2Pixel(geocd, shpmaxX, shpminY)
            # Calculate the pixel size of the new image
            pxWidth = int(lrX - ulX)
            pxHeight = int(lrY - ulY)
            clip = srcImage.ReadAsArray(ulX,ulY,pxWidth,pxHeight)   #***只读要的那块***
            xoffset =  ulX
            yoffset =  ulY
            geoTrans = list(geoTrans)
            geoTrans[0] = shpminX
            geoTrans[3] = shpmaxY
            save_path='dataset/sat_train/'+i[:-4]+'.tif'
            write_img(save_path, geoProj, geoTrans, clip)
            gdal.ErrorReset()
            outData.Release()
            maskData.Release()
            flag=1
            output = mask_train_path + i[:-4] +'.tif'
            nodata=0
            shp2Raster(outLayerName,save_path,output,nodata)
    if flag==0:
        print(raster_path+"没有制作")
    else:
        print(raster_path)
labelData.Release()
    

 做了一半的

import os
import os
import numpy as np
from osgeo import gdal, gdalnumeric, ogr, osr, gdal_array
gdal.UseExceptions()

def world2Pixel(geoMatrix, x, y):
  """
  Uses a gdal geomatrix (gdal.GetGeoTransform()) to calculate
  the pixel location of a geospatial coordinate
  """
  ulX = geoMatrix[0]
  ulY = geoMatrix[3]
  xDist = geoMatrix[1]
  yDist = geoMatrix[5]
  rtnX = geoMatrix[2]
  rtnY = geoMatrix[4]
  pixel = int((x - ulX) / xDist)
  line = int((ulY - y) / xDist)
  return (pixel, line)

#
#  EDIT: this is basically an overloaded
#  version of the gdal_array.OpenArray passing in xoff, yoff explicitly
#  so we can pass these params off to CopyDatasetInfo
#
def OpenArray( array, prototype_ds = None, xoff=0, yoff=0 ):
    # ds = gdal.Open( gdalnumeric.GetArrayFilename(array))
    ds = gdal_array.OpenArray(array)

    if ds is not None and prototype_ds is not None:
        if type(prototype_ds).__name__ == 'str':
            prototype_ds = gdal.Open( prototype_ds )
        if prototype_ds is not None:
            gdalnumeric.CopyDatasetInfo( prototype_ds, ds, xoff=xoff, yoff=yoff )
    return ds


def write_img(filename,im_proj,im_geotrans,im_data):
    if 'int8' in im_data.dtype.name:
        datatype = gdal.GDT_Byte
    elif 'int16' in im_data.dtype.name:
        datatype = gdal.GDT_UInt16
    else:
        datatype = gdal.GDT_Float32

    if len(im_data.shape) == 3:
        im_bands, im_height, im_width = im_data.shape
    else:
        im_bands, (im_height, im_width) = 1,im_data.shape 

    driver = gdal.GetDriverByName("GTiff")
    dataset = driver.Create(filename, im_width, im_height, im_bands, datatype)

    dataset.SetGeoTransform(im_geotrans)
    dataset.SetProjection(im_proj)
    if im_bands == 1:
        dataset.GetRasterBand(1).WriteArray(im_data)
    else:
        for i in range(im_bands):
            dataset.GetRasterBand(i+1).WriteArray(im_data[i])

    del dataset

pre_path='dataset/pre/'
labellist = filter(lambda x: x.find('label')!=-1, os.listdir(pre_path))
list1 = list(map(lambda x: x[:], labellist))
label_name=pre_path +  list1[0]

boundarylist = filter(lambda x: x.find('.shp')!=-1, os.listdir(pre_path+'boundary/'))
list2 = list(map(lambda x: x[:], boundarylist))


imagelist = filter(lambda x: x.find('tif')!=-1, os.listdir(pre_path+'img'))
list3 = list(map(lambda x: x[:], imagelist))
img_path=pre_path +  list3[0]


"""
矢量裁剪
:param label_name: 要裁剪的矢量文件
:param boundary_name: 掩膜矢量文件
img_path: 影像
:param saveFolderPath: 裁剪后的矢量文件保存目录
:return:
"""
print('开始用矢量范围裁剪影像')
ogr.RegisterAll()
gdal.SetConfigOption("GDAL_FILENAME_IS_UTF8", "YES")
# 载入要裁剪的矢量文件

labelData = ogr.Open(label_name)

labelLayer = labelData.GetLayer( os.path.split( os.path.splitext( label_name )[0] )[1] )

spatial = labelLayer.GetSpatialRef()
geomType = labelLayer.GetGeomType()


# 载入掩膜矢量文件

def new_func(outLayerName):
    return outLayerName

for i in list2:
    boundary_name=pre_path+'boundary/'+ i
    maskData = ogr.Open(boundary_name)
    maskLayer = maskData.GetLayer()
    #裁剪shp
    # 生成裁剪后的矢量文件
    save_shp_dir='./dataset/pre/shp/'
    if not os.path.exists(save_shp_dir):
        os.mkdir(save_shp_dir)
    outLayerName = (save_shp_dir+i)
    gdal.SetConfigOption("SHAPE_ENCODING", "GBK")
    driver = ogr.GetDriverByName("ESRI Shapefile")
    outData = driver.CreateDataSource(outLayerName)
    outLayer = outData.CreateLayer(new_func(outLayerName), spatial, geomType)
    labelLayer.Clip(maskLayer, outLayer)
    
    lyr = maskData.GetLayer( os.path.split( os.path.splitext( boundary_name )[0] )[1] )
    shpminX, shpmaxX, shpminY, shpmaxY = lyr.GetExtent()



    #裁剪tif
    flag=0
    for j in list3:
        raster_path = pre_path+'img/'+j
        srcImage = gdal.Open(raster_path)
        geocd = srcImage.GetGeoTransform()
        geoProj = srcImage.GetProjection()
        row = srcImage.RasterXSize
        line = srcImage.RasterYSize
        tifxmin = geocd[0]
        tifymin = geocd[3]
        tifxmax = geocd[0] + (row) * geocd[1] + (line) * geocd[2]
        tifymax = geocd[3] + (row) * geocd[4] + (line) * geocd[5]
        if shpminX>=tifxmin and shpmaxX<=tifxmax and shpminY<=tifymin and shpmaxY>=tifymax:
            ulX, ulY = world2Pixel(geocd, shpminX, shpmaxY)
            lrX, lrY = world2Pixel(geocd, shpmaxX, shpminY)
            # Calculate the pixel size of the new image
            pxWidth = int(lrX - ulX)
            pxHeight = int(lrY - ulY)
            clip = srcImage.ReadAsArray(ulX,ulY,pxWidth,pxHeight)   #***只读要的那块***
            xoffset =  ulX
            yoffset =  ulY
            geocd = list(geocd)
            geocd[0] = shpminX
            geocd[3] = shpmaxY
            save_path='dataset/sat_train/'+i[:-4]+'.tif'
            write_img(save_path, geoProj, geocd, clip)
            gdal.ErrorReset()
            outData.Release()
            maskData.Release()
            flag=1
    if flag==0:
        print(raster_path+"没有制作")
    else:
        print(raster_path)
labelData.Release()