从标准偏差拉伸带创建新栅格（.tif），与 dstack 一起使用，但不编写新文件，Python

Question

如果标题不清楚，我很抱歉，我是python新手，词汇量有限。

我要做的是对 .tif 栅格中的每个波段应用标准偏差拉伸，然后通过使用 GDAL (Python) 堆叠这些波段来创建一个新的栅格 (.tif)。

我能够创建具有不同波段组合的新假色光栅并保存它们，并且我能够使用 dstack（第一个代码块）在 python 中创建我想要的图像，但我无法将该图像保存为地理校正.tif 文件。

所以要使用 dstack 创建拉伸图像，我的代码如下所示：

import os
import numpy as np
import matplotlib.pyplot as plt
import math
from osgeo import gdal

# code from my prof
def std_stretch_data(data, n=2):
    """Applies an n-standard deviation stretch to data."""

    # Get the mean and n standard deviations.
    mean, d = data.mean(), data.std() * n

    # Calculate new min and max as integers. Make sure the min isn't
    # smaller than the real min value, and the max isn't larger than
    # the real max value.
    new_min = math.floor(max(mean - d, data.min()))
    new_max = math.ceil(min(mean + d, data.max()))

    # Convert any values smaller than new_min to new_min, and any
    # values larger than new_max to new_max.
    data = np.clip(data, new_min, new_max)

    # Scale the data.
    data = (data - data.min()) / (new_max - new_min)
    return data

# open the raster
img = gdal.Open(r'/Users/Rebekah/ThesisData/TestImages/OG/OG_1234.tif')

#open the bands
red = img.GetRasterBand(1).ReadAsArray()
green = img.GetRasterBand(2).ReadAsArray()
blue = img.GetRasterBand(3).ReadAsArray()

# create alpha band where a 0 indicates a transparent pixel and 1 is a opaque pixel
# (this is from class and i dont FULLY understand it)
alpha = np.where(red + green + blue ==0, 0, 1).astype(np.byte)

red_stretched = std_stretch_data(red, 1)
green_stretched = std_stretch_data(green, 1)
blue_stretched = std_stretch_data(blue, 1)

data_stretched = np.dstack((red_stretched, green_stretched, blue_stretched, alpha))
plt.imshow(data_stretched)
plt.show()

这让我在一个单独的窗口中准确地看到了我想要的东西。但该代码中没有任何位置可以分配投影，或将其保存为多波段 tif。

所以我采用了它并尽我所能将其应用于我用来创建假彩色图像的代码，但它失败了（代码如下）。如果我使用 alpha 波段创建 4 波段 tif，则输出为空 tif，如果我创建 3 波段 tif 并省略 alpha 波段，则输出是完全黑色的 tif。

import os
import numpy as np
import matplotlib.pyplot as plt
import math
from osgeo import gdal

#code from my professor
def std_stretch_data(data, n=2):
    """Applies an n-standard deviation stretch to data."""

    # Get the mean and n standard deviations.
    mean, d = data.mean(), data.std() * n

    # Calculate new min and max as integers. Make sure the min isn't
    # smaller than the real min value, and the max isn't larger than
    # the real max value.
    new_min = math.floor(max(mean - d, data.min()))
    new_max = math.ceil(min(mean + d, data.max()))

    # Convert any values smaller than new_min to new_min, and any
    # values larger than new_max to new_max.
    data = np.clip(data, new_min, new_max)

    # Scale the data.
    data = (data - data.min()) / (new_max - new_min)
    return data

#open image
img = gdal.Open(r'/Users/Rebekah/ThesisData/TestImages/OG/OG_1234.tif')

# get geotill driver
gtiff_driver = gdal.GetDriverByName('GTiff')

# read in bands
red = img.GetRasterBand(1).ReadAsArray()
green = img.GetRasterBand(2).ReadAsArray()
blue = img.GetRasterBand(3).ReadAsArray()

# create alpha band where a 0 indicates a transparent pixel and 1 is a opaque pixel
# (this is from class and i dont FULLY understand it)
alpha = np.where(red + green + blue ==0, 0, 1).astype(np.byte)

# apply the 1 standard deviation stretch
red_stretched = std_stretch_data(red, 1)
green_stretched = std_stretch_data(green, 1)
blue_stretched = std_stretch_data(blue, 1)

# create empty tif file
NewImg = gtiff_driver.Create('/Users/riemann/ThesisData/TestImages/FCI_tests/1234_devst1.tif', img.RasterXSize, img.RasterYSize, 4, gdal.GDT_Byte)
if NewImg is None:
    raise IOerror('could not create new raster')

# set the projection and geo transform of the new raster to be the same as the original
NewImg.SetProjection(img.GetProjection())
NewImg.SetGeoTransform(img.GetGeoTransform())

# write new bands to the new raster
band1 = NewImg.GetRasterBand(1)
band1.WriteArray(red_stretched)

band2 = NewImg.GetRasterBand(2)
band2.WriteArray(green_stretched)

band3= NewImg.GetRasterBand(3)
band3.WriteArray(blue_stretched)

alpha_band = NewImg.GetRasterBand(4)
alpha_band.WriteArray(alpha)

del band1, band2, band3, img, alpha_band

我不完全确定如何从这里开始并创建一个新文件来显示不同频段的拉伸。

该图像只是从 earthexplorer 下载的 4 波段光栅 (NAIP)，如果需要，我可以上传用于测试的特定图像，但与其他 NAIP 图像相比，该文件本身并没有什么特别之处。

Answer 1

您也应该关闭新数据集 (NewImg)，方法是将其添加到您已有的 del 列表中，或将其设置为 None。

这会正确关闭文件并确保所有数据都写入磁盘。

但是还有另一个问题，您在 0 和 1 之间缩放数据，但将其存储为 Byte。所以要么将输出数据类型从gdal.GDT_Byte 更改为gdal.GDT_Float32。或者将缩放后的数据相乘以适应输出数据类型，在 Byte multiple 为 255 的情况下（不要忘记 alpha），您应该正确舍入以确保准确性，否则 GDAL 将截断为最接近的整数。

您可以使用np.iinfo() 检查数据类型的范围，以防您不确定要为其他数据类型使用什么乘法。

根据您的用例，使用gdal.Translate 进行缩放可能是最简单的。如果您要稍微修改缩放函数以返回缩放参数而不是数据，您可以使用类似：

ds = gdal.Translate(output_file, input_file, outputType=gdal.GDT_Byte, scaleParams=[
    [old_min_r, old_max_r, new_min_r, new_max_r], # red
    [old_min_g, old_max_g, new_min_g, new_max_g], # green
    [old_min_b, old_max_b, new_min_b, new_max_b], # blue
    [old_min_a, old_max_a, new_min_a, new_max_a], # alpha
])
ds = None

您还可以为非线性拉伸添加exponents 关键字。

使用gdal.Translate 将使您免于创建所有标准文件的样板文件，但您仍然需要考虑数据类型，因为与输入文件相比，它可能会发生变化。