问题描述

最近需要对无人机影像中绿色植被信息进行提取，查看相关论文，发现用的比较多的就是HSV色彩转换方法，动手实践一下。

解决思路

解决过程

HSV转换

import os
import numpy as np
from osgeo import gdal
'''
本代码实现将RGB波段的无人机影像转换到HSV色彩空间。
注意无人机影像波段顺序为R,G,B
'''
def RGB2HSV(file):src = gdal.Open(file)col = src.RasterXSizerow = src.RasterYSize_ = 3red = np.array(src.GetRasterBand(1).ReadAsArray().astype(float))green = np.array(src.GetRasterBand(2).ReadAsArray().astype(float))blue = np.array(src.GetRasterBand(3).ReadAsArray().astype(float))arr = np.stack((red, green, blue), axis=2)dst = np.zeros((row, col, 3), dtype=np.float32)for i in range(row):for j in range(col):r, g, b = arr[i, j] / 255.0minn = np.min([r, g, b])maxx = np.max([r, g, b])dst[i, j, 2] = maxx  # Vdelta = maxx - minnh, s = 0, 0if maxx != 0:s = delta / maxxif r == maxx:h = (g - b) / deltaelif g == maxx:h = 2 + (b - r) / deltaelse:h = 4 + (r - g) / deltah *= 60if h < 0:h += 360dst[i, j, 0] = hdst[i, j, 1] = sprint('开始输出')out_name = 'hsv.tif'out_ds = gdal.GetDriverByName('GTiff').Create(out_name, col, row, 3, gdal.GDT_Float32)for i in range(3):# data = out_ds.GetRasterBand(i+1).ReadAsArray()band = out_ds.GetRasterBand(i+1).WriteArray(dst[:,:,i])del bandout_ds.SetProjection(src.GetProjection())out_ds.SetGeoTransform(src.GetGeoTransform())out_ds.FlushCache()return out_dsfile = r'G:\temp\forestVFC\标注影像\ZJM_ZLHF_2308221.dat'# Convert to HSV
dst = RGB2HSV(file)

直方图查看

原结果如下：

HSV结果如下所示：

查看HSV结果波段1的直方图
我这里的需求是分为植被和非植被，查看直方图中第一个波谷的折点即为植被和非植被区域阈值，也就是下图中的28.1115

掩膜提取

在ENVI——bandmath中使用下面这个公式对原始结果进行掩膜即可得到绿色植被区域的影像。

(float(b1) ge 28.0 ) * float(b2)
# b1 为HSV影像的波段1，b2为原始无人机影像

参考：
https://blog.csdn.net/u012294613/article/details/141096007

版权声明：本文为博主原创文章，遵循 CC 4.0 BY-SA 版权协议，转载请附上原文出处链接和本声明。