首页 学习区

在使用 GDAL 读取遥感影像的 dateset 后,通过 Projection 和 GeoTransform 这两个参数即可确定影像的地理位置。对于开源的遥感影像数据集,它们的这两个主要参数一般都会进行“脱敏”处理,即去除投影坐标系信息,修饰 GeoTransform 的范围,从而使影像地理坐标位于笛卡尔坐标系下。去除投影信息的结果是显然意见的,那么,对六个 GeoTransform 参数进行修饰会带来怎样的影响呢?一次批量裁剪数据时出现的影响错位问题困扰了我许久,拨开云雾查明原因后,才发现这个问题与这些参数密切相关......

简介 GeoTransform

在遥感图像中,x 与 y 分别对应行与列,原点对应于左上角的像素坐标。使用 gdal 读取 GeoTransform 的方式如下:

from osgeo import gdal
filename = "D:\\Desktop\a.tif"
ds = gdal.Open(filename)
geo_trans = ds.GetGeoTransform()
print(geo_trans)
del ds

对于未做特殊处理的遥感影像,结果一般如下所示:

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

这六个参数的详细解释如下图所示:

geo_trans

其中需留意的是,这里的 geo_trans[5]负数,也正因为此,图像才能从左上角的 (0, 0) 铺展到右下角。如果这个值脱离一般情形,变为正数的话,那么结果会怎样呢?

反面教材

以下为反面实例。如果不调整 geo_trans[5] 的话,裁剪出来的影像将如下图,虽然保留了正确的图像信息,但位置发生了偏移。

image-20220108003111772

打开 ArcGIS 查看原始图像的基本属性:

image-20220108003607665

从中可以看到,坐标系已做了脱敏处理。试着用 GDAL 来读取一下 geotrans 试试:

image-20220108004150619

发现 geo_trans 也进行了脱敏处理。留意这里显示的左上角坐标与 ArcGIS 中显示的有所不同。就是这两个小细节会影响最终结果的整体表现,所以处理数据一定要细心!

改进后的随机裁剪

随机裁剪常用于批量生成训练集,测试算法性能等。一副遥感影像可以被 GDAL 整体读取为 np.array,随机裁剪正可以由此进行。一个简单的思路是从左上方开始,引入 random 库随机在可行的位置定位小窗口的左上角,然后进行裁剪,按坐标偏移量更新 geotrans,最后对裁剪出的窗口设置与原图像相同的投影并导出。

在上述的基础上,我额外进行了如下改进:

  • 引入全局变量记录左上角坐标值以避免裁剪出重复图像;
  • 对脱敏数据额外进行判断,避免裁剪出的图像出现坐标偏移;
  • 结合实际数据加入 +—0.5 的坐标偏离量(非必需,具体问题具体分析);
  • 确保裁剪出图像左上角坐标位置为 20 的整数倍,以匹配后续图像。

核心代码如下,已经写了简要的注释:

from osgeo import gdal
import numpy as np
import random


coords_done = []  # check duplicated situation

def random_crop(ImagePath, ImageSavePath, CropSize, CutNum):
    """
    ImagePath: The input raster to be clipped
    ImageSavePath: The path to save the clipped images
    CropSize: The size of each clipped image (unit: pixel)
    CutNum: The number of exported images with `CropSize`
    """
    dataset_img = readTif(ImagePath)
    width = dataset_img.RasterXSize
    height = dataset_img.RasterYSize
    proj = dataset_img.GetProjection()
    geotrans = dataset_img.GetGeoTransform()
    img = dataset_img.ReadAsArray(0, 0, width, height)
    
    geotrans_list = []
    geotrans_list.append(geotrans[0])
    geotrans_list.append(geotrans[1])
    geotrans_list.append(geotrans[2])
    geotrans_list.append(geotrans[3])
    geotrans_list.append(geotrans[4])
    if geotrans[5] < 0:
        geotrans_list.append(geotrans[5]*-1)
    else:
        geotrans_list.append(geotrans[5])

    if not os.path.exists(ImageSavePath):
        os.makedirs(ImageSavePath)

    fileNum = len(os.listdir(ImageSavePath))
    new_name = fileNum + 1
    
    global coords_done
    
    while (new_name < CutNum + fileNum + 1):
        # ensure the left upper coordinate is k*20 (k is integeral)
        UpperLeftX = 20 * random.randint(0, (height - CropSize) // 20)
        UpperLeftY = 20 * random.randint(0, (width - CropSize) // 20)
        if (UpperLeftX, UpperLeftY) not in coords_done:
            coords_done.append((UpperLeftX, UpperLeftY))
        else:
            print("Image Duplicate! Skip to next clipped image...")
            continue
        # When there is only one band, then `img.shape` will be 2 rather than 3
        if (len(img.shape) == 2):
            imgCrop = img[UpperLeftX: UpperLeftX + CropSize,
                      UpperLeftY: UpperLeftY + CropSize]
        else:
            imgCrop = img[:,
                      UpperLeftX: UpperLeftX + CropSize,
                      UpperLeftY: UpperLeftY + CropSize]
        x = int(UpperLeftX)  # row
        y = int(UpperLeftY)  # col
        # calculate new position
        px = geotrans[0] + y * geotrans[1] + x * geotrans[2]
        if geotrans[5] < 0:
            py = geotrans[3] + y * geotrans[4] + x * geotrans[5] * -1
        else:
            py = geotrans[3] + y * geotrans[4] + x * geotrans[5]
        # The +- 0.5 is needed to match the image accurately
        geotrans_list[0] = px - 0.5
        geotrans_list[3] = py + 0.5
        geotrans_tuple = tuple(geotrans_list)
        writeTiff(imgCrop, geotrans_tuple, proj, ImageSavePath + "/%d.tif" % new_name)
        new_name = new_name + 1
    coords_done = []

最终取得了心仪的结果(原始底图调整了波段组合与透明度以强调对比):

image-20220108005117915




文章评论

目录