| How to apply a vegetation index on multiple images? [message #89626] |
Sun, 02 November 2014 10:06  |
Muhammad Hasan
Messages: 8
Registered: November 2014
|
Junior Member |
|
|
Hi All,
I have a vegetation index (Tasseled cap transformation). I have 10 folders for different areas and each folder has 20 Landsat images. My friend helped me to write following code:
pro LCTCT
infile=dialog_pickfile(title='Select you all your data!!!',filter='*.tif',get_path=path,/multiple_files)
if N_elements(infile) eq -1 then return
cd, path
filenum=size(infile,/n_elements)
for i=0,filenum-1 do begin
envi_open_file, infile(i), r_fid=fid
ENVI_FILE_QUERY,fid,dims=dims,ns=ns,nl=nl,nb=nb,data_type=da ta_type
map_info=envi_get_map_info(fid=fid)
b1=envi_get_data(fid=fid,dims=dims,pos=0)
b2=envi_get_data(fid=fid,dims=dims,pos=1)
b3=envi_get_data(fid=fid,dims=dims,pos=2)
b4=envi_get_data(fid=fid,dims=dims,pos=3)
b5=envi_get_data(fid=fid,dims=dims,pos=4)
b6=envi_get_data(fid=fid,dims=dims,pos=5)
tct=fltarr(ns,nl,6)
print,ns,nl,nb
file = DIALOG_PICKFILE(FILTER='*.txt')
OPENR, lun, file, /GET_LUN
para=fltarr(6,6)
WHILE NOT EOF(lun) DO BEGIN & $
READF, lun, para
ENDWHILE
; Close the file and free the file unit
FREE_LUN, lun
tct(*,*,0)=float((para(0,0)*b1)+(para(1,0)*b2)+(para(2,0)*b3 )+( para(3,0)*b4)+(para(4,0)*b5)+(para(5,0)*b6))
tct(*,*,1)=float((para(0,1)*b1)+(para(1,1)*b2)+(para(2,1)*b3 )+( para(3,1)*b4)+(para(4,1)*b5)+(para(5,1)*b6))
tct(*,*,2)=float((para(0,2)*b1)+(para(1,2)*b2)+(para(2,2)*b3 )+( para(3,2)*b4)+(para(4,2)*b5)+(para(5,2)*b6))
tct(*,*,3)=float((para(0,3)*b1)+(para(1,3)*b2)+(para(2,3)*b3 )+( para(3,3)*b4)+(para(4,3)*b5)+(para(5,3)*b6))
tct(*,*,4)=float((para(0,4)*b1)+(para(1,4)*b2)+(para(2,4)*b3 )+( para(3,4)*b4)+(para(4,4)*b5)+(para(5,4)*b6))
tct(*,*,5)=float((para(0,5)*b1)+(para(1,5)*b2)+(para(2,5)*b3 )+( para(3,5)*b4)+(para(4,5)*b5)+(para(5,5)*b6))
help,tct
envi_enter_data,tct,map_info=map_info,r_fid=outfid
filename = FILE_BASENAME(infile(i))
outfiledir=file_dirname(infile(i),/MARK_DIRECTORY)
pointPos = STRPOS(filename,'.',/REVERSE_SEARCH)
IF pointPos[0] NE -1 THEN BEGIN
filename= STRMID(filename,0,pointPos)
ENDIF
outfilename = outfiledir+filename+'_NewTCT_fromTOAinsteadofSR.tif'
envi_output_to_external_format,dims = dims,pos=[0,1,2,3,4,5],/tiff,fid = outfid,out_name=outfilename
envi_file_mng,id = fid,/remove
envi_file_mng,id = outfid,/remove,/delete
endfor
print,'ok'
end
What should I do to make this code capable of computing TCT index automatically? (I am not good at programming)
Regards,
|
|
|
|
| Re: How to apply a vegetation index on multiple images? [message #89634 is a reply to message #89626] |
Mon, 03 November 2014 07:46  |
Phillip Bitzer
Messages: 223
Registered: June 2006
|
Senior Member |
|
|
On Sunday, November 2, 2014 12:06:02 PM UTC-6, Muhammad Hasan wrote:
> Hi All,
> I have a vegetation index (Tasseled cap transformation). I have 10 folders for different areas and each folder has 20 Landsat images. My friend helped me to write following code:
I think you need to have a talk with your friend. There are a few things that should be cleaned up. I'll highlight a few lines. One of the nice things about IDL is that it is relatively flexible - you can do things the wrong way and still the right answer. But, it's good to learn the right way now!
> if N_elements(infile) eq -1 then return
I don't think this will ever be true. If this is supposed to catch when the user cancels out the dialog, you should check for the null string.
> cd, path
Honestly, I hate programs that change my working directory (and don't put it back), but it doesn't really hurt you. (But watch out if you start to save things like images/plots/etc - the files will likely end up somewhere you're not looking.)
> envi_open_file, infile(i), r_fid=fid
The first of many like this. Using () to index arrays is very old school and is considered wrong now. Use brackets [] to index arrays. See http://exelisvis.com/docs/Understanding_Array_Subs.html
> file = DIALOG_PICKFILE(FILTER='*.txt')
So, you have to pick a file for every infile? Which means you have to do this 200 times? Crikey, that would drive me crazy.
> WHILE NOT EOF(lun) DO BEGIN & $
Why is the &$ there?
> READF, lun, para
> ENDWHILE
So, your while loop goes through the entire file. At every iteration, it reads in a 6x6 float array...until it reaches the end. The variable para will only be the *last* 6x6 array in the file when you move the next step, Is that what you want?
> tct(*,*,0)=float((para(0,0)*b1)+(para(1,0)*b2)+(para(2,0)*b3 )+( para(3,0)*b4)+(para(4,0)*b5)+(para(5,0)*b6))
I'm quite sure there are more efficient ways to do this. This looks like a matrix multiplication problem.
> filename = FILE_BASENAME(infile(i))
> outfiledir=file_dirname(infile(i),/MARK_DIRECTORY)
> pointPos = STRPOS(filename,'.',/REVERSE_SEARCH)
> IF pointPos[0] NE -1 THEN BEGIN
Hey! Square bracket! Alright! ;-)
> What should I do to make this code capable of computing TCT index automatically? (I am not good at programming)
I haven't the foggiest idea what a TCT index is. You'll have to tell us what it is before we can calculate it.
There are a lot of good resources out there to help you get better at programming. This newsgroup, IDL Data point blog, the online help from Excelsis, Coyote's website (idlcoyote.com), and others I'm not thinking of right now, I'ms ure.
|
|
|
|
| Re: How to apply a vegetation index on multiple images? [message #89667 is a reply to message #89634] |
Sun, 09 November 2014 06:32 |
Muhammad Hasan
Messages: 8
Registered: November 2014
|
Junior Member |
|
|
On Monday, November 3, 2014 11:46:40 PM UTC+8, Phillip Bitzer wrote:
> On Sunday, November 2, 2014 12:06:02 PM UTC-6, Muhammad Hasan wrote:
>> Hi All,
>> I have a vegetation index (Tasseled cap transformation). I have 10 folders for different areas and each folder has 20 Landsat images. My friend helped me to write following code:
>
> I think you need to have a talk with your friend. There are a few things that should be cleaned up. I'll highlight a few lines. One of the nice things about IDL is that it is relatively flexible - you can do things the wrong way and still the right answer. But, it's good to learn the right way now!
>
>> if N_elements(infile) eq -1 then return
>
> I don't think this will ever be true. If this is supposed to catch when the user cancels out the dialog, you should check for the null string.
>
>> cd, path
>
> Honestly, I hate programs that change my working directory (and don't put it back), but it doesn't really hurt you. (But watch out if you start to save things like images/plots/etc - the files will likely end up somewhere you're not looking.)
>
>> envi_open_file, infile(i), r_fid=fid
>
> The first of many like this. Using () to index arrays is very old school and is considered wrong now. Use brackets [] to index arrays. See http://exelisvis.com/docs/Understanding_Array_Subs.html
>
>> file = DIALOG_PICKFILE(FILTER='*.txt')
>
> So, you have to pick a file for every infile? Which means you have to do this 200 times? Crikey, that would drive me crazy.
>
>> WHILE NOT EOF(lun) DO BEGIN & $
>
> Why is the &$ there?
>
>> READF, lun, para
>> ENDWHILE
>
> So, your while loop goes through the entire file. At every iteration, it reads in a 6x6 float array...until it reaches the end. The variable para will only be the *last* 6x6 array in the file when you move the next step, Is that what you want?
>
>> tct(*,*,0)=float((para(0,0)*b1)+(para(1,0)*b2)+(para(2,0)*b3 )+( para(3,0)*b4)+(para(4,0)*b5)+(para(5,0)*b6))
>
> I'm quite sure there are more efficient ways to do this. This looks like a matrix multiplication problem.
>
>> filename = FILE_BASENAME(infile(i))
>> outfiledir=file_dirname(infile(i),/MARK_DIRECTORY)
>> pointPos = STRPOS(filename,'.',/REVERSE_SEARCH)
>> IF pointPos[0] NE -1 THEN BEGIN
>
> Hey! Square bracket! Alright! ;-)
>
>> What should I do to make this code capable of computing TCT index automatically? (I am not good at programming)
>
> I haven't the foggiest idea what a TCT index is. You'll have to tell us what it is before we can calculate it.
>
> There are a lot of good resources out there to help you get better at programming. This newsgroup, IDL Data point blog, the online help from Excelsis, Coyote's website (idlcoyote.com), and others I'm not thinking of right now, I'ms ure.
Dear Philip Bitzer,
I owe debt of gratitude to you for correcting so many mistakes in my code. I am sorry for late reply as I was very busy in writing my thesis which kept me away from frequent use of internet. So nice of you for this kind help.
Regards,
|
|
|
|
comp.lang.idl-pvwave archive
Members
Search
Help
Login
Home
