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 os

filename = "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.float16
    height = dataset.RasterYSize
    width = dataset.RasterXSize
    projection = 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_image

image_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 = 7
        dtype = gdal.GDT_Float32
        DIMENSION_OF_IMAGE = 3
        if(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 image
        del ds_to_save
write(save_path1, image_part1, projection, geotransform)
write(save_path2, image_part2, projection, new_geotransform)
结果

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

### 使用GDAL库进行地理空间数据分析 为了利用Python中的GDAL库执行地理空间数据处理,安装必要的软件包是首要条件。可以从指定链接下载所需的代码和样本数据[^1]。这些资源不仅包含了书本实例还提供了自定义工具以及所有用于示例的数据集。 一旦环境准备就绪,在Python中操作GDAL主要依赖于`osgeo.gdal`模块来读取、写入栅格文件并获取元数据信息;而矢量数据的操作则通过`osgeo.ogr`完成。下面是一个简单的例子展示如何打开一个GeoTIFF图像文件: ```python from osgeo import gdal dataset = gdal.Open('path_to_file.tif') if dataset is None: print("Unable to open file") else: print(f"Driver: {dataset.GetDriver().ShortName}/{dataset.GetDriver().LongName}") print(f"Size is {dataset.RasterXSize} x {dataset.RasterYSize} x {dataset.RasterCount}") print(f"Projection is {dataset.GetProjection()}") geotransform = dataset.GetGeoTransform() if not geotransform is None: print(f"Origin = ({geotransform[0]}, {geotransform[3]})") print(f"Pixel Size = ({geotransform[1]}, {-geotransform[5]})") band = dataset.GetRasterBand(1) print(f"Band Type={gdal.GetDataTypeName(band.DataType)}") min_value = band.ComputeRasterMinMax()[0] max_value = band.ComputeRasterMinMax()[1] print(f"Min={min_value}, Max={max_value}") scanline = band.ReadRaster(xoff=0, yoff=0, xsize=band.XSize, ysize=band.YSize, buf_xsize=band.XSize, buf_ysize=band.YSize, buf_type=gdal.GDT_Float32) import struct tuple_of_floats = struct.unpack('f' * band.XSize*band.YSize, scanline) ``` 此脚本展示了基本功能,比如加载影像、访问其属性(尺寸大小、投影坐标系)、查询单波段统计特性(最小最大值),并且能够逐像素遍历整个图层的内容。 对于更复杂的任务,如重采样、裁剪、镶嵌等,则可能需要用到额外的功能函数或是第三方扩展库的支持。例如,可以借助`rasterio`这样的高级接口简化上述过程,或者采用专门针对特定应用领域开发的附加组件来进行更加专业的分析工作。 此外,一些辅助性的Python库可以帮助更好地管理和可视化地理空间数据,如Pandas用于结构化表格型数据管理[^2],Matplotlib或Seaborn提供图形绘制能力,Missingno用来检测缺失值情况等等。
评论 4
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值