欢迎您访问程序员文章站本站旨在为大家提供分享程序员计算机编程知识!
您现在的位置是: 首页

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

程序员文章站 2024-03-24 22:32:04
...

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

代码如下

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

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