GDAL(python) 之GeoTransform_综合

演示

使用GDAL库读出的dataset带有两个重要的地理参数，分别是Projection和GeoTransform。有了这两个参数，就确定了影像的地理位置。

再GDAL for Python中，GeoTransform是一个六个元素的元组。
例如，我找了一个影像，读取并显示它的GeoTransform，则为如下形式：
在这里插入图片描述

形式

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

六个参数分别为：

左上角x坐标，水平分辨率，旋转参数，左上角y坐标，旋转参数，竖直分辨率。如另一篇博客中所说，满足如下关系式：
在这里插入图片描述
一般来说，旋转参数都为0。

实践

有了上述知识，不妨利用一张图片进行实际的验证。

总体的思路是，首先读取我下载的某Landsat8数据集的全色波段，然后截取一个2000 * 2000的部分。

然后我们手动的把图片分成上下两部分，然后给上下两部分按照他们的情况，定义新的GeoTransform，再将图片写出，看看我们定义的GeoTransform是否正确，以此验证我们的知识。

下面是代码验证：

import numpy as np
from osgeo import gdal, gdal_array
import osfilename = "F:\SRGAN_program\dataset\LC81290352019095LGN00\LC08_L1TP_129035_20190405_20190422_01_T1_B8.tif"
dataset = gdal.Open(filename, gdal.GA_ReadOnly)
if dataset == None:raise Exception("Image name error.")
else:datatype = np.float16height = dataset.RasterYSizewidth = dataset.RasterXSizeprojection = dataset.GetProjection()geotransform = dataset.GetGeoTransform()band_image = np.zeros((height, width, 1), dtype = datatype)band_data = dataset.GetRasterBand(1)band_image[:, :, 0] = band_data.ReadAsArray()del dataset
test_image = band_image[4000:6000, 4000:6000, :]
del band_imageimage_part1 = test_image[:1000, :2000, :] # 取前一千行
image_part2 = test_image[1000:2000, :2000, :] # 取1000 ~ 2000行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 = "F:\\SRGAN_program\\dataset\\LC81290352019095LGN00\\test\\top"
save_path2 = "F:\\SRGAN_program\dataset\\LC81290352019095LGN00\\test\\bottom"def write(save_path, image, projection, geotransform, format = 'ENVI'):LANDSAT8_MULTI_BAND = 7dtype = gdal.GDT_Float32DIMENSION_OF_IMAGE = 3if(len(image.shape) != DIMENSION_OF_IMAGE):raise Exception("The dimension of the image is incorrect.") else:height = image.shape[0]width = image.shape[1]channels = image.shape[2]driver = gdal.GetDriverByName(format)ds_to_save = driver.Create(save_path, width, height, channels, 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 imagedel ds_to_save
write(save_path1, image_part1, projection, geotransform)
write(save_path2, image_part2, projection, new_geotransform)

结果

我们将图片分成了上下两部分，并写入了我们定义的新的地理参考。下面从ENVI打开来显示一下，看我们定义的是否正确。
下半部分：
在这里插入图片描述
上半部分：

总体显示：