GDAL(python) 之GeoTransform（3）

使用GDAL库读出的dataset带有两个重要的地理参数，分别是Projection（暂未搞清这个是干啥的）和GeoTransform（用作仿射变换的）。有了这两个参数，就确定了影像的地理位置。

再GDAL for Python中，GeoTransform是一个六个元素的元组。
例如加载一个影像，读取并显示它的GeoTransform，则为如下形式：

(486892.5, 15.0, 0.0, 4105507.5, 0.0, -15.0)

六个参数分别为：

左上角x坐标， 水平分辨率，旋转参数， 左上角y坐标，旋转参数，竖直分辨率（负值，原点坐标在左上角）。
上代码：

注意：本代码仅测试输出的top和botton可以上下接起来，并不讨论top.tif和bottopm.tif的真实投影位置。

# coding=utf-8
import numpy as np
from osgeo import gdal, gdal_array
import os


#os.chdir(r'D:\DeskTop\learn_py_must\Learn_GDAL\osgeopy-data\osgeopy-data\Washington\dem')

filename = r"D:\DeskTop\learn_py_must\Learn_GDAL\osgeopy-data\osgeopy-data\Landsat\Washington\nat_color.tif"
dataset = gdal.Open(filename, gdal.GA_ReadOnly)
if dataset == None:
    raise Exception("Image name error.")
else:
    datatype = np.float16
    height = dataset.RasterYSize
    width = dataset.RasterXSize
    projection = dataset.GetProjection()
    print("projection:", projection)
    geotransform = dataset.GetGeoTransform()
    print("geotransform:",geotransform)

    band_image = np.zeros((height, width, 1), dtype=datatype)
    print(band_image.shape)
    band_data = dataset.GetRasterBand(1)
    band_image[:, :, 0] = band_data.ReadAsArray()
    del dataset
test_image = band_image[1000:3000, 1000:3000]
print(test_image.shape)
del band_image

image_part1 = test_image[:1000, :2000]  # 取前一千行
print(image_part1.shape)
image_part2 = test_image[1000:2000, :2000]  # 取1000 ~ 2000行
print(image_part2.shape)
new_x_geo = geotransform[0] + geotransform[1] * 0  # 新横坐标起始量
new_y_geo = geotransform[3] + geotransform[5] * 1000  # 新纵坐标起始量
new_geotransform = (new_x_geo, geotransform[1], geotransform[2], new_y_geo, geotransform[4], geotransform[5])

save_path1 = "D:\\DeskTop\\test\\Top.tif"
save_path2 = "D:\\DeskTop\\test\\bottom.tif"


def write(save_path, image, projection, geotransform, format='GTiff'):

    dtype = gdal.GDT_Float32


    height = image.shape[0]
    width = image.shape[1]


    driver = gdal.GetDriverByName(format)
    ds_to_save = driver.Create(save_path, width, height, dtype)
    ds_to_save.SetGeoTransform(geotransform)
    ds_to_save.SetProjection(projection)

    for band in range(1):
        ds_to_save.GetRasterBand(band + 1).WriteArray(image[:, :, band])
        ds_to_save.FlushCache()

    del image
    del ds_to_save


write(save_path1, image_part1, projection, geotransform)
write(save_path2, image_part2, projection, new_geotransform)

参考链接：GDAL(python) 之GeoTransform-CSDN博客