[已解决]请大佬帮忙看看这个问题怎么解决 谢谢

Justheartyoung

下面代码是想对无人机影像（RGB, TIff格式）处理，提取得到影像中的沙丘脊线，但是在对提取结果保存为带有和原始影像一样的投影坐标栅格结果后，用Arcgis软件打开是一个全黑的（结果图放在了下面）。是不是将结果保存时出了问题，请大佬帮我改写为正确的代码！
from tifffile import imread
import numpy as np
from osgeo import gdal, osr

# 读取TIFF格式无人机影像数据
image_path = 'F:\\duneline\\dune\\dune.tif'
output_path = 'overlay_image.tif'

image = imread(image_path)

# 将输入图像转换为灰度图像
gray_image = np.mean(image, axis=2)

# 自适应阈值二值化
binary = np.zeros_like(gray_image, dtype=np.uint8)
window_size = 51 # 窗口大小
k = 0.3 # 控制阈值的参数
for i in range(window_size//2, gray_image.shape[0]-window_size//2):
    for j in range(window_size//2, gray_image.shape[1]-window_size//2):
        window = gray_image[i-window_size//2:i+window_size//2+1, j-window_size//2:j+window_size//2+1]
        threshold = np.mean(window) - k * np.std(window)
        if gray_image[i, j] > threshold:
            binary[i, j] = 255

# 配置输出的空间参考信息
in_ds = gdal.Open(image_path)
projection = in_ds.GetProjection()
geotransform = in_ds.GetGeoTransform()

# 创建新的输出图像
driver = gdal.GetDriverByName('GTiff')
out_ds = driver.Create(output_path, image.shape[1], image.shape[0], 1, gdal.GDT_UInt16)

# 设置输出图像的空间参考信息
out_ds.SetProjection(projection)
out_ds.SetGeoTransform(geotransform)

# 将沙脊线部分写入输出图像
out_band = out_ds.GetRasterBand(1)
out_band.WriteArray(binary * 255)  # 将二值图像数据缩放到0-255的范围

# 设置颜色表
color_table = gdal.ColorTable()
color_table.SetColorEntry(0, (0, 0, 0, 0))  # 黑色，对应值为0
color_table.SetColorEntry(255, (255, 255, 255, 255))  # 白色，对应值为255
out_band.SetColorTable(color_table)

# 释放资源
out_band.FlushCache()
out_ds = None
in_ds = None

print("保存成功！")
结果：

最佳答案

月排行榜 / 总排行榜

isdkz

2023-8-5 13:34:02

你的代码看起来是正确的，但是有一点需要注意。你在二值化处理后，已经将binary数组的值设置为255，然后在写入输出图像时，你又将binary数组的值乘以255，这就导致了值的溢出，因此在Arcgis软件中打开时看到的是全黑的图像。

你应该直接将binary数组写入输出图像，而不是binary * 255。另外，你的输出图像数据类型设置为了gdal.GDT_UInt16，这是16位无符号整数，但是你的binary数组是8位无符号整数，所以你应该将输出图像的数据类型设置为gdal.GDT_Byte。


以下是修改后的代码：

from tifffile import imread
import numpy as np
from osgeo import gdal, osr


# 读取TIFF格式无人机影像数据
image_path = 'F:\duneline\dune\dune.tif'
output_path = 'overlay_image.tif'


image = imread(image_path)


# 将输入图像转换为灰度图像
gray_image = np.mean(image, axis=2)


# 自适应阈值二值化
binary = np.zeros_like(gray_image, dtype=np.uint8)
window_size = 51 # 窗口大小
k = 0.3 # 控制阈值的参数
for i in range(window_size//2, gray_image.shape[0]-window_size//2):
    for j in range(window_size//2, gray_image.shape[1]-window_size//2):
        window = gray_image[i-window_size//2:i+window_size//2+1, j-window_size//2:j+window_size//2+1]
        threshold = np.mean(window) - k * np.std(window)
        if gray_image[i, j] > threshold:
            binary[i, j] = 255


# 配置输出的空间参考信息
in_ds = gdal.Open(image_path)
projection = in_ds.GetProjection()
geotransform = in_ds.GetGeoTransform()


# 创建新的输出图像
driver = gdal.GetDriverByName('GTiff')
out_ds = driver.Create(output_path, image.shape[1], image.shape[0], 1, gdal.GDT_Byte)  # 修改数据类型为gdal.GDT_Byte


# 设置输出图像的空间参考信息
out_ds.SetProjection(projection)
out_ds.SetGeoTransform(geotransform)


# 将沙脊线部分写入输出图像
out_band = out_ds.GetRasterBand(1)
out_band.WriteArray(binary)  # 直接写入binary数组


# 设置颜色表
color_table = gdal.ColorTable()
color_table.SetColorEntry(0, (0, 0, 0, 0))  # 黑色，对应值为0
color_table.SetColorEntry(255, (255, 255, 255, 255))  # 白色，对应值为255
out_band.SetColorTable(color_table)


# 释放资源
out_band.FlushCache()
out_ds = None
in_ds = None


print("保存成功！")

希望这个修改能解决你的问题。

跳转到最佳答案楼层

学习编程中的Ben

问题出在将二值图像保存到输出文件时，没有正确地进行数据类型转换和范围缩放。

在这里，将二值图像乘以255的操作是正确的，但是需要将数据类型转换为无符号字节型（gdal.GDT_Byte）。此外，还需要设置输出图像的颜色表为灰度色彩映射。

请修改以下代码段：

# 将沙脊线部分写入输出图像out_band = out_ds.GetRasterBand(1)out_band.WriteArray((binary * 255).astype(np.uint8))  # 将二值图像数据缩放到0-255的范围并转换为无符号字节型# 设置颜色表color_table = gdal.ColorTable()for i in range(256):    color_table.SetColorEntry(i, (i, i, i, 255))  # 灰度色彩映射out_band.SetColorTable(color_table)

这样，代码会将二值图像数据转换为无符号字节型，并设置输出图像的颜色表为灰度色彩映射。然后，保存输出图像后，用ArcGIS软件打开即可正确显示沙脊线部分的图像。

注意：由于无法访问您的文件系统，上述代码中的output_path和image_path变量需要替换为您自己的路径。

希望这可以帮助您解决问题！如果有任何其他疑问，请随时提问。

如果回答对你有帮助，请给我一个最佳答案！

isdkz

这个问题可能是由于保存结果时出了问题导致的。你可以尝试修改以下代码：

1. 在创建输出图像时，将数据类型设置为gdal.GDT_Byte，而不是gdal.GDT_UInt16。因为二值图像只有0和255两个值，使用8位整数类型即可。

2. 在设置颜色表之前，将输出图像的颜色解释设置为颜色表。可以使用out_band.SetRasterColorInterpretation(gdal.GCI_PaletteIndex)来实现。

修改后的代码如下：

from tifffile import imread
import numpy as np
from osgeo import gdal, osr

# 读取TIFF格式无人机影像数据
image_path = 'F:\\duneline\\dune\\dune.tif'
output_path = 'overlay_image.tif'

image = imread(image_path)

# 将输入图像转换为灰度图像
gray_image = np.mean(image, axis=2)

# 自适应阈值二值化
binary = np.zeros_like(gray_image, dtype=np.uint8)
window_size = 51 # 窗口大小
k = 0.3 # 控制阈值的参数
for i in range(window_size//2, gray_image.shape[0]-window_size//2):
    for j in range(window_size//2, gray_image.shape[1]-window_size//2):
        window = gray_image[i-window_size//2:i+window_size//2+1, j-window_size//2:j+window_size//2+1]
        threshold = np.mean(window) - k * np.std(window)
        if gray_image[i, j] > threshold:
            binary[i, j] = 255

# 配置输出的空间参考信息
in_ds = gdal.Open(image_path)
projection = in_ds.GetProjection()
geotransform = in_ds.GetGeoTransform()

# 创建新的输出图像
driver = gdal.GetDriverByName('GTiff')
out_ds = driver.Create(output_path, image.shape[1], image.shape[0], 1, gdal.GDT_Byte)

# 设置输出图像的空间参考信息
out_ds.SetProjection(projection)
out_ds.SetGeoTransform(geotransform)

# 将沙脊线部分写入输出图像
out_band = out_ds.GetRasterBand(1)
out_band.WriteArray(binary * 255)  # 将二值图像数据缩放到0-255的范围

# 设置颜色表
color_table = gdal.ColorTable()
color_table.SetColorEntry(0, (0, 0, 0, 0))  # 黑色，对应值为0
color_table.SetColorEntry(255, (255, 255, 255, 255))  # 白色，对应值为255
out_band.SetRasterColorInterpretation(gdal.GCI_PaletteIndex)
out_band.SetColorTable(color_table)

# 释放资源
out_band.FlushCache()
out_ds = None
in_ds = None

print("保存成功！")

希望这个修改能够解决你的问题！如果还有其他问题，请随时提问。

Justheartyoung

isdkz 发表于 2023-8-5 13:19
这个问题可能是由于保存结果时出了问题导致的。你可以尝试修改以下代码：

1. 在创建输出图像时，将数据 ...

不行，结果还是那样，全黑的

isdkz

你的代码看起来是正确的，但是有一点需要注意。你在二值化处理后，已经将binary数组的值设置为255，然后在写入输出图像时，你又将binary数组的值乘以255，这就导致了值的溢出，因此在Arcgis软件中打开时看到的是全黑的图像。

你应该直接将binary数组写入输出图像，而不是binary * 255。另外，你的输出图像数据类型设置为了gdal.GDT_UInt16，这是16位无符号整数，但是你的binary数组是8位无符号整数，所以你应该将输出图像的数据类型设置为gdal.GDT_Byte。


以下是修改后的代码：

from tifffile import imread
import numpy as np
from osgeo import gdal, osr


# 读取TIFF格式无人机影像数据
image_path = 'F:\duneline\dune\dune.tif'
output_path = 'overlay_image.tif'


image = imread(image_path)


# 将输入图像转换为灰度图像
gray_image = np.mean(image, axis=2)


# 自适应阈值二值化
binary = np.zeros_like(gray_image, dtype=np.uint8)
window_size = 51 # 窗口大小
k = 0.3 # 控制阈值的参数
for i in range(window_size//2, gray_image.shape[0]-window_size//2):
    for j in range(window_size//2, gray_image.shape[1]-window_size//2):
        window = gray_image[i-window_size//2:i+window_size//2+1, j-window_size//2:j+window_size//2+1]
        threshold = np.mean(window) - k * np.std(window)
        if gray_image[i, j] > threshold:
            binary[i, j] = 255


# 配置输出的空间参考信息
in_ds = gdal.Open(image_path)
projection = in_ds.GetProjection()
geotransform = in_ds.GetGeoTransform()


# 创建新的输出图像
driver = gdal.GetDriverByName('GTiff')
out_ds = driver.Create(output_path, image.shape[1], image.shape[0], 1, gdal.GDT_Byte)  # 修改数据类型为gdal.GDT_Byte


# 设置输出图像的空间参考信息
out_ds.SetProjection(projection)
out_ds.SetGeoTransform(geotransform)


# 将沙脊线部分写入输出图像
out_band = out_ds.GetRasterBand(1)
out_band.WriteArray(binary)  # 直接写入binary数组


# 设置颜色表
color_table = gdal.ColorTable()
color_table.SetColorEntry(0, (0, 0, 0, 0))  # 黑色，对应值为0
color_table.SetColorEntry(255, (255, 255, 255, 255))  # 白色，对应值为255
out_band.SetColorTable(color_table)


# 释放资源
out_band.FlushCache()
out_ds = None
in_ds = None


print("保存成功！")

希望这个修改能解决你的问题。