TransWikia.com

Different results extracting data from the same point from MODIS HDF with R

Geographic Information Systems Asked on December 26, 2020

I am trying to extract NDVI data from MODIS HDF files using R. I have several files downloaded from EarthData.

The file I am using in this example is here to download.

These MODIS HDF files have a sinusoidal projection, so I have to reproject the coordinates of my points of interest from WGS84 to the projection of the .hdf file or viceversa. Trying both I get different results for the same point. How can I know which one is correct?

Here is my code:

library(sp) # SpatialPoints(), CRS()
library(MODIS)
library(rgdal) # readGDAL()

# defining the coordinates of interest
site = cbind(-71, -41)
latlon = CRS('+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0')
site_wgs84 = SpatialPoints(site, latlon)

# I have several files but here I open just one as an example
f <- "MOD13Q1.A2000049.h12v13.006.2015136104516.hdf"

# open the file
sds <- getSds(f)

# in this case I am interested in the "250m 16 days NDVI" dataset
ndvi_sinus <- raster(readGDAL(sds$SDS4gdal[1], as.is=TRUE))

# transform the site coordinates to the sinusoidal projection
sinus <- crs(ndvi_sinus)
site_sinus <- spTransform(site_wgs84, sinus)

# extract the data divided by the scale_factor
data_ndvi_sinus <- extract(ndvi_sinus, site_sinus) / 10000 # I get data_ndvi_sinus=0.1912

# Now I reproject the raster to WGS84 and then extract the data
ndvi_wgs84 <- projectRaster(ndvi_sinus, crs=latlon)
data_ndvi_wgs84 <- extract(ndvi_wgs84, site_wgs84) / 10000 # I get data_ndvi_wgs84=0.2162236

Add your own answers!

Ask a Question

Get help from others!

© 2024 TransWikia.com. All rights reserved. Sites we Love: PCI Database, UKBizDB, Menu Kuliner, Sharing RPP