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

Question

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

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

我可以使用不同的波段组合创建新的伪彩色栅格并将其保存，并且可以使用dstack（第一块代码）在python中创建所需的IMAGE，但是我无法将该图像另存为georectified。 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倍数为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可以使您摆脱所有标准文件创建样板，您仍然需要考虑数据类型，因为与输入文件相比，数据类型可能会发生变化。