pro Part_5
compile_opt idl2
e=envi(/headless)
;建立结果文件夹
file_mkdir,"D:\Lansat_Data_Composition\"
for i=1,91,1 do begin
ii=strtrim(string(i),1)
;if (i lt 23) then continue
;if ((i lt 13) or (i gt 19)) then continue
if (i ne 17) then continue
;读取文件名、含云量、图像完整度
folder="D:\Lansat_Data_Cloud_Result\"+ii+"\"
filearr = file_search(folder,'*.txt',count=num)
arr_1_1_temp = strarr(1,num)
arr_1_2_temp = strarr(1,num)
arr_2_1_temp = fltarr(1,num)
arr_2_2_temp = fltarr(1,num)
arr_3_1_temp = fltarr(1,num)
arr_3_2_temp = fltarr(1,num)
n_1 = 0
n_2 = 0
;打开第i个文件夹的第 j号txt文件
for j=0,num-1,1 do begin
txt_file = filearr[j]
openr,1,txt_file
row = file_lines(txt_file)
box_data = strarr(row)
readf, 1, box_data
free_lun,1
;先剔除不合格数据
if (box_data[5] ne 0) and (box_data[4] lt 80) then begin
if (box_data[5] gt 95) then begin
;完整性好的赋值给arr_n_1_temp
arr_1_1_temp[n_1] = box_data[0]
arr_2_1_temp[n_1] = box_data[4]
arr_3_1_temp[n_1] = box_data[5]
;计数合格元素的数量
n_1 = n_1 + 1
;完整性差的赋值给arr_n_2_temp
endif else if((box_data[5] le 95) and (box_data[5] gt 5)) then begin
arr_1_2_temp[n_2] = box_data[0]
arr_2_2_temp[n_2] = box_data[4]
arr_3_2_temp[n_2] = box_data[5]
n_2 = n_2 + 1
endif
endif
endfor
;依据前方计数的有效元素数量n_1，删除前三个个数组之中的空元素
if (n_1 ge 1) then begin
arr_1_1 = strarr(1,n_1)
arr_2_1 = fltarr(1,n_1)
arr_3_1 = fltarr(1,n_1)
for k=0,n_1-1,1 do begin
arr_1_1[k] = arr_1_1_temp[k]
arr_2_1[k] = arr_2_1_temp[k]
arr_3_1[k] = arr_3_1_temp[k]
endfor
if (n_1 ge 2) then begin
;开始冒泡排序法，对第一组数据，按照图像质量进行排序
temp_1_1 = 'PASS'
temp_2_1 = 0.0
temp_3_1 = 0.0
for m=0,n_1-2,1 do begin
for n=0,n_1-2,1 do begin
if(arr_2_1[n] gt arr_2_1[n+1]) then begin
temp_1_1 = arr_1_1[n]
temp_2_1 = arr_2_1[n]
temp_3_1 = arr_3_1[n]
arr_1_1[n] = arr_1_1[n+1]
arr_2_1[n] = arr_2_1[n+1]
arr_3_1[n] = arr_3_1[n+1]
arr_1_1[n+1] = temp_1_1
arr_2_1[n+1] = temp_2_1
arr_3_1[n+1] = temp_3_1
endif
endfor
endfor
endif
endif
;依据前方计数的有效元素数量n_2，删除前三个个数组之中的空元素
if (n_2 ge 1) then begin
arr_1_2 = strarr(1,n_2)
arr_2_2 = fltarr(1,n_2)
arr_3_2 = fltarr(1,n_2)
for l=0,n_2-1,1 do begin
arr_1_2[l] = arr_1_2_temp[l]
arr_2_2[l] = arr_2_2_temp[l]
arr_2_2[l] = arr_3_2_temp[l]
endfor
;开始冒泡排序法，对第二组数据，按照图像质量进行排序
if (n_2 ge 1) then begin
arr_1_2 = strarr(1,n_2)
arr_2_2 = fltarr(1,n_2)
arr_3_2 = fltarr(1,n_2)
for k=0,n_2-1,1 do begin
arr_1_2[k] = arr_1_2_temp[k]
arr_2_2[k] = arr_2_2_temp[k]
arr_3_2[k] = arr_3_2_temp[k]
endfor
if (n_2 ge 2) then begin
;开始冒泡排序法，按照图像质量进行排序
temp_1_2 = 'PASS'
temp_2_2 = 0.0
temp_3_2 = 0.0
for m=0,n_2-2,1 do begin
for n=0,n_2-2,1 do begin
if(arr_2_2[n] gt arr_2_2[n+1]) then begin
temp_1_2 = arr_1_2[n]
temp_2_2 = arr_2_2[n]
temp_3_2 = arr_3_2[n]
arr_1_2[n] = arr_1_2[n+1]
arr_2_2[n] = arr_2_2[n+1]
arr_3_2[n] = arr_3_2[n+1]
arr_1_2[n+1] = temp_1_2
arr_2_2[n+1] = temp_2_2
arr_3_2[n+1] = temp_3_2
endif
endfor
endfor
endif
endif
endif
;合并两个排序之后的数组得到运算结果
if ((n_1 ge 1) and(n_2 ge 1)) then begin
arr_1 = [[arr_1_1],[arr_1_2]]
arr_2 = [[arr_2_1],[arr_2_2]]
arr_3 = [[arr_3_1],[arr_3_2]]
endif else if ((n_1 ge 1) and (n_2 lt 1)) then begin
arr_1 = arr_1_1
arr_2 = arr_2_1
arr_3 = arr_3_1
endif else if ((n_1 lt 1) and (n_2 ge 1)) then begin
arr_1 = arr_1_2
arr_2 = arr_2_2
arr_3 = arr_3_2
endif else begin
arr_1 = ['PASS']
arr_2 = [0.0]
arr_3 = [0.0]
endelse
;由于上边的图像给出的是云检测结果的文件名，我们需要改成裁切后结果的文件名
;将第一张影像作为底图，后续(n_1+n_2-1)张作为副图
name_part_1 = "D:\Lansat_Data_分区切割"
if i lt 10 then begin
name_part_2 = strmid(arr_1[0],20,56)
endif else if (i ge 10) then begin
name_part_2 = strmid(arr_1[0],20,57)
endif
pic_base = name_part_1 + name_part_2
raster = e.OpenRaster(pic_base)
fid = ENVIRasterToFID(raster)
MAP_INFO_base = ENVI_GET_MAP_INFO(fid=fid)
ENVI_FILE_QUERY, fid, ns=ns, nl=nl, nb=nb, dims=dim, data_type=data_type
base_data = raster.getdata()
;将排序好的数组循环(n_1+n_2-1)次，按照波段比值公式
for p=1,n_1+n_2-1,1 do begin
if i lt 10 then begin
name_part_3 = strmid(arr_1[p],20,56)
endif else if (i ge 10) then begin
name_part_3 = strmid(arr_1[p],20,57)
endif
pic_assistant = name_part_1 + name_part_3
raster = e.OpenRaster(pic_assistant)
fid = ENVIRasterToFID(raster)
MAP_INFO_base = ENVI_GET_MAP_INFO(fid=fid)
ENVI_FILE_QUERY, fid, ns=ns, nl=nl, nb=nb, dims=dim, data_type=data_type
assistant_data = raster.getdata()
;波段比值细化，选择更佳像素
for q=0,1516,1 do begin
for r=0,1016,1 do begin
;首先计算它的DN值总和的范围
num_all_1 = long(base_data[r,q,0])+base_data[r,q,1]+base_data[r,q,2]
num_all_2 = long(assistant_data[r,q,0])+assistant_data[r,q,1]+assistant_data[r,q,2]
;使用CASE…OF…语句比IF…ELSE…更快
;注意这里不是等于0而是不存在
if((num_all_1 + num_all_2) eq 0) then continue
if((num_all_1 ne 0) and (num_all_2 eq 0)) then continue
if ((num_all_1 eq 0) and (num_all_2 gt 0)) then begin
base_data[r,q,0] = assistant_data[r,q,0]
base_data[r,q,1] = assistant_data[r,q,1]
base_data[r,q,2] = assistant_data[r,q,2]
endif
;以下都要求num_all_2大于等于25000
if((num_all_1 ne 0) and (num_all_2 ne 0)) then begin
;0~50000,直接PASS
if ((num_all_1 gt 0) and (num_all_1 lt 50000)) then continue
;50000~60000
if ((num_all_1 ge 50000) and (num_all_1 lt 60000)) then begin
ratio_value_1 = float(base_data[r,q,0])/float(base_data[r,q,2])
ratio_value_2 = float(assistant_data[r,q,0])/float(assistant_data[r,q,2])
if ((ratio_value_1 ge 0.9) and (ratio_value_2 le 1.1)) then begin
if (((num_all_2 ge 27000) and(num_all_2 lt 50000)) or (((num_all_2 ge 50000) and (num_all_2 lt 60000)) and ((ratio_value_2 lt 0.9) or (ratio_value_2 gt 1.1))) or (((num_all_2 ge 60000) and (num_all_2 lt 70000)) and ((ratio_value_2 lt 0.8) or (ratio_value_2 gt 1.75))) or (((num_all_2 ge 70000) and (num_all_2 lt 80000)) and ((ratio_value_2 lt 0.9) or (ratio_value_2 gt 1.9))) or (((num_all_2 ge 80000) and (num_all_2 lt 90000)) and ((ratio_value_2 lt 0.9) or (ratio_value_2 gt 2.5))) or (((num_all_2 ge 90000) and (num_all_2 lt 100000)) and ((ratio_value_2 lt 0.95) or (ratio_value_2 gt 3.5))) and (((num_all_2 ge 100000) or (num_all_2 lt 200000)) and ((ratio_value_2 lt 0.95) or (ratio_value_2 gt 4.0)))) then begin
base_data[r,q,0] = assistant_data[r,q,0]
base_data[r,q,1] = assistant_data[r,q,1]
base_data[r,q,2] = assistant_data[r,q,2]
endif
endif
endif
;60000~70000
if ((num_all_1 ge 60000) and (num_all_1 lt 70000)) then begin
ratio_value_1 = float(base_data[r,q,0])/float(base_data[r,q,2])
ratio_value_2 = float(assistant_data[r,q,0])/float(assistant_data[r,q,2])
if (ratio_value_1 ge 0.8) and (ratio_value_2 le 1.75) then begin
if (((num_all_2 ge 27000) and(num_all_2 lt 50000)) or (((num_all_2 ge 50000) and (num_all_2 lt 60000)) and ((ratio_value_2 lt 0.9) or (ratio_value_2 gt 1.1))) or (((num_all_2 ge 60000) and (num_all_2 lt 70000)) and ((ratio_value_2 lt 0.8) or (ratio_value_2 gt 1.75))) or (((num_all_2 ge 70000) and (num_all_2 lt 80000)) and ((ratio_value_2 lt 0.9) or (ratio_value_2 gt 1.9))) or (((num_all_2 ge 80000) and (num_all_2 lt 90000)) and ((ratio_value_2 lt 0.9) or (ratio_value_2 gt 2.5))) or (((num_all_2 ge 90000) and (num_all_2 lt 100000)) and ((ratio_value_2 lt 0.95) or (ratio_value_2 gt 3.5))) and (((num_all_2 ge 100000) or (num_all_2 lt 200000)) and ((ratio_value_2 lt 0.95) or (ratio_value_2 gt 4.0)))) then begin
base_data[r,q,0] = assistant_data[r,q,0]
base_data[r,q,1] = assistant_data[r,q,1]
base_data[r,q,2] = assistant_data[r,q,2]
endif
endif
endif
;70000~80000
if ((num_all_1 ge 70000) and (num_all_1 lt 80000)) then begin
ratio_value_1 = float(base_data[r,q,0])/float(base_data[r,q,2])
ratio_value_2 = float(assistant_data[r,q,0])/float(assistant_data[r,q,2])
if (ratio_value_1 ge 0.9) and (ratio_value_2 le 1.9) then begin
if (((num_all_2 ge 27000) and(num_all_2 lt 50000)) or (((num_all_2 ge 50000) and (num_all_2 lt 60000)) and ((ratio_value_2 lt 0.9) or (ratio_value_2 gt 1.1))) or (((num_all_2 ge 60000) and (num_all_2 lt 70000)) and ((ratio_value_2 lt 0.8) or (ratio_value_2 gt 1.75))) or (((num_all_2 ge 70000) and (num_all_2 lt 80000)) and ((ratio_value_2 lt 0.9) or (ratio_value_2 gt 1.9))) or (((num_all_2 ge 80000) and (num_all_2 lt 90000)) and ((ratio_value_2 lt 0.9) or (ratio_value_2 gt 2.5))) or (((num_all_2 ge 90000) and (num_all_2 lt 100000)) and ((ratio_value_2 lt 0.95) or (ratio_value_2 gt 3.5))) and (((num_all_2 ge 100000) or (num_all_2 lt 200000)) and ((ratio_value_2 lt 0.95) or (ratio_value_2 gt 4.0)))) then begin
base_data[r,q,0] = assistant_data[r,q,0]
base_data[r,q,1] = assistant_data[r,q,1]
base_data[r,q,2] = assistant_data[r,q,2]
endif
endif
endif
;80000~90000
if ((num_all_1 ge 80000) and (num_all_1 lt 90000)) then begin
ratio_value_1 = float(base_data[r,q,0])/float(base_data[r,q,2])
ratio_value_2 = float(assistant_data[r,q,0])/float(assistant_data[r,q,2])
if (ratio_value_1 ge 0.9) and (ratio_value_2 le 2.5) then begin
if (((num_all_2 ge 27000) and(num_all_2 lt 50000)) or (((num_all_2 ge 50000) and (num_all_2 lt 60000)) and ((ratio_value_2 lt 0.9) or (ratio_value_2 gt 1.1))) or (((num_all_2 ge 60000) and (num_all_2 lt 70000)) and ((ratio_value_2 lt 0.8) or (ratio_value_2 gt 1.75))) or (((num_all_2 ge 70000) and (num_all_2 lt 80000)) and ((ratio_value_2 lt 0.9) or (ratio_value_2 gt 1.9))) or (((num_all_2 ge 80000) and (num_all_2 lt 90000)) and ((ratio_value_2 lt 0.9) or (ratio_value_2 gt 2.5))) or (((num_all_2 ge 90000) and (num_all_2 lt 100000)) and ((ratio_value_2 lt 0.95) or (ratio_value_2 gt 3.5))) and (((num_all_2 ge 100000) or (num_all_2 lt 200000)) and ((ratio_value_2 lt 0.95) or (ratio_value_2 gt 4.0)))) then begin
base_data[r,q,0] = assistant_data[r,q,0]
base_data[r,q,1] = assistant_data[r,q,1]
base_data[r,q,2] = assistant_data[r,q,2]
endif
endif
endif
;90000~100000
if ((num_all_1 ge 90000) and (num_all_1 lt 100000)) then begin
ratio_value_1 = float(base_data[r,q,0])/float(base_data[r,q,2])
ratio_value_2 = float(assistant_data[r,q,0])/float(assistant_data[r,q,2])
if (ratio_value_1 ge 0.95) and (ratio_value_2 le 3.5) then begin
if (((num_all_2 ge 27000) and(num_all_2 lt 50000)) or (((num_all_2 ge 50000) and (num_all_2 lt 60000)) and ((ratio_value_2 lt 0.9) or (ratio_value_2 gt 1.1))) or (((num_all_2 ge 60000) and (num_all_2 lt 70000)) and ((ratio_value_2 lt 0.8) or (ratio_value_2 gt 1.75))) or (((num_all_2 ge 70000) and (num_all_2 lt 80000)) and ((ratio_value_2 lt 0.9) or (ratio_value_2 gt 1.9))) or (((num_all_2 ge 80000) and (num_all_2 lt 90000)) and ((ratio_value_2 lt 0.9) or (ratio_value_2 gt 2.5))) or (((num_all_2 ge 90000) and (num_all_2 lt 100000)) and ((ratio_value_2 lt 0.95) or (ratio_value_2 gt 3.5))) and (((num_all_2 ge 100000) or (num_all_2 lt 200000)) and ((ratio_value_2 lt 0.95) or (ratio_value_2 gt 4.0)))) then begin
base_data[r,q,0] = assistant_data[r,q,0]
base_data[r,q,1] = assistant_data[r,q,1]
base_data[r,q,2] = assistant_data[r,q,2]
endif
endif
endif
;100000~150000
if ((num_all_1 gt 100000) and (num_all_2 lt 200000)) then begin
ratio_value_1 = float(base_data[r,q,0])/float(base_data[r,q,2])
ratio_value_2 = float(assistant_data[r,q,0])/float(assistant_data[r,q,2])
if (ratio_value_1 ge 0.95) and (ratio_value_2 le 4.0) then begin
if (((num_all_2 ge 27000) and(num_all_2 lt 50000)) or (((num_all_2 ge 50000) and (num_all_2 lt 60000)) and ((ratio_value_2 lt 0.9) or (ratio_value_2 gt 1.1))) or (((num_all_2 ge 60000) and (num_all_2 lt 70000)) and ((ratio_value_2 lt 0.8) or (ratio_value_2 gt 1.75))) or (((num_all_2 ge 70000) and (num_all_2 lt 80000)) and ((ratio_value_2 lt 0.9) or (ratio_value_2 gt 1.9))) or (((num_all_2 ge 80000) and (num_all_2 lt 90000)) and ((ratio_value_2 lt 0.9) or (ratio_value_2 gt 2.5))) or (((num_all_2 ge 90000) and (num_all_2 lt 100000)) and ((ratio_value_2 lt 0.95) or (ratio_value_2 gt 3.5))) and (((num_all_2 ge 100000) or (num_all_2 lt 200000)) and ((ratio_value_2 lt 0.95) or (ratio_value_2 gt 4.0)))) then begin
base_data[r,q,0] = assistant_data[r,q,0]
base_data[r,q,1] = assistant_data[r,q,1]
base_data[r,q,2] = assistant_data[r,q,2]
endif
endif
endif
endif
endfor
endfor
endfor
;保存图像
out_path = "D:\Lansat_Data_Composition\"+ii+".dat"
ENVI_WRITE_ENVI_FILE, base_data, out_name=out_path, /nocopy, ns=ns, nl=nl, nb=nb, offset=offset, bnames=bnames, map_info=map_info_base
print,"OK",ii
endfor
end
1135
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
