python代码对遥感图像进行裁剪并有坐标系

1、对地图的裁剪需要知道一个公式

• geoTransform[0]：左上角像素经度
• geoTransform[1]：影像宽度方向上的分辨率(经度范围/像素个数)
• geoTransform[2]：旋转, 0表示上面为北方
• geoTransform[3]：左上角像素纬度
• geoTransform[4]：旋转, 0表示上面为北方
• geoTransform[5]：影像宽度方向上的分辨率(纬度范围/像素个数)

之后对进行裁剪， 令px=geoTransform[0]   py=geoTransform[3]

#new_x是裁剪图像左上角距离原始图像的距离
#new_y是裁剪图像左上角距离原始图像
#上面的geoTransform就是这里的trans：
px = trans[0] + new_y* trans[1] + new_x * trans[2] 
py = trans[3] + new_y* trans[4] + new_x* trans[5]

之后设置新得坐标系就可里

new_geotransform = (px, trans[1], trans[2], py, trans[4], trans[5])

完整代码

from GetTiff import Tiff
from osgeo import gdal
import numpy as np
#这个就是对图像进行读取的类库
A=Tiff()

area=A.read_tif("C:\\Users\\Rare\\Desktop\\DSM.tif")
area1=area[2]
print(area1.shape)
row=0
new_x=int(area1.shape[0]/3.0*2)
new_y=int(area1.shape[1]/3*2)
row=0
mat=np.zeros((area1.shape[0]-new_x,area1.shape[1]-new_y),dtype=np.float32)
print(mat.shape)
for i in range(area1.shape[0]):
    line=0
    for j in range(area1.shape[1]):
        if(i>=new_x and j >=new_y):
            mat[row][line]=area1[i][j]
            line+=1
    if(line!=0):
        row+=1
print(row,line)
trans=area[1]
px = trans[0] + new_y* trans[1] + new_x * trans[2]
py = trans[3] + new_y* trans[4] + new_x* trans[5]
new_geotransform = (px, trans[1], trans[2], py, trans[4], trans[5])
A.write_tif("C:\\Users\\Rare\\Desktop\\area.tif",area[0],new_geotransform,mat)

 有些图像可能是黑色的这个跟拉伸比例有关，可以通过更改拉伸进行图像显示。