请大哥们帮弟弟看看，问题出在哪？并麻烦给出改正代码，谢谢！

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 = 11 # 窗口大小
k = 0.1 # 控制阈值的参数
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)

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

print("保存成功！")

在图片预览加载和在Arcgis软件中加载情况如下：

isdkz

问题出在代码最后将二值化结果写入输出图像的部分。在使用 out_band.WriteArray(binary) 将数据写入输出图像时，需要指定数据的数据类型。在这里，二值化结果 binary 的数据类型是 np.uint8 ，但是在写入输出图像时，指定的数据类型是 gdal.GDT_Byte 。这两者并不一致，导致输出图像无法正确显示。

要解决这个问题，需要将二值化结果的数据类型转换为与输出图像的数据类型一致。修改代码如下：

1. # 将沙脊线部分写入输出图像
   
2. out_band = out_ds.GetRasterBand(1)
   
3. out_band.WriteArray(binary.astype(np.uint8))
   
4. 
   
5. # 释放资源
   
6. out_band.FlushCache()
   
7. out_ds = None
   
8. in_ds = None
   

复制代码

这样就可以将二值化结果正确地写入输出图像了。重新运行代码，输出的图像就可以在Arcgis软件中正确显示了。

学习编程中的Ben

在你提供的代码中，有几个可能导致输出结果在ArcGIS软件中无法显示的问题。以下是需要修改的部分：

1. 输出数据类型问题：你使用了gdal.GDT_Byte来指定输出图像的数据类型，但这可能会限制图像的灰度范围。根据你的代码逻辑，你应该使用gdal.GDT_UInt16作为数据类型，以便能够保持二值化结果的灰度范围。

1. out_ds = driver.Create(output_path, image.shape[1], image.shape[0], 1, gdal.GDT_UInt16)

复制代码

2. 设置颜色表：对于二值图像，在ArcGIS中显示时，你可能需要设置颜色表以正确显示沙脊线。你可以使用SetColorTable函数来定义一个颜色表，然后将其应用于输出图像的波段。

1. color_table = gdal.ColorTable()
2. color_table.SetColorEntry(0, (0, 0, 0, 0))  # 黑色，对应值为0
3. color_table.SetColorEntry(255, (255, 255, 255, 255))  # 白色，对应值为255
4. out_band.SetColorTable(color_table)

复制代码

3. 写入数据范围问题：在写入二值图像数据之前，你需要确保将数据缩放到正确的范围（0-255）。否则，ArcGIS可能无法正确解释灰度值。

1. out_band.WriteArray(binary * 255)  # 将二值图像数据缩放到0-255的范围

复制代码

请按照上述修改，重新运行代码并尝试在ArcGIS中加载结果。希望这能解决你的问题！如果还有其他问题，请随时提问。

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

Justheartyoung

isdkz 发表于 2023-8-5 10:44
问题出在代码最后将二值化结果写入输出图像的部分。在使用 out_band.WriteArray(binary) 将数据写入输出图 ...

运行了一下您的代码，还是不行，还是空白

Justheartyoung

学习编程中的Ben 发表于 2023-8-5 10:44
在你提供的代码中，有几个可能导致输出结果在ArcGIS软件中无法显示的问题。以下是需要修改的部分：

1. 输 ...

可以帮我和上述原代码整理合并一下吗？谢谢

Mike_python小

Justheartyoung 发表于 2023-8-5 08:14
可以帮我和上述原代码整理合并一下吗？谢谢

1. from tifffile import imread
   
2. import numpy as np
   
3. from osgeo import gdal, osr
   
4. 
   
5. # 读取TIFF格式无人机影像数据
   
6. image_path = 'F:\\duneline\\dune\\dune.tif'
   
7. output_path = 'overlay_image.tif'
   
8. 
   
9. image = imread(image_path)
   
10. 
    
11. # 将输入图像转换为灰度图像
    
12. gray_image = np.mean(image, axis=2)
    
13. 
    
14. # 自适应阈值二值化
    
15. binary = np.zeros_like(gray_image, dtype=np.uint8)
    
16. window_size = 11 # 窗口大小
    
17. k = 0.1 # 控制阈值的参数
    
