I've downloaded a set of modis mod13q1 HDF files, e.g.
I'd like to load the
250m_16_days_pixel_reliability bands into R as rasters. My understanding is that I first need to use the modis reprojection tool (MRT) to create TIF files from the HDF. [Incorrect! The HDF file can be loaded directly into R; see answer and code below.] I'm using 64-bit linux and I have the MRT installed in
MRT gives several options for the output projection type. I know very little about projections. What parameters should I choose if my objectives are
- I'd like the conversion to run quickly
- I want to lose as little information as possible, i.e. I want to keep the 250 meter resolution of the original HDF file
- I need to be able to load the rasters into R
The ultimate goal is to extract raster values at a few thousand (say, 10,000) spatial points.
Edit: additional information, in case it's helpful:
- I'm getting my HDF data from https://search.earthdata.nasa.gov/search. As far as I can tell there is no option to download TIF files, or any other format that can be loaded directly into R
- I have two sets of HDF files for which I'd like to run this procedure. The first are in the USA (lower 48), and the second set covers Brazil.
Here's a simple example that seems to work. This borrows from Cannot properly import MODIS MOD07_L2 .hdf files into R using rgdal and mdsumner's answer. It skips the MRT step entirely.
library(raster) library(rgdal) hdf_path <- "/home/my_username/my_modis_dir" hdf_files <- list.files(hdf_path, pattern="^MOD13Q1.*\\.hdf$") stopifnot(length(hdf_files) > 0) infile <- hdf_files # MOD13Q1.A2000049.h11v04.006.2015136104609.hdf ## Modis h11v04 includes Wisconsin, parts of Michigan, parts of Illinois, ... gdal_info <- GDALinfo(infile, returnScaleOffset=FALSE) # Two warnings ## Two warnings: 1: In dim(x) : no bands in dataset; 2: In GDALinfo GeoTransform values not available subdataset_metadata <- attr(gdal_info, "subdsmdata") length(subdataset_metadata) # 24 strings starting with SUBDATASET_1_NAME=, SUBDATASET_1_DESC=, ... ## Extract substring following SUBDATASET_2_NAME= evi_subdataset <- subdataset_metadata[grep("250m 16 days EVI", subdataset_metadata)] evi_subdataset_name <- sapply(strsplit(evi_subdataset, "="), "[[", 2) # "HDF4_EOS:EOS_GRID: ... " spatial_grid_df <- readGDAL(evi_subdataset_name, as.is=TRUE) # Careful, need as.is=TRUE ## readGDAL call prints "... has GDAL driver HDF4Image and has 4800 rows and 4800 columns" str(spatial_grid_df) ## Formal class 'SpatialGridDataFrame' [package "sp"] with 4 slots ## @ data : 'data.frame': 23040000 obs. of 1 variable ## @ grid : Formal class 'GridTopology' [package "sp"] with 3 slots ## @ bbox : num [1:2, 1:2] -7783654 4447802 -6671703 5559753 ## @ proj4string:Formal class 'CRS' [package "sp"] with 1 slots ## projargs: chr "+proj=sinu +lon_0=0 +x_0=0 +y_0=0 +a=6371007.181 +b=6371007.181 +units=m +no_defs" summary(spatial_grid_df@data) ## band1 Min: -2000 Mean: 1249 Max: 9675 NAs:374180, plausible values for EVI ## See https://lpdaac.usgs.gov/dataset_discovery/modis/modis_products_table/mod13q1_v006 evi_raster <- raster(spatial_grid_df) identical(crs(evi_raster), crs(spatial_grid_df)) # True plot(evi_raster) # Extent [-7783654, -6671703] for x, [4447802, 5559753] for y ## Sanity check: extract EVI values for points in Lake Michigan n_points <- 1000 lat <- seq(42.5, 44, length.out=n_points) lon <- -runif(n_points, min=86.6, max=87.4) wgs84 <- "+init=epsg:4326 +proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs +towgs84=0,0,0" points_in_lake_michigan <- SpatialPoints(coords=cbind(lon, lat), proj4string=CRS(wgs84)) evi_in_lake_michigan <- extract(evi_raster, spTransform(points_in_lake_michigan, crs(evi_raster))) summary(evi_in_lake_michigan) # 21 NAs, otherwise median around -200, plausible for water ## Points in Wisconsin lat <- seq(44, 45, length.out=n_points) lon <- -runif(n_points, min=89, max=90) points_in_wisconsin <- SpatialPoints(coords=cbind(lon, lat), proj4string=CRS(wgs84)) evi_in_wisconsin <- extract(evi_raster, points_in_wisconsin) summary(evi_in_wisconsin) # Zero NAs, median around 1750 plot(evi_raster) points(spTransform(points_in_lake_michigan, crs(evi_raster)), col="red") points(spTransform(points_in_wisconsin, crs(evi_raster)), col="blue", pch=2) ## Plot above looks correct (but gets messed up if I resize the window) usa <- getData("GADM",country="USA",level=1) plot(subset(usa, NAME_1 %in% c("Wisconsin", "Illinois", "Indiana", "Michigan")), border="black") points(points_in_lake_michigan, col="red") points(points_in_wisconsin, col="blue", pch=2) # Looks correct