将 HDF5 转换为光栅格式的问题

Question

我正在尝试将我的 h5 每日文件转换为光栅格式。我转换成光栅格式。当我提取我感兴趣的领域时。我无法从光栅图像中提取我感兴趣的区域请任何人指导我如何解决这个问题。 R 代码和 hf5 文件以及转换后的光栅图像存在于链接中（附加）。谢谢

library(rhdf5)
library(sp)
library(raster)
h5ls("reconstruction_indus_CY2001.h5")
h5readAttributes(file = "reconstruction_indus_CY2001.h5", name = "Grid")
h5f = H5Fopen("reconstruction_indus_CY2001.h5")
# h5f
# h5f&'Grid'

#system.time( swe <- h5f$Grid$swe )
system.time( melt <- h5f$Grid$melt )

locations <- data.frame(
  lon=c(74.86764,73.48753, 74.87066 , 73.37798 , 78.82102 ,75.85160 ,75.78263 , 78.46446 ), 
  lat = c(35.16700, 36.25674, 36.49362, 35.21188, 34.20916, 34.48459, 35.76965, 33.23380)
)

coordinates(locations) <- ~lon+lat
proj4string(locations) <-  CRS("+proj=longlat")

  swe180 <- melt[,,180]
  b <- swe180 == 65535
  # table(b)
  swe180[b] <- -1
  
  b <- swe180 > 200
  # table(b)
  swe180[b] <- 200
  
  b <- swe180 < 0
  # table(b)
  swe180[b] <- 20
  
  # image(swe180)
  
  # image(swe180)
  # str(swe180)
  
  # h5readAttributes(file = "reconstruction_Sierra_2016.h5", name = "Grid")$ReferencingMatrix
  
  RM <- h5readAttributes(file = "reconstruction_indus_CY2001.h5", name = "Grid")$ReferencingMatrix
  #GT <- GridTopology(c(RM[3,1], RM[3,2]+RM[1,2]*dim(swe)[1]), c(RM[2,1], -RM[1,2]), c(dim(swe)[2],dim(swe)[1]))
  GT <- GridTopology(c(RM[3,1], RM[3,2]+RM[1,2]*dim(melt)[1]), c(RM[2,1], -RM[1,2]), c(dim(melt)[2],dim(melt)[1]))
  # GT <- GridTopology(c(-1.088854e+07, 4718608.3619-463.3127*1978), c(463.3127, 463.3127), c(2171,1978))
  # GT
  SG = SpatialGrid(GT)
  # str(SG)
  # proj4string(SG) <- CRS("+proj=sinu")
  # str(SG)
  proj4string(SG) <- CRS("+proj=utm +zone=43 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0")
  locations_aea <- spTransform(locations, CRS(proj4string(SG)))
  SGDF = SpatialGridDataFrame(SG, data.frame(melt = as.numeric(t(swe180))))
  gridded(SGDF)<- TRUE
  r = raster(SGDF)
  plot(SGDF, axes=T)
  writeRaster(r,"test_2001.tif",overwrite=TRUE)
  
## Open Raster Files and Extract Area of Interest

shp= readOGR("Hunza.shp")
e = extent(shp)
   r1  = raster("test_2001.tif")
  crs(r1) = "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 "
  plot(r1)
  r1_mask = raster::mask(r1,shp)
  plot(r1_mask,axes = TRUE,ext = extent(shp))
  # Extracting Values as Data Frame
  r1_extract = raster::extract(r1,shp, df=TRUE,na.rm = TRUE)
  # Stroing as Raster
  writeRaster(r1_mask,paste0('/shared/MODIS/shastaH5SWEinR/2001_swe/Hunza/','hunza.tif'))
  c = cbind(r1_extract,y)
  c1=t(c)
  write.csv(c1,file = 'Hunza_SWE_2001.csv')

https://drive.google.com/drive/folders/18-hj2LEYWBN-uIDDTdqZ-x-WUxpCJu7H?usp=sharing

Answer 1

您可以使用terra 包（raster 替换）来快捷方式。 terra 可以直接读取 hdf5 文件。

该文件有多个子数据集，因此最容易将其读取为 SpatDataSet

library(terra)
f <- "reconstruction_indus_CY2001.h5"
s <- sds(f)
s
#class       : SpatDataSet 
#subdatasets : 3 
#dimensions  : 3651, 1641 (nrow, ncol)
#nlyr        : 1, 365, 365 
#resolution  : 0.2193784, 0.04930156  (x, y)
#extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
#coord. ref. : +proj=longlat +datum=WGS84 +no_defs 
#names       : maxswedates, melt, swe

现在得到感兴趣的变量

r <- s$swe
r
#class       : SpatRaster 
#dimensions  : 3651, 1641, 365  (nrow, ncol, nlyr)
#resolution  : 0.2193784, 0.04930156  (x, y)
#extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
#coord. ref. : +proj=longlat +datum=WGS84 +no_defs 
#data source : swe 
#names       : swe_1, swe_2, swe_3, swe_4, swe_5, swe_6, ... 

获得相同结果的更直接的方法是

r <- rast(f, "//Grid/swe")

您可以通过运行发现 HDF5 文件的内容

sds_info(f)

绘制第一层

plot(r, 1)

提取感兴趣的区域，例如这样

v <- vect("Hunza.shp")
x <- crop(r, v)
y <- mask(x, v)

要保存为光栅文件，您可以为上述函数添加文件名。或者你可以稍后再这样做

y <- writeRaster(y, "hunza.tif")

要将值保存到 csv 文件：

vy <- values(y)       
write.csv(vy, 'Hunza_SWE_2001.csv', row.names=FALSE)

terra 中的大多数函数名称与 raster 中的相同。有关差异，请参阅?terra。如果你想继续raster，你可以这样做

library(raster)
b <- brick(y)