非AI大佬 能不能帮忙看下输出结果的栅格值为什么那么奇怪 谢谢

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 = 131 # 窗口大小
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("保存成功！")
结果值：

isdkz

输出结果的栅格值奇怪的原因是代码中没有对二值化后的图像进行数据类型转换。在创建新的输出图像时，应该将数据类型设置为gdal.GDT_Byte，即8位无符号整数。而在写入沙脊线部分时，应该将二值化后的图像数组binary转换为与输出图像相同的数据类型。修改代码如下：

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 = 131 # 窗口大小
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.astype(np.uint8))  # 将二值化后的图像数组转换为与输出图像相同的数据类型

# 设置颜色表
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

我已经对你提供的代码进行修改，以解决输出栅格值异常的问题。以下是修改后的代码：

from tifffile import imreadimport numpy as npfrom 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 = 131  # 窗口大小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)# 设置颜色表color_table = gdal.ColorTable()color_table.SetColorEntry(0, (0, 0, 0, 0))  # 黑色，对应值为0color_table.SetColorEntry(255, (255, 255, 255, 255))  # 白色，对应值为255out_band.SetColorTable(color_table)# 释放资源out_band.FlushCache()out_ds = Nonein_ds = Noneprint("保存成功！")

以上是修改后的完整代码。我注意到你在原始代码中使用了反斜杠\作为文件路径的分隔符，但在Python中应该使用正斜杠/或双反斜杠\\。我将路径中的反斜杠替换为正斜杠以进行修正。

请尝试使用这段修改后的代码，并检查输出栅格的值是否符合预期。如果还有其他问题，请随时提问。

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

学习编程中的Ben

主要是不是AI我们不会啊

Justheartyoung

isdkz 发表于 2023-8-5 15:13
输出结果的栅格值奇怪的原因是代码中没有对二值化后的图像进行数据类型转换。在创建新的输出图像时，应该将 ...

我原本就是gdal.GDT_Byte啊，换成astype(np.uint8)还是不行

Justheartyoung

学习编程中的Ben 发表于 2023-8-5 15:14
主要是不是AI我们不会啊

您不会吗 大佬 这是您给的代码

学习编程中的Ben

Justheartyoung 发表于 2023-8-5 15:17
您不会吗 大佬 这是您给的代码

AI呀!主要是我们技术再好也不可能懂得每一个库

学习编程中的Ben

Justheartyoung 发表于 2023-8-5 15:16
我原本就是gdal.GDT_Byte啊，换成astype(np.uint8)还是不行

非常抱歉，我之前的回答有误。从您提供的代码来看，您在创建输出图像时确实已经正确地将数据类型设置为gdal.GDT_Byte，并且在写入沙脊线部分也进行了数据类型转换使用astype(np.uint8)。

根据您的描述，如果即使在进行数据类型转换后仍然存在问题，可能是其他原因导致的。以下是一些可能导致输出结果错误的常见问题和建议的解决方法：

数据范围不正确：请确保输入的灰度图像、二值化后的图像以及最终的输出图像的数据范围都在0到255之间。可以尝试通过打印相关数组的最小值和最大值来检查数据范围是否正确。

影像坐标系或投影信息不匹配：确认输入图像和输出图像的空间参考信息是否一致，包括投影和地理转换参数。

颜色表设置有误：在设置颜色表时，请确保使用的颜色表与数据类型一致。在您的代码中，数据类型是gdal.GDT_Byte，对应的颜色表应该使用SetColorEntry设置8位无符号整数。

其他错误：如果以上方法都没有解决问题，可能需要更详细地检查代码逻辑、调试代码或查看完整的错误信息来确定问题的具体原因。

希望这些建议能对您有所帮助。如果问题仍然存在，请提供更多信息，我将尽力协助您解决问题。谢谢！

学习编程中的Ben

Justheartyoung 发表于 2023-8-5 15:17
您不会吗 大佬 这是您给的代码

我告诉你,前三个回答的一定都是AI

Justheartyoung

学习编程中的Ben 发表于 2023-8-5 15:22
我告诉你,前三个回答的一定都是AI

