遙感影像按尺寸裁剪成若干個python的gdal


遙感影像 按照尺寸裁剪

因為我項目要求是1200*1200

但是這樣子還會有邊邊角角的地方我索性就還是把邊邊角角也一起弄進去了

如圖

 

所以! 上代碼!

首先! 感謝這個博主   https://blog.csdn.net/zsc201825/article/details/90721718

修改后的代碼!!!

 

# -*- coding: utf-8 -*-
import gdal
in_ds = gdal.Open(r'D:\graduate\GF1_WFV1_E120.0_N34.7_20180711_L1A0003315461.tiff')              # 讀取要切的原圖
print("open tif file succeed")
width = in_ds.RasterXSize                         # 獲取數據寬度   
height = in_ds.RasterYSize                        # 獲取數據高度
outbandsize = in_ds.RasterCount                   # 獲取數據波段數
im_geotrans = in_ds.GetGeoTransform()             # 獲取仿射矩陣信息
im_proj = in_ds.GetProjection()                   # 獲取投影信息
datatype = in_ds.GetRasterBand(1).DataType
im_data = in_ds.ReadAsArray()                     #獲取數據 

# 讀取原圖中的每個波段
in_band1 = in_ds.GetRasterBand(1)
in_band2 = in_ds.GetRasterBand(2)
in_band3 = in_ds.GetRasterBand(3)
in_band4 = in_ds.GetRasterBand(4)

# 定義切圖的起始點坐標
offset_x = 0  
offset_y = 0

# 定義切圖的大小(矩形框)
size = 1200
col_num = int(width / size)  #寬度可以分成幾塊
row_num = int(height / size) #高度可以分成幾塊
if(width % size != 0):
    col_num += 1
if(height % size != 0):
    row_num += 1
#這邊就知道我們一共是分成了多少個 如果說有多余的 那我們就讓那些也自己一小塊好吧
num = 1    #這個就用來記錄一共有多少塊的
# 現在我們知道的是寬度是1304  高度是666
print("row_num:%d   col_num:%d" %(row_num,col_num))
for i in range(row_num):    #從高度下手!!! 可以分成幾塊!
    for j in range(col_num):
        offset_x = i * size 
        offset_y = j * size
        ## 從每個波段中切需要的矩形框內的數據(注意讀取的矩形框不能超過原圖大小)
        b_ysize = min(width - offset_y, size)
        b_xsize = min(height - offset_x, size)

        print("width:%d     height:%d    offset_x:%d    offset_y:%d     b_xsize:%d     b_ysize:%d" %(width,height,offset_x,offset_y, b_xsize, b_ysize))
        # print("\n")
        out_band1 = in_band1.ReadAsArray(offset_y, offset_x, b_ysize, b_xsize)
        out_band2 = in_band2.ReadAsArray(offset_y, offset_x, b_ysize, b_xsize)
        out_band3 = in_band3.ReadAsArray(offset_y, offset_x, b_ysize, b_xsize)
        out_band4 = in_band4.ReadAsArray(offset_y, offset_x, b_ysize, b_xsize)
        # 獲取Tif的驅動,為創建切出來的圖文件做准備
        gtif_driver = gdal.GetDriverByName("GTiff")
        file = r'C:\Users\Administrator\Desktop\last\%04d.tiff' % num
        num += 1 
        # 創建切出來的要存的文件
        out_ds = gtif_driver.Create(file, b_ysize, b_xsize, outbandsize, datatype)
        print("create new tif file succeed")

        # 獲取原圖的原點坐標信息
        ori_transform = in_ds.GetGeoTransform()
        if ori_transform:
                print (ori_transform)
                print("Origin = ({}, {})".format(ori_transform[0], ori_transform[3]))
                print("Pixel Size = ({}, {})".format(ori_transform[1], ori_transform[5]))

        # 讀取原圖仿射變換參數值
        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,out_band1,out_band2,out_band3,out_band4

 


免責聲明!

本站轉載的文章為個人學習借鑒使用,本站對版權不負任何法律責任。如果侵犯了您的隱私權益,請聯系本站郵箱yoyou2525@163.com刪除。



 
粵ICP備18138465號   © 2018-2025 CODEPRJ.COM