sentinel数据使用IDL语言出有地理坐标的彩色图

前面提到的使用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