In R, how do I reproject from the MODIS sinusoidal projection to latlong (ellps=WGS84) projection

Question

I have managed to extract MODIS land cover data from an HDF file and put it into a raster.

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)

I want to put it into a dataframe but reproject it from the original sinusoidal which is +proj=sinu +lon_0=0 +x_0=0 +y_0=0 +a=6371007.181 +b=6371007.181 +units=m +no_defs, I think.

To a normal ellipsoid/WGS84 one that is compatible with other datasets I'm analyzing.

This is what I've tried and seemed to work:

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

However, when I then put my my landcover data in this new projection into a data frame all the landcover values have gone to NA.

This is what the dataframe looked like in a sinu projection (the 4's being a land cover classifications)

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

and with my reprojection it looks like this:

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

Any advice on what I'm doing wrong or how I can get the land cover coordinates to reproject properly?

Cheers!!!!

Also I read there was a NASA MODIS reprojection tool but that that no long exists / is available. Anyone know anything about that?

Answer 1

# 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

Answer 2

All you are doing seems to be OK. Except that in projectRaster you should use method="ngb" because the numbers represent classes (I suppose).

Are you sure reprojected_DF has no values. Can you show(reprojected_DF) (it should show the min and max values if there are any. That is, it could very well be that there are some NA values around the values. In that case you may do na.omit(reprojected_DF). Can you also show(r) to check the coordinate reference system?

Another, possibly better, option, would be to project the coordinates in DF, like this

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

Note that this is in the same area as the coordinates you showed --- again suggesting that there are in fact values in your data.

An easier route might be to use the terra package (very similar to raster, but simpler and faster), and do

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

And take it from there with pretty much the same code (see ?terra for the differences)

Answer 3

As you said, MRT is no longer available. Nevertheless, you can install NASA's HDF-EOS to GeoTIFF Conversion Tool (HEG) (see here) to do what you need. I never used MRT, but I believe it's kind of the same thing. I've used HEG several times, and you can run in batch (in a java-based GUI or in command line/terminal) reprojection (considering nearest neighbour or bilinear interpolation) and save output data as *.tif.

Answer 4

library(terra)
f<-list.files(pattern="hdf")
s <- sds(f[1])
r<-s[1]
r<-project(r, "EPSG:4326")

The terra package is definitely the choice. Quicker and Safer.