前面提到的使用IDL代码可以将sentinel数据做一个加坐标系的功能,这个代码是自己编写的,因为刚开始已经把数据导入进了envi,并且顺便导出成了没有坐标系的tiff格式。(实际上并不用,可以直接不用导出,全自动的读取,只是后边的读取t,read_tiff编程对应的read_JP2000就行)这个程序我是使用了10m的8波段和20m的11、12波段,你们可以自己选择。
由于此次出图的量比较少,所以代码不是很完善,能偷懒的就偷懒了,如果你们有需要用到的,需要改几个地方:
1、如果你们想直接读取JP2,不想做导出这一步,需要更改路径,具体可以查看我的上一篇文章,路径在哪里。
2、读取xml数据的时候我其实只读取了角点的信息,没有读取对应的分辨率、中心经线信息,如果有需要,可以按照读取角点信息更改代码,把几种数据读全
3、输出的是彩色图像,需要可以自行输出,输出数据为三维,则为彩色图像,二维则是灰度图像
pro read_tiff1
dir='H:\sentinel2_L2\geo\'
file_list=file_search(dir,'*.tif')
file_n=n_elements(file_list)
file_list1=[]
for file_i=0,file_n-1 do begin
file_name=file_basename(file_list[file_i],'.tif')
file_name=strsplit(file_name,'_',/extract)
name_now=file_name[0]+'_'+file_name[1]
if where(file_list1 eq name_now) eq -1 then file_list1=[file_list1,name_now]
endfor
data_n=n_elements(file_list1)
for data_i=0,data_n-1 do begin
;data_need=[]
data_name=file_list1[data_i]
data_b11_dir=file_search(dir,data_name+'_B11'+'*.tif')
data_b12_dir=file_search(dir,data_name+'_B12'+'*.tif')
data_b8_dir=file_search(dir,data_name+'_B08'+'*.tif')
data_b11=read_tiff(data_b11_dir)
data_b11=data_b11*0.0001
data_b12=read_tiff(data_b12_dir)
data_b12=data_b12*0.0001
data_b8=read_tiff(data_b8_dir)
data_b8=data_b8*0.0001
size_data=size(data_b8)
need_col=size_data[1]
need_row=size_data[2]
data_b11=congrid(data_b11,need_col,need_row)
data_b12=congrid(data_b12,need_col,need_row)
out_data=fltarr(3,need_col,need_row)
out_data[0,*,*]=data_b12
out_data[1,*,*]=data_b11
out_data[2,*,*]=data_b8
data_b12=!null
data_b11=!null
data_b8=!null
;no_cloud=data_b2 le 0.2
;data_b4=data_b4*no_cloud
;data_b8=data_b8*no_cloud
;ndvi_now=float(data_b8-data_b4)/float(data_b8+data_b4)
;ndvi_now[where(~finite(ndvi_now))] = -9999
;ndvi_now=(ndvi_now ne -9999)*ndvi_now
name_dir=strmid(data_name,0,12)
file_xml_dir=file_search('H:\sentinel2_L2\','*'+string(name_dir)+'*'+'.SAFE',/test_directory);
;xml_dir_pos=where(file_xml_dir eq '*'+data_name+'.SAFE')
txt_dir=file_search(file_xml_dir,'MTD_TL.xml');file_xml_dir+'\'+'GRANULE\L2A_T47'+name_dir+'*'+'\MTD_TL.xml'
openr,1,txt_dir
lines=file_lines(txt_dir)
text1=strarr(1,lines)
readf,1,text1
for linei=0L,lines[0]-1 do begin
temp=text1[0,linei]
temp=temp[0]
result=strmatch(temp,'*<ULX>*',/FOLD_CASE);strmatch使用*,为模糊查询
if result eq 1 then begin;[where(strmatch(strtrim(temp,13),'*<ULX>') eq 1)]
line=linei
break
endif
;print,11
endfor
line_text=text1[0,line]
line_text1=text1[0,line+1]
x_text=line_text[0];先读取所在的行
y_text=line_text1[0]
x_num=float(strmid(x_text,13,6));读取x的值
y_num=float(strmid(y_text,13,7))
; tiaodai=text1[*,16]
; tiaodai=string(tiaodai[0])
; num=float(strmid(tiaodai,31,5))
; x_text=where(text1 eq '*<HORIZONTAL_CS_CODE>EPSG*')
; x_text=x_text[0]
; y_text=text1[*,31]
; y_text=y_text[0]
; ;resoulution_text=text1[*,29]
; ;resoulution_text=resoulution_text[0]
; ;resolution=float(strmid(resoulution_text,31,2))
; x_num=float(strmid(x_text,13,6))
; y_num=float(strmid(y_text,13,7))
;print,x_text
out_dir='H:\sentinel2_L2\out_TC\'+data_name
;print,tiaodai[8:11]
free_lun,1
; print,resolution
; print,x_num,y_num
; print,num
geo_info={MODELPIXELSCALETAG:[10.0,10.0,0.0],$
MODELTIEPOINTTAG:[0.0,0.0,0.0,x_num,y_num,0.0] ,$
GTMODELTYPEGEOKEY:1,$
GTRASTERTYPEGEOKEY:1,$
GEOGLINEARUNITSGEOKEY:9001,$
GEOGANGULARUNITSGEOKEY:9102,$
PROJECTEDCSTYPEGEOKEY:32647}
; write_tiff,out_dir+'_b4'+'.tiff',data_b4,geotiff=geo_info,/float
; write_tiff,out_dir+'_b8'+'.tiff',data_b8,geotiff=geo_info,/float
write_tiff,out_dir+'.tiff',out_data,geotiff=geo_info,/float
print,111
endfor
end