Python 将tif文件分割成多个小tif

个人整理的python地理信息库,本文代码也在其中，欢迎使用: https://github.com/Adam0429/geo-py

可能安装依赖或者使用以下程序时大家会遇到问题，本人提供有偿远程调试，需要的可以私信我

顺便帮公司打个广告，欢迎遥感类专家的加入～～

投简历戳这里，英视睿达欢迎您！

原始图像是这样的
目标是剪裁出指定大小的若干小tif

代码很简单，直接放了

from tqdm import tqdm
from osgeo import gdal,ogr,osr
import numpy as np
from glob import glob

def read_tif(tif_path):
    ds = gdal.Open(tiff_path)
    row = ds.RasterXSize
    col = ds.RasterYSize
    band = ds.RasterCount

    for i in range(band):
        data = ds.GetRasterBand(i+1).ReadAsArray()

        data = np.expand_dims(data , 2)
        if i == 0:
            allarrays = data
        else:
            allarrays = np.concatenate((allarrays, data), axis=2)
    return {'data':allarrays,'transform':ds.GetGeoTransform(),'projection':ds.GetProjection(),'bands':band,'width':row,'height':col}
    # 左上角点坐标 GeoTransform[0],GeoTransform[3] Transform[1] is the pixel width, and Transform[5] is the pixel height

def write_tif(fn_out, im_data, transform,proj=None):
    '''
    功能:
    ----------
    将矩阵按某种投影写入tif，需指定仿射变换矩阵，可选渲染为rgba
    
    参数:
    ----------
    fn_out:str
        输出tif图的绝对文件路径
    im_data: np.array
        tif图对应的矩阵
    transform: list/tuple 
        gdal-like仿射变换矩阵，若im_data矩阵起始点为左上角且投影为4326，则为
            (lon_x.min(), delta_x, 0, 
             lat_y.max(), 0, delta_y)
    proj: str（wkt格式）
        投影，默认投影坐标为4326，可用osr包将epsg转化为wkt格式，如
            srs = osr.SpatialReference()# establish encoding
            srs.ImportFromEPSG(4326)    # WGS84 lat/lon
            proj = srs.ExportToWkt()    # create wkt fromat of proj


    # 设置投影，proj为wkt format
    if proj is None:
        proj = 'GEOGCS["WGS 84",\
                     DATUM["WGS_1984",\
                             SPHEROID["WGS 84",6378137,298.257223563, \
                                    AUTHORITY["EPSG","7030"]], \
                             AUTHORITY["EPSG","6326"]], \
                     PRIMEM["Greenwich",0, \
                            AUTHORITY["EPSG","8901"]], \
                     UNIT["degree",0.0174532925199433, \
                            AUTHORITY["EPSG","9122"]],\
                     AUTHORITY["EPSG","4326"]]'
    # 渲染为rgba矩阵
    # 设置数据类型
    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
    # 将(通道数、高、宽)顺序调整为(高、宽、通道数)
    # print('shape of im data:', im_data.shape)
    im_bands = min(im_data.shape)
    im_shape = list(im_data.shape)
    im_shape.remove(im_bands)
    im_height, im_width = im_shape
    band_idx = im_data.shape.index(im_bands)
    # 找出波段是在第几个

    # 创建文件
    driver = gdal.GetDriverByName("GTiff")
    dataset = driver.Create(fn_out, im_width, im_height, im_bands, datatype)

    # if dataset is not None:
    dataset.SetGeoTransform(transform)  # 写入仿射变换参数
    dataset.SetProjection(proj)  # 写入投影

    if im_bands == 1:

        # print(im_data[:, 0,:].shape)
        if band_idx == 0:
            dataset.GetRasterBand(1).WriteArray(im_data[0, :, :])
        elif band_idx == 1:
            dataset.GetRasterBand(1).WriteArray(im_data[:, 0, :])
        elif band_idx == 2:
            dataset.GetRasterBand(1).WriteArray(im_data[:, :, 0])

    else:

        for i in range(im_bands):
            if band_idx == 0:
                dataset.GetRasterBand(i + 1).WriteArray(im_data[i, :, :])
            elif band_idx == 1:
                dataset.GetRasterBand(i + 1).WriteArray(im_data[:, i, :])
            elif band_idx == 2:
                dataset.GetRasterBand(i + 1).WriteArray(im_data[:, :, i])

    dataset.FlushCache()
    del dataset
    driver = None



def split(origin_data,origin_transform,output_size):
    origin_size = origin_data.shape
    x = origin_transform[0]
    y = origin_transform[3]
    x_step = origin_transform[1]
    y_step = origin_transform[5]
    output_x_step = x_step
    output_y_step = y_step
    for i in range(origin_size[0]//output_size[0]):
        for j in range(origin_size[1]//output_size[1]):
            output_data = origin_data[i*output_size[0]:(i+1)*output_size[0],j*output_size[1]:(j+1)*output_size[1],:]
            output_transform = (x+j*output_x_step*output_size[0],output_x_step,0,y+i*output_y_step*output_size[0],0,output_y_step) 
            write_tif(f'test/{i}_{j}.tif', output_data, output_transform)

data = read_tif('/Users/0_20210508GF2.tif')
split(origin_data=data['data'],origin_transform=data['transform'],output_size=(255,255))

放下最终效果，中心部分就是截取出的小tif，覆盖在原tif之上会有色差，这是因为软件对不同tif文件显示不同色域的原因。