18. for i in range(window_size//2, gray_image.shape[0]-window_size//2):
    
19.     for j in range(window_size//2, gray_image.shape[1]-window_size//2):
    
20.         window = gray_image[i-window_size//2:i+window_size//2+1, j-window_size//2:j+window_size//2+1]
    
21.         threshold = np.mean(window) - k * np.std(window)
    
22.         if gray_image[i, j] > threshold:
    
23.             binary[i, j] = 255
    
24. 
    
25. # 配置输出的空间参考信息
    
26. in_ds = gdal.Open(image_path)
    
27. projection = in_ds.GetProjection()
    
28. geotransform = in_ds.GetGeoTransform()
    
29. 
    
30. # 创建新的输出图像
    
31. driver = gdal.GetDriverByName('GTiff')
    
32. out_ds = driver.Create(output_path, image.shape[1], image.shape[0], 1, gdal.GDT_UInt16)
    
33. 
    
34. # 设置输出图像的空间参考信息
    
35. out_ds.SetProjection(projection)
    
36. out_ds.SetGeoTransform(geotransform)
    
37. 
    
38. # 设置颜色表
    
39. color_table = gdal.ColorTable()
    
40. color_table.SetColorEntry(0, (0, 0, 0, 0))  # 黑色，对应值为0
    
41. color_table.SetColorEntry(255, (255, 255, 255, 255))  # 白色，对应值为255
    
42. out_band.SetColorTable(color_table)
    
43. 
    
44. # 将沙脊线部分写入输出图像
    
45. out_band = out_ds.GetRasterBand(1)
    
46. out_band.WriteArray(binary * 255)  # 将二值图像数据缩放到0-255的范围
    
47. 
    
48. # 释放资源
    
49. out_band.FlushCache()
    
50. out_ds = None
    
51. in_ds = None
    
52. 
    
53. print("保存成功！")
    

复制代码

如果回答对你有帮助，请设置最佳答案！

Justheartyoung

Mike_python小 发表于 2023-8-5 12:08
如果回答对你有帮助，请设置最佳答案！

大佬，报错了：
Traceback (most recent call last):
  File "F:/duneline/duneline.py", line 145, in <module>
    out_band.SetColorTable(color_table)
NameError: name 'out_band' is not defined

Mike_python小

Justheartyoung 发表于 2023-8-5 09:32
大佬，报错了：
Traceback (most recent call last):
  File "F:/duneline/duneline.py", line 145, in ...

1. from tifffile import imread
   
2. import numpy as np
   
3. from osgeo import gdal, osr
   
4. 
   
5. # 读取TIFF格式无人机影像数据
   
6. image_path = 'F:\\duneline\\dune\\dune.tif'
   
7. output_path = 'overlay_image.tif'
   
8. 
   
9. image = imread(image_path)
   
10. 
    
11. # 将输入图像转换为灰度图像
    
12. gray_image = np.mean(image, axis=2)
    
13. 
    
14. # 自适应阈值二值化
    
15. binary = np.zeros_like(gray_image, dtype=np.uint8)
    
16. window_size = 11 # 窗口大小
    
17. k = 0.1 # 控制阈值的参数
    
18. for i in range(window_size//2, gray_image.shape[0]-window_size//2):
    
19.     for j in range(window_size//2, gray_image.shape[1]-window_size//2):
    
20.         window = gray_image[i-window_size//2:i+window_size//2+1, j-window_size//2:j+window_size//2+1]
    
21.         threshold = np.mean(window) - k * np.std(window)
    
22.         if gray_image[i, j] > threshold:
    
23.             binary[i, j] = 255
    
24. 
    
25. # 配置输出的空间参考信息
    
26. in_ds = gdal.Open(image_path)
    
27. projection = in_ds.GetProjection()
    
28. geotransform = in_ds.GetGeoTransform()
    
29. 
    
30. # 创建新的输出图像
    
31. driver = gdal.GetDriverByName('GTiff')
    
32. out_ds = driver.Create(output_path, image.shape[1], image.shape[0], 1, gdal.GDT_UInt16)
    
33. 
    
34. # 设置输出图像的空间参考信息
    
35. out_ds.SetProjection(projection)
    
36. out_ds.SetGeoTransform(geotransform)
    
37. 
    
38. # 将沙脊线部分写入输出图像
    
39. out_band = out_ds.GetRasterBand(1)
    
40. out_band.WriteArray(binary * 255)  # 将二值图像数据缩放到0-255的范围
    
41. 
    
42. # 设置颜色表
    
43. color_table = gdal.ColorTable()
    
44. color_table.SetColorEntry(0, (0, 0, 0, 0))  # 黑色，对应值为0
    
45. color_table.SetColorEntry(255, (255, 255, 255, 255))  # 白色，对应值为255
    
46. out_band.SetColorTable(color_table)
    
47. 
    
48. # 释放资源
    
49. out_band.FlushCache()
    
50. out_ds = None
    
51. in_ds = None
    
52. 
    
53. print("保存成功！")
    

复制代码