根据标准差拉伸波段创建新栅格(.tif)，适用于数据堆栈，但不能写入新文件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

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

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

Answer 1

您还应该通过将新数据集(NewImg)添加到已有的del列表中或将其设置为None来关闭它。

这将正确地关闭该文件，并确保所有数据都写入磁盘。

然而，还有另一个问题，您在0和1之间缩放数据，但将其存储为Byte。因此，要么将输出数据类型从gdal.GDT_Byte更改为类似于gdal.GDT_Float32的类型。或者将缩放数据相乘以适合输出数据类型，在字节倍数与255 (不要忘记字母)的情况下，您应该正确地对其进行舍入以提高精度，否则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将使您从所有标准文件创建样板文件中解脱出来，您仍然需要考虑数据类型，因为与输入文件相比，数据类型可能会发生变化。