云检测-IDL实现2 IDL

一，前言
之前我的博客里面简单介绍了针对高分一高分辨率卫星遥感影像的厚云监测方法，主要是利用云的特征，基于阈值法进行提取的云，该方法简单，但是只能提取出绝大部分的厚云，对薄云提取不是很友好，后续我对该算法进行了改进，使得精度更高，具体流程图如下

文章图片

该算法主要是在原先阈值算法的基础上增加了影像处理，通过特定的系数，对RGB色彩影像进行处理，生成的影像更加突出云特征，这样可以大大提高云检测的精度
二，算法实现
【云检测-IDL实现2】主要通过IDL语言实现的，以下是具体的实现过程

PRO cloud_snow_dection,input_DEM=input_DEM,file=file COMPILE_OPT IDL2 e=envi(/headless) file='文件路径' raster=e.OpenRaster(file) ; 辐射定标为表观反射率 Task = ENVITask('RadiometricCalibration') Task.INPUT_RASTER =raster Task.CALIBRATION_TYPE='Top-of-Atmosphere Reflectance'; GF6为‘Radioance’ Task.OUTPUT_DATA_TYPE = 'Double' Task.Execute raster=Task.Output_Raster; raster为辐射定标的TOA数据 ; 潜在云像元估计 band1=raster.GetData(band=0); 蓝波段 band2=raster.GetData(band=1); 绿波段 band3=raster.GetData(band=2); 红波段 band4=raster.GetData(band=3); 近红波段 fid=ENVIRASTERTOFID(raster) envi_file_query,fid,dims=dims,ns=ns,nl=nl,nb=nb ; 统计蓝波段各个方法获取的自动阈值 iso=IMAGE_THRESHOLD(band1,threshold=thresh_value1,/ISODATA) otsu=IMAGE_THRESHOLD(band1, threshold=thresh_value2,/OTSU) moments=IMAGE_THRESHOLD(band1, threshold=thresh_value3,/MOMENTS) me=IMAGE_THRESHOLD(band1, threshold=thresh_value4,/MAXENTROPY) mean=IMAGE_THRESHOLD(band1,threshold=thresh_value5, /MEAN) iso_band3=IMAGE_THRESHOLD(band3,threshold=band3_threshold,/MAXENTROPY) a=[thresh_value1[0],thresh_value2[0],thresh_value3[0],thresh_value4[0],thresh_value5[0]] value=https://www.it610.com/article/max(a); 取最大阈值作为最后的阈值 condition1=band1 gt value condition2=band3 gt band3_threshold[0] condition=condition1 and condition2; 云象元初提取 ; 基于灰度图像方法进行云象元二次提取 ; 按照特定的方式对RGB进行重整，生成特定的灰度图像参考文献Saranya M.cloud removal from satellite images using information cloning index_data=[[0.229,0.587,0.114],[0.595716,-0.274453,-0.321263],[0.211456,-0.522591,0.311135]] R=reform(band3,ns*nl,1) G=reform(band2,ns*nl,1) B=reform(band1,ns*nl,1) temp=make_array(ns*nl,3) temp[*,0]=R temp[*,1]=G temp[*,2]=B ; A##B表示矩阵相乘，用A行乘以B的列，A#B表示用A列乘以B的行 Y=index_data[*,0]##temp I=index_data[*,1]##temp Q=index_data[*,2]##temp y_data=reform(Y,ns,nl) I_data=reform(I,ns,nl) Q_data=reform(Q,ns,nl) temp_raster=make_array(ns,nl,3) temp_raster[*,*,0]=y_data temp_raster[*,*,1]=I_data temp_raster[*,*,2]=Q_data Y_I_Q_raster=ENVIRASTER(temp_raster,uri=e.GetTemporaryFilename()) Y_I_Q_raster.Save y_data=Y_I_Q_raster.GetData(band=0) I_data=Y_I_Q_raster.GetData(band=1) Q_data=Y_I_Q_raster.GetData(band=2) index1=IMAGE_THRESHOLD(y_data, threshold=y_value,/MAXENTROPY) index3=IMAGE_THRESHOLD(Q_data, threshold=q_value,/MAXENTROPY) index2=IMAGE_THRESHOLD(I_data, threshold=i_value,/OTSU) condition3=y_data gt y_value[0] condition4=Q_data gt q_value[0] condition5=I_data gt i_value[0] condition1_1=condition3 and condition4 and condition5 cloud_mask=condition and condition1_1; 生成的云的掩膜的二值文件 END

三、结果展示