Python遥感影像分块（gdal）

分享一下利用python的gdal库把大的遥感影像分成小块（如500*500）的代码

import numpy as np
from osgeo import gdal
import pandas as pd
import os
import re
def get_year_from_filename(filename):
    """
    Extracts the year from the filename using regex
    """
    match = re.search(r'\d{4}', filename)
    return match.group(0) if match else None

def split_image(image_path, output_dir, tile_size, year):
    """
    Splits the image into smaller tiles of size tile_size x tile_size
    """
    ds = gdal.Open(image_path)
    band_count = ds.RasterCount
    width = ds.RasterXSize
    height = ds.RasterYSize
    geotransform = ds.GetGeoTransform()
    projection = ds.GetProjection()

    for i in range(0, width, tile_size):
        for j in range(0, height, tile_size):
            w = min(tile_size, width - i)
            h = min(tile_size, height - j)
            tile = ds.ReadAsArray(i, j, w, h)
            tile_output_dir = os.path.join(output_dir, f'tile_{i}_{j}')
            os.makedirs(tile_output_dir, exist_ok=True)
            tile_output_path = os.path.join(tile_output_dir, f'{year}_tile_{i}_{j}.tif')
            driver = gdal.GetDriverByName('GTiff')
            out_ds = driver.Create(tile_output_path, w, h, band_count, gdal.GDT_Float32)
            out_ds.SetGeoTransform((geotransform[0] + i * geotransform[1], geotransform[1], 0, geotransform[3] + j * geotransform[5], 0, geotransform[5]))
            out_ds.SetProjection(projection)
            for k in range(band_count):
                out_ds.GetRasterBand(k + 1).WriteArray(tile[k])
            out_ds.FlushCache()
            out_ds = None

def main(input_directory, output_directory, tile_size):
    # Get all image files in the input directory
    image_files = [f for f in os.listdir(input_directory) if f.endswith('.tif')]
    
    # Process each image file
    for image_file in image_files:
        year = get_year_from_filename(image_file)
        if year:
            image_path = os.path.join(input_directory, image_file)
            split_image(image_path, output_directory, tile_size, year)
            print(f'Processed image {image_file} for year {year}')

需要改的只有文件夹路径和切片的大小：

if __name__ == "__main__":
    input_directory = "Long_time_Landsat_dataset"
    output_directory = "Split_zone"
    tile_size = 4000  # or 500
    main(input_directory, output_directory, tile_size)