获取HDF文件中的波段平均值并另存为geotiff

Question

我想要得到一个HDF文件的平均波段，并保存为geotiff。具体地说，MODIS MCD18A1有15个波段)，我想要“均值(波段4:波段11)”。

MCD18A1.A2001045.h24v06.006.20171801358.hdf:MODISRAD:GMT_0000_DSR“MCD18A1.A2001045.h24v06.006.20171801358.hdf:MODISRAD:GMT_0300_DSR”MCD18A1.A2001045.h24v06.006.20171801358.hdf:MODISRAD:GMT_0600_DSR“MCD18A1.A2001045.h24v06.006.20171801358.hdf:MODISRAD:GMT_0900_DSR”MCD18A1.A2001045.h24v06.006.20171801358.hdf:MODISRAD:GMT_1200_DSR“MCD18A1.A2001045.h24v06.006.20171801358.hdf:MODISRAD:GMT_1500_DSR”MCD18A1.A2001045.h24v06.006.20171801358.hdf:MODISRAD:GMT_1800_DSR“MCD18A1.A2001045.h24v06.006.20171801358.hdf:MODISRAD:GMT_2100_DSR

这些是MODIS格式的MCD18A1短波辐射产品中的波段，我想取这些波段的平均值，然后提取/保存为"geotiff“。

我可以提取单波段信息，但如何处理HDF栅格中的多波段均值？

=============================================================文件对应于MCD18A1 HDF磁贴。

for (i in 1:length(files)){
  sds <- get_subdatasets(files[i])
  gdal_translate(sds[1],
    dst_dataset = paste0(basename(files[i]),".tif"),  
    of = "GTiff")  
} 

==============================================================

谢谢

Answer 1

我用luna (github)下载了一个示例文件，并展示了如何用terra (CRAN)处理它

library(luna)
e <- c(6.04,  6.23, 49.72, 49.91)
# for this you need to sign up with NASA for a username and pwd. It is free
f <- getModis("MCD18A1", "2020-01-01", "2020-01-01", e, path=".",
    download=TRUE, username="*****", password="*****")
f
#"./MCD18A1.A2020001.h18v04.006.2020050045103.hdf"

假设您已经有了一个文件f，您可以跳过上面的步骤

library(terra)
r <- rast(f)
#class       : SpatRaster 
#dimensions  : 240, 240, 21  (nrow, ncol, nlyr)
#resolution  : 4633.127, 4633.127  (x, y)
#extent      : 0, 1111951, 4447802, 5559753  (xmin, xmax, ymin, ymax)
#coord. ref. : +proj=sinu +lon_0=0 +x_0=0 +y_0=0 +ellps=clrk66 +units=m +no_defs 
#sources     : MCD18A1.A2020001.h18v04.006.2020050045103.hdf:MODISRAD:DSR  (4 layers) 
#              MCD18A1.A2020001.h18v04.006.2020050045103.hdf:MODISRAD:Direct  (4 layers) 
#              MCD18A1.A2020001.h18v04.006.2020050045103.hdf:MODISRAD:Diffuse  (4 layers) 
#              ... and 9 more source(s)
#names       : MODIS~:DSR_1, MODIS~:DSR_2, MODIS~:DSR_3, MODIS~:DSR_4, MODIS~rect_1, MODIS~rect_2, ... 

names(r)
# [1] "MODISRAD:DSR_1"        "MODISRAD:DSR_2"        "MODISRAD:DSR_3"        "MODISRAD:DSR_4"        "MODISRAD:Direct_1"     "MODISRAD:Direct_2"    
# [7] "MODISRAD:Direct_3"     "MODISRAD:Direct_4"     "MODISRAD:Diffuse_1"    "MODISRAD:Diffuse_2"    "MODISRAD:Diffuse_3"    "MODISRAD:Diffuse_4"   
#[13] "MODISRAD:GMT_0000_DSR" "MODISRAD:GMT_0300_DSR" "MODISRAD:GMT_0600_DSR" "MODISRAD:GMT_0900_DSR" "MODISRAD:GMT_1200_DSR" "MODISRAD:GMT_1500_DSR"
#[19] "MODISRAD:GMT_1800_DSR" "MODISRAD:GMT_2100_DSR" "MODISRAD:DSR_Quality" 

选择你想要的层，例如

s <- r[[4:11]] 

并计算平均值

x <- mean(s)
x
#class       : SpatRaster 
#dimensions  : 240, 240, 1  (nrow, ncol, nlyr)
#resolution  : 4633.127, 4633.127  (x, y)
#extent      : 0, 1111951, 4447802, 5559753  (xmin, xmax, ymin, ymax)
#coord. ref. : +proj=sinu +lon_0=0 +x_0=0 +y_0=0 +ellps=clrk66 +units=m +no_defs 
#source      : memory 
#name        :     mean 
#min value   : 10.92438 
#max value   : 70.07254 

保存到geotiff

writeRaster(x, filename="test.tif")

或一步完成

 m <- app(r[[4:11]], mean, filename="test.tif", overwrite=TRUE)

您还可以将文件作为子数据集的集合进行读取

 s <- sds(f)
 s
 #class       : SpatDataSet 
 #subdatasets : 12 
 #dimensions  : 240, 240 (nrow, ncol)
 #nlyr        : 4, 4, 4, 1, 1, 1, 1, 1, 1 
 #resolution  : 4633.127, 4633.127  (x, y)
 #extent      : 0, 1111951, 4447802, 5559753  (xmin, xmax, ymin, ymax)
 #coord. ref. : +proj=sinu +lon_0=0 +x_0=0 +y_0=0 +ellps=clrk66 +units=m +no_defs 
 #source(s)   : MCD18A1.A2020001.h18v04.006.2020050045103.hdf 
 #names       : DSR, Direct, Diffuse, GMT_0000_DSR, GMT_0300_DSR, GMT_0600_DSR, GMT_0900_DSR, GMT_1200_DSR, GMT_1500_DSR, GMT_1800_DSR, GMT_2100_DSR, DSR_Quality 

并像这样提取子数据集

s1 <- s[1]
s1
#class       : SpatRaster 
#dimensions  : 240, 240, 4  (nrow, ncol, nlyr)
#resolution  : 4633.127, 4633.127  (x, y)
#extent      : 0, 1111951, 4447802, 5559753  (xmin, xmax, ymin, ymax)
#coord. ref. : +proj=sinu +lon_0=0 +x_0=0 +y_0=0 +ellps=clrk66 +units=m +no_defs 
#data source : MCD18A1.A2020001.h18v04.006.2020050045103.hdf:MODISRAD:DSR 
#names       : MCD_1, MCD_2, MCD_3, MCD_4