IDL读取GLASS12B02.V02GPP数据，获取地面站点处的值

数据：GLASS12B02.V02.A2009001.2019093.hdf
语言：IDL
任务：读取1982-2017文件夹内的所有8天时间分辨率的GPP遥感图像数据，并找出地面站点经纬度处的像元值，用txt存储，第0-2列为地面站点txt原来的年份，纬度和经度，第3-48列存储遥感图像值。如图：

站点txt文件内容：

代码如下

pro read_GLASS_GPP

  time_0 = systime()
  print,'Start at: ',time_0
  start_time = systime(1)
  compile_opt idl2
  
  input1='E:\碳循环课题组组会\GPP\site_FLUXNET.txt'
  site=read_txt_data_file(input1)
  size_data=size(site,/dimensions)
  result=FLTARR(49,size_data[1,0])
  result[0,*]=site[0,*]
  result[1,*]=site[5,*]
  result[2,*]=site[6,*]
  indir='E:\碳循环课题组组会\GPP'
  ;按照txt中年份遍历，遍历完第一行的46景图像数据，再启动第二行，再去遍历第二行对应年份的46景图像
  for ii = 0,size_data[1,0]-1 do begin
    year = string(uint(site[0,ii]))
    index = 3;第3列-第48列存储第0-第45景图像上的值
    for j = 1,365,8 do begin
      if j lt 10 then begin
        doy='00'+strtrim(string(j),2)
        input=indir+'\'+strtrim(year,2)+'\GLASS12B02.V02.A'+strtrim(year,2)+strtrim(doy,2)+'.2019093.hdf'
      endif
      if (j GT 9) and (j LT 100) then begin
        doy='0'+strtrim(string(j),2)
        input=indir+'\'+strtrim(year,2)+'\GLASS12B02.V02.A'+strtrim(year,2)+strtrim(doy,2)+'.2019093.hdf'
      endif
      if J GT 99 then begin
        doy=strtrim(string(j),2)
        input=indir+'\'+strtrim(year,2)+'\GLASS12B02.V02.A'+strtrim(year,2)+strtrim(doy,2)+'.2019093.hdf'
      endif
      GRID_NAME='GLASS12B02'
      FIELD_NAME='GPP'
      work_dir = 'E:\碳循环课题组组会\GPP\'+strtrim(year,2)
      cd,work_dir
      input = 'GLASS12B02.V02.A'+strtrim(year,2)+strtrim(doy,2)+'.2019093.hdf'
      ;注意，IDL需要及时切换目录，尤其文件不能和根目录放一起，应该先切换文件前的根目录，然后再将文件单独作为一个变量去进行函数或者进程运算。
      ;先前没有进行这样的操作时，在READ_MYD11C1_CMG函数部分会读取文件错误。就是说在input = 'E:\碳循环课题组组会\GPP\1982\GLASS12B02.V02.A1982001.2019093.hdf'时，就会报错
      r_data = READ_MYD11C1_CMG(input, GRID_NAME, FIELD_NAME)
      ;如果没有该景遥感图像，则返回-1，赋值为-9999
      ;如果有该景遥感图像，则函数返回正常的值
      if r_data[0,0] eq -1 then begin
        result[index,ii] = -9999
        index += 1
        ;print,input,' is an invalid file.'
      endif
      if r_data[0,0] ne -1 then begin
        ;x4为像元纬度对应的行号
        x4=fix(abs(90-site[5,ii])/(360.0/7200.0))
        ;y4为像元经度对应的列号
        y4=fix(abs(site[6,ii]+180)/(360.0/7200.0))
        ;第index列（index从3-48）第ii行（ii从0-1351）的值为：
        result[index,ii] = r_data[y4,x4]
        ;print,result[index,ii]
        index += 1;该行该列的下一列存储下一个遥感图像值
      endif
    endfor
  endfor
  ;print,result[3,0]
  ;注意，IDL需要及时切换目录，尤其文件不能和根目录放一起，应该先切换文件前的根目录，然后再将文件单独作为一个变量去进行函数或者进程运算。
  work_dir = 'E:\碳循环课题组组会\GPP\'
  cd,work_dir
  output1='fluxnet_site_Corresponding_to_GPP_extraction.txt'
  openw,SI1,output1,/get_lun
  printf,SI1,result
  free_lun,SI1
  end_time = systime(1)
  print,'耗时：',(end_time-start_time)/60,'  分钟！'
  !null=dialog_message('the task is finished!',/information)
end
FUNCTION read_txt_data_file, infilename
    ;Get the number of lines
    nlines = FILE_LINES(infilename)
    ;指定一个内存号lun1
    OPENR, lun1, infilename, /GET_LUN
    ;存储文本文档的全部行到temp_str
    tmp_str = ''
    READF, lun1, tmp_str
    tmp = STRSPLIT(tmp_str, COUNT = col_count)
    POINT_LUN, lun1, 0
    
    ;Allocate memory
    data = DBLARR(col_count, nlines)
    
    row_count = 0L
    WHILE ~EOF(lun1) DO BEGIN
        READF, lun1, tmp_str
        IF ~STRCMP(tmp_str, '') THEN BEGIN
            tmp_str_split = STRSPLIT(tmp_str, /EXTRACT)
            data_line = DOUBLE(tmp_str_split)
            data[*, row_count] = data_line
            row_count = row_count + 1
        ENDIF        
    ENDWHILE
    FREE_LUN, lun1
    RETURN, data[*, 0 : (row_count - 1)]
END
FUNCTION READ_MYD11C1_CMG, IN_FNAME, GRID_NAME, FIELD_NAME
  compile_opt idl2
  status = EOS_GD_QUERY(in_fname, GRID_NAME, grid_info)
  IF (status NE 1) THEN BEGIN
    ;PRINT, 'Invalid hdf file.'
    RETURN, -1
  ENDIF
  ;IDL函数只有一个return，如果前面执行了return，则后面的全部语句都不再执行
  hd_id = HDF_SD_START(string(in_fname), /READ)
  sd_index = HDF_SD_NAMETOINDEX(hd_id, FIELD_NAME);返回数据集对应的索引号index
  sd_id = HDF_SD_SELECT(hd_id, sd_index);返回该文件该索引数据集的标识符sd_id
  hdf_sd_getdata,sd_id,field_data;原始数据存储在field_data

  ;描述文件整体属性
  ; hdf_sd_fileinfo,hd_id,datasets,attributes;有几个数据集和几个外围属性集
  ;hdf_sd_attrinfo,hd_id,attrindex,name=a,data=b,count=c,hdf_type=d,type=e

  attr_index1 = HDF_SD_ATTRFIND(sd_id, 'scale_factor')
  attr_index2 = HDF_SD_ATTRFIND(sd_id, 'add_offset')
  ;HDF_SD_ATTRINFO,数据集标识符，属性索引，属性值  = 变量名
  ;HDF_SD_ATTRINFO, SD_id, Attr_Index [, COUNT=variable] [, DATA=variable] [, HDF_TYPE=variable] [, NAME=variable] [, TYPE=variable]
  HDF_SD_ATTRINFO, sd_id, attr_index1, DATA = scale_factor
  HDF_SD_ATTRINFO, sd_id, attr_index2, DATA = add_offset
  HDF_SD_ENDACCESS, sd_id
  HDF_SD_END, hd_id
  r_data = field_data * scale_factor[0] + add_offset[0]
  return, r_data
END

---

人的选择有千千万，
我还是那个我，善良，真诚，正直的我! ~