鱼C论坛

 找回密码
 立即注册
查看: 1241|回复: 4

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

[复制链接]
发表于 2023-8-5 13:19:12 | 显示全部楼层 |阅读模式

马上注册,结交更多好友,享用更多功能^_^

您需要 登录 才可以下载或查看,没有账号?立即注册

x
下面代码是想对无人机影像(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("保存成功!")
结果:

最佳答案
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("保存成功!")


希望这个修改能解决你的问题。
空白.jpg
想知道小甲鱼最近在做啥?请访问 -> ilovefishc.com
回复

使用道具 举报

发表于 2023-8-5 13:19:24 | 显示全部楼层
问题出在将二值图像保存到输出文件时,没有正确地进行数据类型转换和范围缩放。

在这里,将二值图像乘以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_pathimage_path变量需要替换为您自己的路径。

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

如果回答对你有帮助,请给我一个最佳答案!
想知道小甲鱼最近在做啥?请访问 -> ilovefishc.com
回复 支持 反对

使用道具 举报

发表于 2023-8-5 13:19:37 | 显示全部楼层
这个问题可能是由于保存结果时出了问题导致的。你可以尝试修改以下代码:

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("保存成功!")

希望这个修改能够解决你的问题!如果还有其他问题,请随时提问。
想知道小甲鱼最近在做啥?请访问 -> ilovefishc.com
回复 支持 反对

使用道具 举报

 楼主| 发表于 2023-8-5 13:30:27 | 显示全部楼层
isdkz 发表于 2023-8-5 13:19
这个问题可能是由于保存结果时出了问题导致的。你可以尝试修改以下代码:

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

不行,结果还是那样,全黑的
想知道小甲鱼最近在做啥?请访问 -> ilovefishc.com
回复 支持 反对

使用道具 举报

发表于 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("保存成功!")


希望这个修改能解决你的问题。
想知道小甲鱼最近在做啥?请访问 -> ilovefishc.com
回复 支持 反对

使用道具 举报

您需要登录后才可以回帖 登录 | 立即注册

本版积分规则

小黑屋|手机版|Archiver|鱼C工作室 ( 粤ICP备18085999号-1 | 粤公网安备 44051102000585号)

GMT+8, 2024-12-25 21:24

Powered by Discuz! X3.4

© 2001-2023 Discuz! Team.

快速回复 返回顶部 返回列表