前言

    近期开展基于GEDI数据的研究，以期得到基于星载激光雷达的森林蓄积量等林分参数的反演结果，在此对近期工作进行记录。

图一 GEDI硬件

数据下载

    有关GEDI数据的介绍可通过官网（https://gedi.umd.edu/）进行查询，在此不多做赘述。小道消息：在Earthdata进行批量数据下载时若出现‘登录不成功’、‘未知错误’、‘网络报错’等情况，可通过手机热点进行下载，若不然只能通过一次次的单击下载（笔者因流量问题未曾尝试使用该方法进行数据下载）。

L2B 数据提取植被指数

        下载完成后的原始 .h5文件通过IO操作写成csv文件。（批量压缩包解压可通过：右键---解压到单独文件夹）

        该部分的源码可通过GitHub搜索‘rGEDI’获取。本次实验所使用的源码如下：

# (GEDI Level2B)

library(rGEDI)
library(sp)

filedir = 'E:/ws/r/GEDI/L2B'
outdir = 'E:/ws/r/GEDI/L2B_O'
filelist = list.files(path = filedir, pattern = '*.h5')

# Do not change anything below unless ur brave
# h5 file to metrics
h5toMetrix = function(x){
  level2Bpath = paste0(filedir,'\\',x,'.h5')
  gedilevel2b = readLevel2B(level2Bpath)
  
  #Get GEDI Vegetation Biophysical Variables
  level2BVPM<-getLevel2BVPM(gedilevel2b)
  #head(level2BVPM[,c("beam","shot_number","pai","fhd_normal","omega","pgap_theta","cover")])
  # Converting shot_number as "integer64" to "character"
  level2BVPM$shot_number<-paste0(level2BVPM$shot_number)
  # Converting GEDI Vegetation Profile Biophysical Variables as data.table to SpatialPointsDataFrame
  #level2BVPM_spdf<-SpatialPointsDataFrame(cbind(level2BVPM$longitude_lastbin,level2BVPM$latitude_lastbin),data=level2BVPM)
  # Exporting GEDI Vegetation Profile Biophysical Variables as ESRI Shapefile
  #raster::shapefile(level2BVPM_spdf,paste0(outdir,'\\',filename))
  
  
  #Get Plant Area Index (PAI)
  level2BPAIProfile<-getLevel2BPAIProfile(gedilevel2b)
  # Converting shot_number as "integer64" to "character"
  level2BPAIProfile$shot_number<-paste0(level2BPAIProfile$shot_number)
  
  #Get Plant Area Volume Density (PAVD) Profiles
  level2BPAVDProfile<-getLevel2BPAVDProfile(gedilevel2b)
  # Converting shot_number as "integer64" to "character"
  level2BPAVDProfile$shot_number<-paste0(level2BPAVDProfile$shot_number)
  
  # save the result to a new dataframe
  re_head = data.frame(
    beam = level2BPAIProfile$beam,
    shot_number = level2BPAIProfile$shot_number,
    long = level2BPAIProfile$lon_lowestmode,
    lat = level2BPAIProfile$lat_lowestmode,
    first_return = level2BVPM$elev_highestreturn,
    last_return = level2BVPM$elev_lowestmode
  )
  
  VPM = level2BVPM[,14:19]
  PAI = level2BPAIProfile[,12:26]
  PAVD = level2BPAVDProfile[,12:26]
  
  re = cbind(re_head,VPM,PAI,PAVD)
  
  write.csv(re,file = paste0(outdir,'/',x,'_Metrix','.csv'))
}
er_file = data.frame()
for (i in 1:length(filelist)) {
  file = substring(filelist[i],1,nchar(filelist[1])-3)
  print(paste0('====',file,' Start running===='))
  
  tryCatch(
    {
      h5toMetrix(file)
      
    },error = function(e) {
      err = file
      rbind(er_file,err)
    }
    )
}

悄悄咪咪的把源文件放在最后面 ： h5tometrix.R 嘿嘿嘿