我设法从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。
这是数据框在sinu投影中的样子（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
有什么建议可以告诉我，我做错了什么 或者我怎么才能让土地覆盖坐标正确地重投？
Cheers!!!
此外，我读到有一个美国宇航局MODIS重投工具，但这早已不存在是可用的。任何人都知道任何有关？
你所做的一切似乎都没问题。除了在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包（非常类似于栅格，但更简单和更快），然后做
library(terra)
r <- rast(".......myfileplath/file.hdf")
然后用几乎相同的代码从那里开始(不同之处见?terra)
# 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
正如你所说 MRT 已不再可用。不过，您可以安装 NASA's HDF-EOS to GeoTIFF Conversion Tool (HEG) (见此)来做你需要的事情。我从来没有用过 MRT但我相信这是同样的事情。我已经用 HEG 多次，你可以批量运行（在一个基于java的GUI或在命令行中）重投（考虑到 nearest neighbour 或 bilinear 插值）并将输出数据保存为 *.tif.