| 08 Mar 2013 01:42 PM |
|
Greetings, I wrote some code to create a cloud-free mask via the model in Oreopoulos et al. (2011). I got the basic function, CalcLandsat57cloudmask, to work find from command line, but am unable to add the function anywhere on the ENVI menu (4.x and ENVI Classic). Basically it involves three PRO files:
1. CalcLandsat57cloudmask_buttons_1 - adds the button to the menu and calls the second PRO. The code currently puts it in the "Filter" menu, even though I'd rather put it elsewhere, but I know that said code worked for a different function;
2. CalcLandsat57cloudmask_widget - gets input and output file names, and sends them to the working PRO file, CalcLandsat57cloudmask;
3. CalcLandsat57cloudmask - creates a cloud-free mask, where cloud = 0, cloud-free = 1.
Does anyone know how to fix this?
Here's the code:
;------------------------------------------------------------------------------
pro CalcLandsat57cloudmask_buttons_1, buttonInfo
compile_opt idl2
;Create button so that the program can be called from within ENVI.
;Button will appear under the Filter menu
envi_define_menu_button, buttonInfo, value = 'Cloud Mask Landsat', $
event_pro='CalcLandsat57cloudmask_widget', ref_uvalue = 'convolution tool', $
ref_index=0, /sibling, uvalue='Cloud Mask Landsat'
end
;------------------------------------------------------------------------------
pro CalcLandsat57cloudmask_widget, ev
compile_opt idl2
infilename=envi_pickfile(title='Input Landsat Top-of-Atmosphere Reflectance File', filter = '*_toaref.dat')
if (infilename eq '') then return
base = widget_auto_base(title='Select output parameters and file name')
wo = widget_outf(base, uvalue='outf', /auto)
result = auto_wid_mng(base)
if (result.accept eq 0) then return
outfilename=result.outf
CalcLandsat57cloudmask, infilename, outfilename
end
;------------------------------------------------------------------------------
pro CalcLandsat57cloudmask, infilename, outfilename
compile_opt idl2, hidden
; This file will read in a top-of-atmosphere (TOA) reflectance image and produce a cloud-free mask image.
; Please note that this is NOT the same as surface reflectance.
; This routine is based upon the algorithm from Oreopoulos et al. 2011 ( http://dx.doi.org/10.1109/LGRS.2010.2095409)
; Variables:
; infilename - name of file including full path of input TOA reflectance file.
; outfilename - name of file including full path of output cloudmask file.
; This code was shamelessly adapted from an example in the ENVI help file.
; Open the datasets and query the file IDs
ENVI_OPEN_FILE, infilename, r_fid=r_fid
; Set up the bands to process
ENVI_FILE_QUERY, r_fid, dims = dims, $
nb = nb, ns = ns, nl = nl, sname = sname, data_type=data_type, xstart=xstart, ystart=ystart
out_proj=envi_get_projection(fid=r_fid,pixel_size=out_ps)
inherit = envi_set_inheritance(r_fid, dims, /spatial)
; Now create cloud mask using algorithm from Oreopoulos et al. 2011 in GRSL. This masks both zero-value pixels and clouds.
bname = 'Cloud and zero-pixel mask ('+file_basename(outfilename)+')'
ostr='Output file: ' + outfilename
rstr=["Input file: " + infilename, ostr]
envi_report_init, rstr, title='Making zero-pixel and cloud mask.', base=base, /interrupt
openw, unit, outfilename, /get_lun
tile_id1 = envi_init_tile(r_fid, 0, num_tiles=num_tiles, $
interleave=0, xs=dims[1], xe=dims[2], $
ys=dims[3], ye=dims[4])
tile_id3 = envi_init_tile(r_fid, 2, match_id=tile_id1)
tile_id4 = envi_init_tile(r_fid, 3, match_id=tile_id1)
tile_id5 = envi_init_tile(r_fid, 4, match_id=tile_id1)
tile_id6 = envi_init_tile(r_fid, 5, match_id=tile_id1)
print, "Number of samples: "+string(ns)
print, "Number of lines: "+string(nl)
print, "Number of tiles: "+string(num_tiles)
envi_report_inc, base, num_tiles
for h=0L, num_tiles-1 do begin
envi_report_stat, base, h, num_tiles, cancel=cancel
if cancel eq 1 then begin
envi_report_init, base=base, /finish
return
endif
b1 = envi_get_tile(tile_id1, h, band_index=0, ye=ye, ys=ys)
tilesize=size(b1)
b3 = envi_get_tile(tile_id3, h, band_index=2, ye=ye, ys=ys)
b4 = envi_get_tile(tile_id4, h, band_index=3, ye=ye, ys=ys)
b5 = envi_get_tile(tile_id5, h, band_index=4, ye=ye, ys=ys)
b6 = envi_get_tile(tile_id6, h, band_index=5, ye=ye, ys=ys)
b7 = make_array(ns,tilesize[2], /byte)
tilesizeb7=size(b7)
if h eq 0L then begin
print, "Input tile dimensions: " + string(tilesize[0])+" x "+string(tilesize[1])
print, tilesize
print, "output tile dimensions: " + string(tilesizeb7[0])+" x "+string(tilesizeb7[1])
print, tilesizeb7
endif
for i=0L, ns-1 do begin
for j=0L, tilesize[2]-1 do begin
if ((b1[i,j] le 0.) or (b1[i,j] gt 1.)) or((b3[i,j] le 0.) or (b3[i,j] ge 1.)) or ((b4[i,j] le 0.) or (b4[i,j] ge 1.)) or ((b5[i,j] le 0.) or (b5[i,j] ge 1.)) or ((b6[i,j] le 0.) or (b6[i,j] ge 1.)) then begin ; see if pixels match zero mask
b7[i,j]=0
endif else begin ; Non-vegetated land test
if (b1[i,j] LT b3[i,j] AND b3[i,j] LT b4[i,j] AND b4[i,j] lt (b5[i,j]*1.07) AND b5[i,j] lt 0.65) then begin
b7[i,j]=1
endif else begin ; Snow/ ice test
if (b3[i,j] GT 0.24 AND b5[i,j] LT 0.16 AND b3[i,j] GT b4[i,j]) then begin
b7[i,j]=1
endif else begin ; Water bodies test
if (b3[i,j] GT b4[i,j] AND b3[i,j] GT (b5[i,j] * 0.67) AND b1[i,j] lt 0.30 and b3[i,j] LT 0.20) then begin
b7[i,j]=1
endif else begin ; Cloud test
if ((b1[i,j] gt 0.15 OR b3[i,j] gt 0.18) AND b5[i,j] gt 0.12 AND max([b1[i,j], b3[i,j]]) gt (b5[i,j] * 0.67)) then b7[i,j]=0 else b7[i,j]=1
endelse
endelse
endelse
endelse
endfor
endfor
writeu, unit, b7
endfor
free_lun, unit
envi_setup_head, fname=outfilename, ns=ns, nl=nl, nb=1, inherit=inherit, interleave=0,$
data_type=1, offset=0, bnames=bname, descrip=bname, xstart=xstart, ystart=ystart, /write, /open
envi_tile_done, tile_id1
envi_tile_done, tile_id3
envi_tile_done, tile_id4
envi_tile_done, tile_id5
envi_tile_done, tile_id6
envi_report_init, base=base, /finish
end |
|
|
|
MariM
Veteran Member
Posts:2396  
| 11 Mar 2013 04:30 PM |
|
I didn't see any reason for making your file input/output separate from the function so I incorporated the function into the main .pro where you select the files and it will work. The button definition was not the correct syntax so I edited that. It runs now, creates a button, and produces the expected dialogs. Have a look at the below and let me know if this helps.
;==================================
pro CalcLandsat57cloudmask_define_buttons, buttonInfo
compile_opt idl2
;Create button so that the program can be called from within ENVI.
;Button will appear under the Filter menu
envi_define_menu_button, buttonInfo, value = 'Cloud Mask Landsat', $
event_pro='CalcLandsat57cloudmask', ref_uvalue = 'convolution tool', $
ref_index=0, /sibling, uvalue='Cloud Mask Landsat'
end
;------------------------------------------------------------------------------
pro CalcLandsat57cloudmask, ev
compile_opt idl2
infilename=envi_pickfile(title='Input Landsat Top-of-Atmosphere Reflectance File', filter = '*.img')
if (infilename eq '') then return
print, infilename
base = widget_auto_base(title='Select output parameters and file name')
wo = widget_outf(base, uvalue='outf', /auto)
result = auto_wid_mng(base)
if (result.accept eq 0) then return
outfilename=result.outf
print, outfilename
;CalcLandsat57cloudmask, infilename, outfilename
;pro CalcLandsat57cloudmask, infilename, outfilename
;compile_opt idl2, hidden
; This file will read in a top-of-atmosphere (TOA) reflectance image and produce a cloud-free mask image.
; Please note that this is NOT the same as surface reflectance.
; This routine is based upon the algorithm from Oreopoulos et al. 2011 ( http://dx.doi.org/10.1109/LGRS.2010.2095409)
; Variables:
; infilename - name of file including full path of input TOA reflectance file.
; outfilename - name of file including full path of output cloudmask file.
; This code was shamelessly adapted from an example in the ENVI help file.
; Open the datasets and query the file IDs
ENVI_OPEN_FILE, infilename, r_fid=r_fid
; Set up the bands to process
ENVI_FILE_QUERY, r_fid, dims = dims, $
nb = nb, ns = ns, nl = nl, sname = sname, data_type=data_type, xstart=xstart, ystart=ystart
out_proj=envi_get_projection(fid=r_fid,pixel_size=out_ps)
inherit = envi_set_inheritance(r_fid, dims, /spatial)
; Now create cloud mask using algorithm from Oreopoulos et al. 2011 in GRSL. This masks both zero-value pixels and clouds.
bname = 'Cloud and zero-pixel mask ('+file_basename(outfilename)+')'
ostr='Output file: ' + outfilename
rstr=["Input file: " + infilename, ostr]
envi_report_init, rstr, title='Making zero-pixel and cloud mask.', base=base, /interrupt
openw, unit, outfilename, /get_lun
tile_id1 = envi_init_tile(r_fid, 0, num_tiles=num_tiles, $
interleave=0, xs=dims[1], xe=dims[2], $
ys=dims[3], ye=dims[4])
tile_id3 = envi_init_tile(r_fid, 2, match_id=tile_id1)
tile_id4 = envi_init_tile(r_fid, 3, match_id=tile_id1)
tile_id5 = envi_init_tile(r_fid, 4, match_id=tile_id1)
tile_id6 = envi_init_tile(r_fid, 5, match_id=tile_id1)
print, "Number of samples: "+string(ns)
print, "Number of lines: "+string(nl)
print, "Number of tiles: "+string(num_tiles)
envi_report_inc, base, num_tiles
for h=0L, num_tiles-1 do begin
envi_report_stat, base, h, num_tiles, cancel=cancel
if cancel eq 1 then begin
envi_report_init, base=base, /finish
return
endif
b1 = envi_get_tile(tile_id1, h, band_index=0, ye=ye, ys=ys)
tilesize=size(b1)
b3 = envi_get_tile(tile_id3, h, band_index=2, ye=ye, ys=ys)
b4 = envi_get_tile(tile_id4, h, band_index=3, ye=ye, ys=ys)
b5 = envi_get_tile(tile_id5, h, band_index=4, ye=ye, ys=ys)
b6 = envi_get_tile(tile_id6, h, band_index=5, ye=ye, ys=ys)
b7 = make_array(ns,tilesize[2], /byte)
tilesizeb7=size(b7)
if h eq 0L then begin
print, "Input tile dimensions: " + string(tilesize[0])+" x "+string(tilesize[1])
print, tilesize
print, "output tile dimensions: " + string(tilesizeb7[0])+" x "+string(tilesizeb7[1])
print, tilesizeb7
endif
for i=0L, ns-1 do begin
for j=0L, tilesize[2]-1 do begin
if ((b1[i,j] le 0.) or (b1[i,j] gt 1.)) or((b3[i,j] le 0.) or (b3[i,j] ge 1.)) or ((b4[i,j] le 0.) or (b4[i,j] ge 1.)) or ((b5[i,j] le 0.) or (b5[i,j] ge 1.)) or ((b6[i,j] le 0.) or (b6[i,j] ge 1.)) then begin ; see if pixels match zero mask
b7[i,j]=0
endif else begin ; Non-vegetated land test
if (b1[i,j] LT b3[i,j] AND b3[i,j] LT b4[i,j] AND b4[i,j] lt (b5[i,j]*1.07) AND b5[i,j] lt 0.65) then begin
b7[i,j]=1
endif else begin ; Snow/ ice test
if (b3[i,j] GT 0.24 AND b5[i,j] LT 0.16 AND b3[i,j] GT b4[i,j]) then begin
b7[i,j]=1
endif else begin ; Water bodies test
if (b3[i,j] GT b4[i,j] AND b3[i,j] GT (b5[i,j] * 0.67) AND b1[i,j] lt 0.30 and b3[i,j] LT 0.20) then begin
b7[i,j]=1
endif else begin ; Cloud test
if ((b1[i,j] gt 0.15 OR b3[i,j] gt 0.18) AND b5[i,j] gt 0.12 AND max([b1[i,j], b3[i,j]]) gt (b5[i,j] * 0.67)) then b7[i,j]=0 else b7[i,j]=1
endelse
endelse
endelse
endelse
endfor
endfor
writeu, unit, b7
endfor
free_lun, unit
envi_setup_head, fname=outfilename, ns=ns, nl=nl, nb=1, inherit=inherit, interleave=0,$
data_type=1, offset=0, bnames=bname, descrip=bname, xstart=xstart, ystart=ystart, /write, /open
envi_tile_done, tile_id1
envi_tile_done, tile_id3
envi_tile_done, tile_id4
envi_tile_done, tile_id5
envi_tile_done, tile_id6
envi_report_init, base=base, /finish
end |
|
|
|
Deleted User
New Member
Posts: 
| 12 Mar 2013 03:44 AM |
|
|
[QUOTE]
for i=0L, ns-1 do begin
for j=0L, tilesize[2]-1 do begin
if ((b1[i,j] le 0.) or (b1[i,j] gt 1.)) or((b3[i,j] le 0.) or (b3[i,j] ge 1.)) or ((b4[i,j] le 0.) or (b4[i,j] ge 1.)) or ((b5[i,j] le 0.) or (b5[i,j] ge 1.)) or ((b6[i,j] le 0.) or (b6[i,j] ge 1.)) then begin ; see if pixels match zero mask
b7[i,j]=0
endif else begin ; Non-vegetated land test
if (b1[i,j] LT b3[i,j] AND b3[i,j] LT b4[i,j] AND b4[i,j] lt (b5[i,j]*1.07) AND b5[i,j] lt 0.65) then begin
b7[i,j]=1
endif else begin ; Snow/ ice test
if (b3[i,j] GT 0.24 AND b5[i,j] LT 0.16 AND b3[i,j] GT b4[i,j]) then begin
b7[i,j]=1
endif else begin ; Water bodies test
if (b3[i,j] GT b4[i,j] AND b3[i,j] GT (b5[i,j] * 0.67) AND b1[i,j] lt 0.30 and b3[i,j] LT 0.20) then begin
b7[i,j]=1
endif else begin ; Cloud test
if ((b1[i,j] gt 0.15 OR b3[i,j] gt 0.18) AND b5[i,j] gt 0.12 AND max([b1[i,j], b3[i,j]]) gt (b5[i,j] * 0.67)) then b7[i,j]=0 else b7[i,j]=1
endelse
endelse
endelse
endelse
endfor
[/QUOTE]
It is a very slow solution for IDL. You can significantly increase the speed of processing.
[QUOTE]
data_b1 = envi_get_tile(tile_id_b1, h, band_index=band_index)
data_b3 = envi_get_tile(tile_id_b3, h, band_index=band_index)
data_b4 = envi_get_tile(tile_id_b4, h, band_index=band_index)
data_b5 = envi_get_tile(tile_id_b5, h, band_index=band_index)
data_b6 = envi_get_tile(tile_id_b6, h, band_index=band_index)
data = byte(envi_get_tile(tile_id_b1, h, band_index=band_index))
data[*] = 5b
;Non-vegetated land
mask_ind=where((data_b1 lt data_b3 and $
data_b3 lt data_b4 and $
data_b4 lt 1.07*data_b5 and $
data_b5 lt 0.65), count)
if mask_ind ne [-1] then data[mask_ind]=1
;Snow/ice
mask_ind=where((data eq 5 and $
data_b3 gt 0.24 and $
data_b5 lt 0.16 and $
data_b3 gt data_b4), count)
if mask_ind ne [-1] then data[mask_ind]=1
;Water
mask_ind=where((data eq 5 and $
data_b3 gt data_b4 and $
data_b3 gt 0.67*data_b5 and $
data_b1 lt 0.30 and $
data_b3 lt 0.20), count)
if mask_ind ne [-1] then data[mask_ind]=1
;Cloud
mask_ind=where((data eq 5 and $
(data_b1 gt 0.15 or data_b3 gt 0.18) and $
data_b5 gt 0.12 and $
(data_b1 gt 0.67*data_b5 or data_b3 gt 0.67*data_b5)), count)
if mask_ind ne [-1] then data[mask_ind]=0
;Other
mask_ind=where((data eq 5), count)
if mask_ind ne [-1] then data[mask_ind]=1
;Background
mask_ind=where(data_b1 le 0. or data_b1 gt 1. or $
data_b3 le 0. or data_b3 gt 1. or $
data_b4 le 0. or data_b4 gt 1. or $
data_b5 le 0. or data_b5 gt 1. or $
data_b6 le 0. or data_b6 gt 1., count)
if mask_ind ne [-1] then data[mask_ind]=0
writeu, unit, data
[/QUOTE]
|
|
|
|
Deleted User
New Member
Posts:
| 18 Apr 2013 10:05 AM |
|
|
Hi Mari,
Thanks for your help. I am still having problems getting the button into the Menu. However, I am wondering if it's a compilation issue as I can't see the function in either Classic or 5.0 even though the files are saved to both the save_add and extension directories, respectively.
Also, do you recommend keeping separate PRO files for command line vs. menu functions? I do a lot of work via command line or via envipy, so this effectively doubles the code that I need to maintain.
Thanks,
Guy
|
|
|
|
Deleted User
New Member
Posts:
| 18 Apr 2013 10:05 AM |
|
|
Alex, thanks for the heads up- I greatly appreciate any advice that will speed up my code. :)
|
|
|
|
MariM
Veteran Member
Posts:2396
| 18 Apr 2013 02:16 PM |
|
|
Hi Guy,
I just copied this code and put it in my ENVI 5 Classic save_add directory and it creates the button under the Filter menu. It was designed to work in Classic so may need some changes to work in ENVi 5.
I don't think there should be a problem keeping your code in save_add but it typically is used for adding user functions to ENVI. If you remove the other code from save_add does the menu item show up?
-Mari
|
|
|
|
Deleted User
New Member
Posts:
| 18 Apr 2013 08:40 PM |
|
|
Hi Mari,
I compiled the code using the following four steps:
1. > .full_reset_session
2. Compile PRO file
3. > resolve_all, skip_routines='envi', /continue_on_error, /quiet
4. > save, /routines, filename='C:\Program Files\Exelis\ENVI
\save_add\Cloudmask.sav'
Please note that this compilation procedure works fine for code that I call from Python using envipy.RunTool(), but I wonder if it improperly compiles for functions that need to be in ENVI Classic. I tried using the SAV file on two different machines, one of which was a new install of ENVI 5.0 Classic with nothing new in the save_add directory, and both with the same result.
Thanks,
Guy
|
|
|
|