那你能做第四个 帮我回答一下吗

学习编程中的Ben

Justheartyoung 发表于 2023-8-5 15:23
那你能做第四个 帮我回答一下吗

em.....你这些库我没用过QAQ
不熟悉啊

Justheartyoung

学习编程中的Ben 发表于 2023-8-5 15:25
em.....你这些库我没用过QAQ
不熟悉啊

哈哈 好的大佬

sfqxx

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 = 131 # 窗口大小
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, binary.shape[1], binary.shape[0], 1, gdal.GDT_Byte,
                      options=['COMPRESS=NONE'])  # 修改数据类型为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()
for i in range(256):  # 对所有可能的像素值进行映射
    if i == 0:
        color_table.SetColorEntry(i, (0, 0, 0, 255))  # 黑色，对应值为0
    elif i == 255:
        color_table.SetColorEntry(i, (255, 255, 255, 255))  # 白色，对应值为255
    else:
        color_table.SetColorEntry(i, (i, i, i, 255))  # 将其他的值映射为相应的灰度色
out_band.SetColorTable(color_table)


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


print("保存成功！")

Justheartyoung

sfqxx 发表于 2023-8-5 17:35

好的大佬  谢谢您 我去试试

Justheartyoung

sfqxx 发表于 2023-8-5 17:35

还是不行大佬，结果栅格值还是从0-255一堆

陶远航

我已经阅读了您的代码，发现了一些可能的错误和改进的地方。以下是我的建议：

1.您的代码中，您使用了tifffile库来读取和写入TIFF格式的影像数据。这个库主要用于处理多页或多通道的TIFF文件，而您的影像数据只有三个通道（RGB）。因此，您可以使用其他库，如PIL或OpenCV，来更简单地读取和写入影像数据。
2.您的代码中，您使用了np.mean函数来将输入图像转换为灰度图像。这个函数会对每个像素的三个通道求平均值，得到一个灰度值。然而，这种方法并不是最佳的灰度化方法，因为它没有考虑到人眼对不同颜色的敏感度。一种更好的方法是使用加权平均法，根据不同颜色的亮度系数来计算灰度值。
例如，您可以使用以下公式：

gray=0.299∗R+0.587∗G+0.114∗B

3.您的代码中，您使用了一个双重循环来实现自适应阈值二值化。这个方法非常耗时，并且容易受到噪声的影响。您可以使用OpenCV库中的adaptiveThreshold函数来更快地实现自适应阈值二值化。这个函数可以根据每个像素周围的区域来计算一个局部阈值，并且可以选择不同的阈值算法和参数。
4.您的代码中，您使用了gdal库来设置输出图像的空间参考信息。这个库是一个强大的地理数据处理库，但是它也比较复杂和底层。如果您只是想要保存一个简单的TIFF文件，并不需要设置空间参考信息，您可以直接使用tifffile或PIL或OpenCV等库来保存输出图像。

综上所述，我为您修改了一下代码，如下所示：

# 导入所需的库
from PIL import Image
import numpy as np
import cv2

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

# 使用PIL库打开图像并转换为numpy数组
image = np.array(Image.open(image_path))

# 使用加权平均法将输入图像转换为灰度图像
gray_image = 0.299 * image[:,:,0] + 0.587 * image[:,:,1] + 0.114 * image[:,:,2]

# 使用OpenCV库实现自适应阈值二值化
window_size = 131 # 窗口大小
k = 0.3 # 控制阈值的参数
binary = cv2.adaptiveThreshold(gray_image.astype(np.uint8), 255, cv2.ADAPTIVE_THRESH_MEAN_C, cv2.THRESH_BINARY, window_size, k)

# 使用PIL库保存输出图像
Image.fromarray(binary).save(output_path)

print("保存成功！")

我希望这些修改能够帮助您提高无人机影像处理沙丘脊线提取的效果。

高山

Justheartyoung 发表于 2023-8-6 10:45
还是不行大佬，结果栅格值还是从0-255一堆

ai的，看风格