|
|
发表于 2023-8-5 15:13:46
|
显示全部楼层
输出结果的栅格值奇怪的原因是代码中没有对二值化后的图像进行数据类型转换。在创建新的输出图像时,应该将数据类型设置为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("保存成功!")
复制代码
这样修改后,输出结果的栅格值应该就不会出现奇怪的情况了。 |
|