[Python]使用GDAL对遥感影像进行图像分块

在使用python进行遥感图像处理时，有时因为研究需要或内存限制需要将大研究区切分成一个个小研究区域，这时候就可以使用python代码批量分块。

本代码参考了http://t.csdn.cn/5c87N

十分感谢大佬们代码，我只是在原先的基础上进行了稍许完善

本代码主要完善了以下功能：

1.实现块与块不重叠的矩形分割

2.修改了原先分块后地理位置重叠的问题，该代码分块后各块地理位置投影完全正确

3.如果原影像长宽不是分块后小影像的长宽的整数倍，该代码也能完美解决

分块逻辑如下：

 

import os
import numpy as np
from osgeo import gdal


class GRID:

    def load_img(self, filename):
        image = gdal.Open(filename)
        img_width = image.RasterXSize
        img_height = image.RasterYSize
        img_geotrans = image.GetGeoTransform()
        img_proj = image.GetProjection()
        img_data = image.ReadAsArray(0, 0, img_width, img_height)
        del image
        return img_proj, img_geotrans, img_data

    def write_img(self, filename, img_proj, origin_x, origin_y, pixel_width, pixel_height,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((origin_x, pixel_width, 0, origin_y, 0, pixel_height))  # 写入仿射变换参数
        dataset.SetProjection(img_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


if __name__=='__main__':
    img_name='F:/gaobiao/7_1test/cropland3011_2.tif'
    save_path='F:/gaobiao/test/'
    proj,geotrans,data=GRID().load_img(img_name)
    dataset = gdal.Open(img_name)
    minx, xres, xskew, maxy, yskew, yres = dataset.GetGeoTransform()
    print(proj)
    print(geotrans)
    print(data.shape)
    channel,width,height=data.shape
    #x_size为沿x轴切分长度，y_size为沿y轴切分长度
    x_size = 900
    y_size = 600

    for j in range((height // x_size)+1):  # 按照像元数切割小图
        for i in range((width // y_size)+1):
            if j == height // x_size and i == width // y_size:
                cur_image = data[:, i * y_size:width, j * x_size: height]
                if cur_image.size == 0:
                    continue
                else:
                    lon = minx + xres * j * x_size
                    lat = maxy + yres * i * y_size
                    GRID().write_img(save_path+'test{}_{}.tif'.format(i, j), proj,
                                     lon, lat, xres, yres, cur_image)  ##写数据
            elif j == height // x_size:
                cur_image = data[:, i * y_size:width, j * x_size: (j + 1) * x_size]
                if cur_image.size == 0:
                    continue
                else:
                    lon =  minx + xres * x_size * j
                    lat = maxy + yres * i * y_size
                    GRID().write_img(save_path+'test{}_{}.tif'.format(i, j), proj,
                                     lon, lat, xres, yres, cur_image)  ##写数据
            elif i == width // y_size:
                cur_image = data[:, i * y_size:(i + 1) * y_size,j * x_size: height]
                if cur_image.size == 0:
                    continue
                else:
                    lon = minx + xres * j * x_size
                    lat = maxy + yres * (i * y_size)
                    GRID().write_img(save_path+'test{}_{}.tif'.format(i, j), proj,
                                     lon, lat, xres, yres, cur_image)  ##写数据
            #如果height、width恰好是x_size y_size的整数倍，就只用得着这一部分代码
            else:
                cur_image = data[:, i * y_size:(i + 1) * y_size, j * x_size: (j + 1) * x_size]
                #更改起始坐标位置
                lon = minx + xres * x_size * j
                lat = maxy + yres * (i * y_size)
                GRID().write_img(save_path+'test{}_{}.tif'.format(i, j), proj,
                                 lon, lat, xres, yres, cur_image)  ##写数据
    print('finish')

GRID类里面一个load_img读图像一个write_img写图像

write_img里面根据不同图像的起始坐标位置和长宽写仿射变换

最后的for循环对分割块进行判断，满足长宽要求的块直接走最后一部分语句，3_2那类的走第一个if，0_2那类的走第二个if，3_0那类的走第三个if。保证全部分割，不留一点黑边。