Final extract
乐正秦斩
pro tib_bt_tiqu
COMPILE_OPT IDL2,hidden
e=envi(/headless)
path='E:\data\BT\bt_500'
;path='E:\data\BT\te_st'
wendu=indgen(20)+1;温度
;判断日期是第几天转化为温度数组下标 写入
fn='C:\Users\Administrator\Desktop\tibtation.csv';
data=READ_CSV(fn,count=nsta,header=header);
point_fid=data.(0);
Lat=data.(2);站点纬度
Lon=data.(3);站点经度
station_id=data.(1);站点编号数组
files = file_search(path,'*.tif',count=n)
;print,files
;二维数组存储
;nb=10波段数 n个文件 nsta站点数 一张图要提取102个站点的bt n张图*nsta站点数
h=fltarr(10,nsta*n)
;编号
id=indgen(nsta*n)+1
;温度
wd=fltarr(nsta*n)
;dem 地表起伏度
dem=fltarr(nsta*n)
;日期
datetime=strarr(nsta*n)
;站点编号
station=lonarr(nsta*n)
;实测雪深
sd_true=fltarr(nsta*n)
;LC
lc=intarr(nsta*n)
;经纬度
jingdu=fltarr(nsta*n)
weidu=fltarr(nsta*n)
;坡度
podu=fltarr(nsta*n)
;坡向
poxiang=intarr(nsta*n)
;粗糙度
cucaodu=fltarr(nsta*n)
position = 0
for k=0,n-1 do begin
file = files[k]
;print,file
;实际需要的文件是中国区域裁剪重采样500m的数据
;image_fname='E:\data\xinout\FY3B_MWRID_GBAL_L1_20141214_1950_010KM_MS_result.tif';文件路径
;打开文件
ENVI_OPEN_FILE,file,r_fid=fid
;print ,fid
ENVI_FILE_QUERY, fid, dims=dims, nb=nb,bnames = bnames,DATA_TYPE = dt,ns = ns, nl = nl;nb波段数
map_info=ENVI_GET_MAP_INFO(fid=fid)
;print,map_info
i_proj=ENVI_PROJ_CREATE(/geographic,datum='WGS-84')
o_proj=map_info.proj
envi_convert_projection_coordinates,lon,lat,i_proj,$
xmap,ymap,o_proj
;print,lon,lat,i_proj,xmap,ymap,o_proj
;将站点的地图坐标转换为文件坐标
envi_convert_file_coordinates,fid,xf,yf,xmap,ymap
xf=floor(xf)
yf=floor(yf)
;E:\data\BT\bt_500\FY3B_MWRID_GBAL_L1_20130102_2020_010KM_MS_result.tif
;E:\data\wendu_lc_500_dem\A20121220.tif
;根据bt的日期构造温度文件名，打开，站点坐标转温度文件坐标 xf_t yf_t 文件号fid_t
;下面站点循环构造dim_t
t_name = strsplit(file,'_',/extract)
t_file = 'E:\data\wendu_lc_500_dem\A'+t_name[5]+'.tif'
envi_open_file,t_file,r_fid=fid_t
envi_file_query,fid_t,dims=dims_t
map_info_t=ENVI_GET_MAP_INFO(fid=fid_t)
;print,map_info
;i_proj=ENVI_PROJ_CREATE(/geographic,datum='WGS-84')
o_proj_t=map_info_t.proj
envi_convert_projection_coordinates,lon,lat,i_proj,$
xmap_t,ymap_t,o_proj_t
;print,lon,lat,i_proj,xmap,ymap,o_proj
;将站点的地图坐标转换为文件坐标
envi_convert_file_coordinates,fid_t,xf_t,yf_t,xmap_t,ymap_t
xf_t=floor(xf_t)
yf_t=floor(yf_t)
;打开lc、dem，站点坐标分别转文件坐标 xf_lc yf_lc xf_dem yf_dem 文件号fid_lc fid_dem
;下面站点循环构造dim_lc dim_dem
;这里只提取地形起伏度
dem_file = 'C:\Users\Administrator\Desktop\arc\fluctuation.tif'
envi_open_file,dem_file,r_fid=fid_dem
envi_file_query,fid_dem,dims=dims_dem
map_info_dem=ENVI_GET_MAP_INFO(fid=fid_dem)
;print,map_info
;i_proj=ENVI_PROJ_CREATE(/geographic,datum='WGS-84')
o_proj_dem=map_info_dem.proj
envi_convert_projection_coordinates,lon,lat,i_proj,$
xmap_dem,ymap_dem,o_proj_dem
;print,lon,lat,i_proj,xmap,ymap,o_proj
;将站点的地图坐标转换为文件坐标
envi_convert_file_coordinates,fid_dem,xf_dem,yf_dem,xmap_dem,ymap_dem
xf_dem=floor(xf_dem)
yf_dem=floor(yf_dem)
;构造podu、pixiang、cucaodu的dims
podu_file = 'C:\Users\Administrator\Desktop\arc\touyingslope.tif'
poxiang_file = 'C:\Users\Administrator\Desktop\arc\poxiang_chongfenlei.tif'
cucaodu_file = 'C:\Users\Administrator\Desktop\arc\cucaodu.tif'
envi_open_file,podu_file,r_fid = fid_podu
envi_file_query,fid_podu,dims=dims_podu
envi_open_file,poxiang_file,r_fid = fid_poxiang
envi_file_query,fid_poxiang,dims=dims_poxiang
envi_open_file,cucaodu_file,r_fid = fid_cucaodu
envi_file_query,fid_cucaodu,dims=dims_cucaodu
map_info_podu=envi_get_map_info(fid=fid_podu)
o_proj_podu=map_info_podu.proj
envi_convert_projection_coordinates,lon,lat,i_proj,$
xmap_podu,ymap_podu,o_proj_podu
envi_convert_file_coordinates,fid_podu,xf_podu,yf_podu,xmap_podu,ymap_podu
xf_podu=floor(xf_podu)
yf_podu=floor(yf_podu)
map_info_poxiang=envi_get_map_info(fid=fid_poxiang)
o_proj_poxiang=map_info_poxiang.proj
envi_convert_projection_coordinates,lon,lat,i_proj,$
xmap_poxiang,ymap_poxiang,o_proj_poxiang
envi_convert_file_coordinates,fid_poxiang,xf_poxiang,yf_poxiang,xmap_poxiang,ymap_poxiang
xf_poxiang=floor(xf_poxiang)
yf_poxiang=floor(yf_poxiang)
map_info_cucaodu=envi_get_map_info(fid=fid_cucaodu)
o_proj_cucaodu=map_info_cucaodu.proj
envi_convert_projection_coordinates,lon,lat,i_proj,$
xmap_cucaodu,ymap_cucaodu,o_proj_cucaodu
envi_convert_file_coordinates,fid_cucaodu,xf_cucaodu,yf_cucaodu,xmap_cucaodu,ymap_cucaodu
xf_cucaodu=floor(xf_cucaodu)
yf_cucaodu=floor(yf_cucaodu)
;构造lc的dims
lc_file = 'E:\data\lc.tif'
envi_open_file,lc_file,r_fid = fid_lc
envi_file_query,fid_lc,dims=dims_lc
map_info_lc=envi_get_map_info(fid=fid_lc)
o_proj_lc=map_info_lc.proj
envi_convert_projection_coordinates,lon,lat,i_proj,$
xmap_lc,ymap_lc,o_proj_lc
envi_convert_file_coordinates,fid_lc,xf_lc,yf_lc,xmap_lc,ymap_lc
xf_lc=floor(xf_lc)
yf_lc=floor(yf_lc)
;循环读取各个站点dn
;print,i bug在此
for i=0, nsta-1 do begin;站点
;bt的dims
dims1=[-1, xf[i], xf[i], yf[i], yf[i]] ;构建空间范围数组DIMS
for j=0 , nb-1 do begin;bt有十个波段
temp = envi_get_data(fid=fid,dims=dims1,pos=j)
if temp eq 0 or temp eq -999 then begin;-999无效值
;h[j,i*k]=0 ;不能乘 0-0 0-102 0 2 4 - 204 0 3 6 9 - 306只能计数
h[j,position]=0
endif else begin
h[j,position]=temp*0.01+327.68
endelse
endfor
;温度的dims
dims2=[-1, xf_t[i], xf_t[i], yf_t[i], yf_t[i]]
temp_t = envi_get_data(fid=fid_t,dims=dims2,pos=0);只有一个波段
;这里可以对温度无效值进行判断
wd[position]=temp_t
;地形的dims
dims3=[-1, xf_dem[i], xf_dem[i], yf_dem[i], yf_dem[i]]
temp_dem = envi_get_data(fid=fid_dem,dims=dims3,pos=0);只有一个波段
;这里可以对地形无效值进行判断
dem[position]=temp_dem
;日期数组
datetime[position]=t_name[5]
;strmid(t_name[5],0);A20121220
;LC数组
dims4=[-1, xf_lc[i], xf_lc[i], yf_lc[i], yf_lc[i]]
temp_lc = envi_get_data(fid=fid_lc,dims=dims4,pos=0);只有一个波段
;这里可以对地形无效值进行判断
lc[position]=temp_lc
;podu数组
dims5=[-1, xf_podu[i], xf_podu[i], yf_podu[i], yf_podu[i]]
temp_podu = envi_get_data(fid=fid_podu,dims=dims5,pos=0);只有一个波段
;这里可以对地形无效值进行判断
podu[position]=temp_podu
;poxiang数组
dims6=[-1, xf_poxiang[i], xf_poxiang[i], yf_poxiang[i], yf_poxiang[i]]
temp_poxiang = envi_get_data(fid=fid_poxiang,dims=dims6,pos=0);只有一个波段
;这里可以对地形无效值进行判断
poxiang[position]=temp_poxiang
;cucaodu数组
dims7=[-1, xf_cucaodu[i], xf_cucaodu[i], yf_cucaodu[i], yf_cucaodu[i]]
temp_cucaodu = envi_get_data(fid=fid_cucaodu,dims=dims7,pos=0);只有一个波段
;这里可以对地形无效值进行判断
cucaodu[position]=temp_cucaodu
;经纬度数组
jingdu[position]=Lon[i]
weidu[position]=Lat[i]
;站点编号数组
station[position]=station_id[i]
;print,station_id[i]
;实测雪深数组
;根据站点编号构造文件名
sd_file = 'C:\Users\Administrator\Desktop\青藏高原雪深数据集1961-2013\青藏高原雪深数据集1961-2013\Tibetan_Plateau_Snow_Depth_1961-2013\'+string(station_id[i],format='(i5)')+'.txt'
;print,sd_file
;根据时间找出那一行的雪深测量值 20121220
time = t_name[5]
year = fix(strmid(time,0,4))
month = fix(strmid(time,4,2))
day = fix(strmid(time,6,2))
;读取txt文件 逐行判断（效率好低）
cnt = file_lines(sd_file)
openr,lun,sd_file,/get_lun
;bug w
line=''
for w = 0 , cnt-1 do begin
;openr,lun,sd_file,/get_lun
; readf,lun,a,y,m,d,dd
;
; if fix(y) eq year and fix(m) eq month and fix(d) eq day then begin
; print,fix(y),fix(m),fix(d),dd
; sd_true[position]=dd
; endif
;free_lun,lun
readf,lun,line
out = strsplit(line[0],/extract)
if fix(out[1]) eq year and fix(out[2]) eq month and fix(out[3]) eq day then begin
if n_elements(out) eq 5 then begin
sd_true[position]=float(out[4])
print,float(out[4])
endif else begin
sd_true[position]=-999
print,-999
endelse
endif
endfor
;print,lun
free_lun,lun
; riqi = strsplit(file,'_',/extract)
; riqitonum = long(riqi[4])
; ;print,riqitonum
; wd[pos]=wendu[riqitonum-20121215]
; print,riqitonum-20121215
; datatime[pos]=riqi[4]+riqi[5]
position +=1
endfor
print,file+"-----提取结束！"
endfor
;print,h[0,*];这一天所有站点第一个波段的值
;print,h[1,*];这一天所有站点第二个波段的值
;print,"********************"
;print,h[*,0];这一天第一个站点所有波段的值
;***********保存结果*********
o_fn='E:\data\tiqu.csv';指定存储结果路径
;header=[header,'h1','h2','h3','h4','h5','h6','h7','h8','h9','h10']
header=['id','station_id','jingdu','weidu','datetime','10V','10H','18V','18H','23V','23H','36V','36H','89V','89H','10H23H','10V23H','10V23V','18H23H','18V23H','23V23H','10V36H','36V36H','wendu','fluctuation','podu','poxiang','cucaodu','LC','sd']
;data=create_struct(data,'H1',h[0,*],'H2',h[1,*],'H3',h[2,*],'H4',h[3,*],'H5',h[4,*],'H6',h[5,*],'H7',h[6,*],'H8',h[7,*],'H9',h[8,*],'H10',h[9,*]);在data中增加一个value域
; help,h[1,*]-h[5,*]
; help,abs(h[0,*]-h[5,*])
; help,h[0,*]
data=create_struct('id',id,'station_id',station,'jingdu',jingdu,'weidu',weidu,'dt',datetime,'H1',h[0,*],'H2',h[1,*],'H3',h[2,*],'H4',h[3,*],'H5',h[4,*],'H6',h[5,*],'H7',h[6,*],'H8',h[7,*],'H9',h[8,*],'H10',h[9,*],'_10H23H',abs(h[1,*]-h[5,*]),'_10V23H',abs(h[0,*]-h[5,*]),'_10V23V',abs(h[0,*]-h[4,*]),'_18H23H',abs(h[3,*]-h[5,*]),'_18V23H',abs(h[2,*]-h[5,*]),'_23V23H',abs(h[4,*]-h[5,*]),'_10V36H',abs(h[0,*]-h[7,*]),'_36V36H',abs(h[6,*]-h[7,*]),'wd',wd,'fluctuation',dem,'podu',podu,'poxiang',poxiang,'cucaodu',cucaodu,'lc',lc,'sd',sd_true);在data中增加一个value域
write_csv,o_fn,data,header=header;写入csv
;e.close
end
类似资料：