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
Get help from others!
Recent Questions
Recent Answers
© 2024 TransWikia.com. All rights reserved. Sites we Love: PCI Database, UKBizDB, Menu Kuliner, Sharing RPP