1

我设法从 HDF 文件中提取 MODIS 土地覆盖数据并将其放入栅格中。

library(gdalUtils) #? gdal_translate?
library(raster)
sds <- get_subdatasets(".......myfileplath/file.hdf")
# Isolate the name of the first sds
name <- sds[5]
filename <- '2009test.tif'
gdal_translate(sds[5], dst_dataset = filename)
# Load the Geotiff created into R
r <- raster(filename)

我想把它放到一个数据框中，但是从原来的正弦曲线中重新投影+proj=sinu +lon_0=0 +x_0=0 +y_0=0 +a=6371007.181 +b=6371007.181 +units=m +no_defs，我认为。

对于与我正在分析的其他数据集兼容的普通椭圆体/WGS84。

这是我尝试过并且似乎有效的方法：

sr <- '+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0' 
 projected_raster <- projectRaster(r, crs = sr)

但是，当我随后将我的土地覆盖数据放入这个新投影中的数据框时，所有土地覆盖值都变为 NA。

这就是数据框在正弦投影中的样子（4 是土地覆盖分类）

DF <-  as.data.frame(r, xy=TRUE)
head(DF)

#           x       y       X2009test
#1   -1111718.9 1111719         4
#2   -1111255.6 1111719         4
#3   -1110792.2 1111719         4
#4   -1110328.9 1111719         4
#5   -1109865.6 1111719         4

通过我的重新投影，它看起来像这样：

reprojected_DF <-  as.data.frame(projected_raster, xy=TRUE)
head(reprojected_DF)

#                x        y X2009test
#1   -10.173076 10.01876        NA
#2   -10.168896 10.01876        NA
#3   -10.164716 10.01876        NA
#4   -10.160536 10.01876        NA
#5   -10.156356 10.01876        NA
#6   -10.152176 10.01876        NA
#7   -10.147996 10.01876        NA

关于我做错了什么或如何让土地覆盖坐标正确重新投影的任何建议？

干杯！！！！

我还读到有一个 NASA MODIS 重投影工具，但它不再存在/可用。有人知道吗？

4

3 回答 3

2

你所做的一切似乎都很好。除了在 projectRaster 中你应该使用method="ngb"，因为数字代表类（我想）。

你确定reprojected_DF没有值。你能吗show(reprojected_DF)（如果有的话，它应该显示最小值和最大值。也就是说，很可能这些值周围有一些 NA 值。在这种情况下，你可以这样做na.omit(reprojected_DF)。你也show(r)可以检查坐标参考系吗?

另一个可能更好的选择是将坐标投影到 中DF，像这样

DF <- data.frame(x=c(-1111718.9, -1111255.6, -1110792.2, -1110328.9, -1109865.6), y=c(1111719, 1111719, 1111719, 1111719, 1111719), X2009test=c(4, 4, 4, 4, 4))

library(rgdal)  
sincrs <- "+proj=sinu +lon_0=0 +x_0=0 +y_0=0 +a=6371007.181 +b=6371007.181 +units=m"
s <- SpatialPoints(DF, proj4string=CRS(sincrs))

lonlat <- '+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0' 

x <- spTransform(s, lonlat)
as.data.frame(x)
#          x        y X2009test
#1 -10.15209 9.997918         4
#2 -10.14786 9.997918         4
#3 -10.14362 9.997918         4
#4 -10.13939 9.997918         4
#5 -10.13516 9.997918         4

请注意，这与您显示的坐标位于同一区域 --- 再次表明您的数据中有实际值。

更简单的方法可能是使用 terra 包（与 raster 非常相似，但更简单、更快），然后执行 

library(terra)
r <- rast(".......myfileplath/file.hdf")

并使用几乎相同的代码从那里获取它（请参阅 ?terra 了解差异）

于 2020-05-17T04:19:17.860 回答
0

正如你所说，MRT不再可用。不过，您可以安装NASA's HDF-EOS to GeoTIFF Conversion Tool (HEG) （请参阅此处）来执行您需要的操作。我从来没有用过MRT，但我相信它是一回事。我已经使用HEG过几次，您可以批量运行（在基于 java 的 GUI 或命令行/终端中）重新投影（考虑nearest neighbour或bilinear插值）并将输出数据保存为*.tif.

于 2020-05-19T19:08:38.200 回答
0 
# Get a list of sds names
sds <- get_subdatasets(".......myfileplath/file.hdf")

sds

# Isolate the name of the first sds
name <- sds[5]
filename <- '2009test.tif'
gdal_translate(sds[5], dst_dataset = filename)
# Load the Geotiff created into R

r <- raster(filename)

# Define the Proj.4 spatial reference 

sr <- '+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0' 

projected_raster <- projectRaster(r, crs = sr)

DF <-  as.data.frame(r, xy=TRUE)
DF


reprojected_DF <-  as.data.frame(projected_raster, xy=TRUE)
reprojected_DF
于 2020-05-16T09:21:26.207 回答