批量读取Landsat九宫格像元DN值的ENVI程序

pro txtpoints

e = ENVI(/headless)

envi,/restore_base_save_files;恢复ENVI sav文件

your_tif=dialog_pickfile(/read,/multiple_files,$

filter=’*_sr_band2.tif’,title=’pleasechoose your tif’,$

path=’J:\Anderson_lab\for_Xudong\Landsat8_surface_reflectance’);读取影像文件路径

n= N_ELEMENTS(your_tif)

data=make_array(n)

data1=make_array(n)

data2=make_array(n)

data3=make_array(n)

data4=make_array(n)

data5=make_array(n)

data6=make_array(n)

data7=make_array(n)

data8=make_array(n)

 

FOR i=0,n-1 do begin;对每个影像，读取指定坐标下的DN值

envi_open_data_file,your_tif(i),/tiff,r_fid=fid  ;  打开影像文件读取DN

print, fid

print,your_tif(i)

xmap=268828.5          ;Easting  Peach

ymap=4037910.5         ;Nothing  Peach

;xmap=268775.3           ;Easting  Almond

;ymap=4037731.3         ;Nothing  Almond

;Landsat是UTM投影的，需要将其度转化为m，http://www.uwgb.edu/dutchs/usefuldata/ConvertUTMNoOZ.HTM

ENVI_CONVERT_FILE_COORDINATES,fid,xf,yf,xmap,ymap;像素坐标到经纬度坐标的快速计算或逆计算

dims=[-1,xf,xf,yf,yf]

data(i)=ENVI_GET_DATA(FID=fid,DIMS=dims,pos=0)

dims1=[-1,xf+1,xf+1,yf,yf]

data1(i)=ENVI_GET_DATA(FID=fid,DIMS=dims1,pos=0)

dims2=[-1,xf-1,xf-1,yf,yf]

data2(i)=ENVI_GET_DATA(FID=fid,DIMS=dims2,pos=0)

dims3=[-1,xf,xf,yf+1,yf+1]

data3(i)=ENVI_GET_DATA(FID=fid,DIMS=dims3,pos=0)

dims4=[-1,xf+1,xf+1,yf+1,yf+1]

data4(i)=ENVI_GET_DATA(FID=fid,DIMS=dims4,pos=0)

dims5=[-1,xf-1,xf-1,yf+1,yf+1]

data5(i)=ENVI_GET_DATA(FID=fid,DIMS=dims5,pos=0)

dims6=[-1,xf,xf,yf-1,yf-1]

data6(i)=ENVI_GET_DATA(FID=fid,DIMS=dims6,pos=0)

dims7=[-1,xf+1,xf+1,yf-1,yf-1]

data7(i)=ENVI_GET_DATA(FID=fid,DIMS=dims7,pos=0)

dims8=[-1,xf-1,xf-1,yf-1,yf-1]

data8(i)=ENVI_GET_DATA(FID=fid,DIMS=dims8,pos=0)

print,i,data(i),data1(i),data2(i),data3(i),data4(i),data5(i),data6(i),data7(i),data8(i)

ENDFOR

;;将DN值写出来

fname=’F:\test.txt’

openw,1,fname

FOR i=0,n-1 do begin;对每个影像，读取指定坐标下的DN值

printf,1,data(i),data1(i),data2(i),data3(i),data4(i),data5(i),data6(i),data7(i),data8(i),format='(9a)’

ENDFOR

print,’completed’

close,1

end