python gdal 基于栅格shp文件裁剪geotif图

import numpy as np
import gdal
import operator
from functools import reduce
import os
import shapefile

def geo2imagexy(dataset, x, y):
‘’’
根据GDAL的六 参数模型将给定的投影或地理坐标转为影像图上坐标
:param dataset: GDAL地理数据
:param x: 投影或地理坐标x
:param y: 投影或地理坐标y
:return: 影坐标或地理坐标(x, y)对应的影像图上坐标
‘’’
trans = dataset.GetGeoTransform()
a = np.array([[trans[1], trans[2]], [trans[4], trans[5]]])
b = np.array([x - trans[0], y - trans[3]])
return np.linalg.solve(a, b) # 使用numpy的linalg.solve进行二元一次方程的求解

#读取要裁剪的原图
in_ds = gdal.Open(“raw.tif”)
print(“open tif file succeed”)

#读取原图中的每个波段
in_band1 = in_ds.GetRasterBand(1)
in_band2 = in_ds.GetRasterBand(2)
in_band3 = in_ds.GetRasterBand(3)
in_band4 = in_ds.GetRasterBand(4)

#导入栅格shp文件
shp_path=‘D:/’
r = shapefile.Reader(shp_path+‘1.shp’)
minX, minY, maxX, maxY = r.bbox

#定义切图的起始点和终点坐标(相比原点的横坐标和纵坐标偏移量)
offset_x,offset_y=geo2imagexy(in_ds,minX,minY)
endset_x,endset_y=geo2imagexy(in_ds,maxX,maxY)

#定义切图的大小（矩形框）
block_xsize=int(endset_x-offset_x+1)
block_ysize=int(abs(endset_y-offset_y)+1)

#从每个波段中裁剪需要的矩形框内的数据
out_band1 = in_band1.ReadAsArray(offset_x, offset_y, block_xsize, block_ysize)
out_band2 = in_band2.ReadAsArray(offset_x, offset_y, block_xsize, block_ysize)
out_band3 = in_band3.ReadAsArray(offset_x, offset_y, block_xsize, block_ysize)
out_band4 = in_band4.ReadAsArray(offset_x, offset_y, block_xsize, block_ysize)

#获取Tif的驱动，为创建切出来的图文件做准备
gtif_driver = gdal.GetDriverByName(“GTiff”)

#创建切出来的要存的文件（4代表4个不都按，最后一个参数为数据类型，跟原文件一致）
out_ds = gtif_driver.Create(‘clip.tif’, block_xsize, block_ysize, 4, in_band1.DataType)
print(“create new tif file succeed”)

#获取原图的原点坐标信息
ori_transform = in_ds.GetGeoTransform()

#读取原图仿射变换参数值
top_left_x = ori_transform[0] # 左上角x坐标
w_e_pixel_resolution = ori_transform[1] # 东西方向像素分辨率
top_left_y = ori_transform[3] # 左上角y坐标
n_s_pixel_resolution = ori_transform[5] # 南北方向像素分辨率

#根据反射变换参数计算新图的原点坐标
top_left_x = top_left_x + offset_x * w_e_pixel_resolution
top_left_y = top_left_y + offset_y * n_s_pixel_resolution

#将计算后的值组装为一个元组，以方便设置
dst_transform = (top_left_x, ori_transform[1], ori_transform[2], top_left_y, ori_transform[4], ori_transform[5])

#设置裁剪出来图的原点坐标
out_ds.SetGeoTransform(dst_transform)

#设置SRS属性（投影信息）
out_ds.SetProjection(in_ds.GetProjection())

#写入目标文件
out_ds.GetRasterBand(1).WriteArray(out_band1)
out_ds.GetRasterBand(2).WriteArray(out_band2)
out_ds.GetRasterBand(3).WriteArray(out_band3)
out_ds.GetRasterBand(4).WriteArray(out_band4)

#将缓存写入磁盘
out_ds.FlushCache()
print(“FlushCache succeed”)
del out_ds
print(“End!”)