MODIStsp in R on MacOS: Error in gdal_utils("buildvrt", gdalfile, output.vrt, opts) : gdal_utils buildvrt: an error occured

Question

My project is on sheep in California and I'm trying to get the vegetation measures of the region over time to see if there's any correlation between disease prevalence and increased vegitation (measured as NDVI). To do this I found a website that went through how to use MODIStsp: https://flograttarola.com/post/modis-downloads/. I have coordinates for the presence of disease but I need the NDVI data.

I am trying to get the NDVI data for California (shapefile downloaded from here: https://data.ca.gov/dataset/ca-geographic-boundaries) into R using the MODIStsp package and inputting my conditions below:

#Not all of these packages are necessary but I have included those that I have tried. The necessary code has stars (**) highlighting it:

**install.packages("MODIStsp")**
**library(MODIStsp)**
**install.packages('terra')**
**library("terra")**
**library(sf)**
**library(tidyverse)**
install.packages("gdalUtils")
library(gdalUtils)
install.packages('rgdal', type='source')
library(rgdal)
library(raster)
library(here)
library(ggplot2)
install.packages("hdf4r")
library("hdf4r")

**spatial_file <-("~/Downloads/ca-state-boundary/CA_State_TIGER2016.shp")**

**MODIStsp(gui             = FALSE,
         out_folder      = "Downloads", #Where I want to store my data
         out_folder_mod  = "Downloads",
         selprod         = 'Vegetation Indexes_16Days_250m (M*D13Q1)',
         sensor          = 'Terra',
         bandsel         = 'NDVI', #Get NDVI data
         spatmeth        = 'file',
         spafile         = spatial_file, #Spatial file of California
         indexes_bandsel = "SR", 
         user            = 'EXAMPLE', # your username for NASA http server
         password        = 'EXAMPLE',  # your password for NASA http server
         start_date      = '2000.01.01', 
         end_date        = '2022.12.31', 
         verbose         = TRUE,
         out_format      = 'GTiff',
         compress        = 'LZW',
         out_projsel     = 'User Defined',
         output_proj     = '+proj=latlong +datum=WGS84 +units=m +no_defs', #I want to get the NDVI data for coordinates
         delete_hdf      = TRUE, #Delete hdf files after making them into GTiff
         parallel        = TRUE
)**

But I keep running into the same error:

<
Error in gdal_utils("buildvrt", gdalfile, output.vrt, opts) : 
  gdal_utils buildvrt: an error occured
In addition: Warning messages:
1: In CPL_gdalbuildvrt(source, destination, options, oo) :
  GDAL Error 4: `Downloads/MOD13Q1.A2000049.h08v04.006.2015136104517.hdf' not recognized as a supported file format.
2: In CPL_gdalbuildvrt(source, destination, options, oo) :
  GDAL Message 1: Can't open Downloads/MOD13Q1.A2000049.h08v04.006.2015136104517.hdf. Skipping it
3: In CPL_gdalbuildvrt(source, destination, options, oo) :
  GDAL Error 4: `Downloads/MOD13Q1.A2000049.h09v04.006.2015136104603.hdf' not recognized as a supported file format.
4: In CPL_gdalbuildvrt(source, destination, options, oo) :
  GDAL Message 1: Can't open Downloads/MOD13Q1.A2000049.h09v04.006.2015136104603.hdf. Skipping it
5: In CPL_gdalbuildvrt(source, destination, options, oo) :
  GDAL Error 4: `Downloads/MOD13Q1.A2000049.h07v05.006.2015136104501.hdf' not recognized as a supported file format.
6: In CPL_gdalbuildvrt(source, destination, options, oo) :
  GDAL Message 1: Can't open Downloads/MOD13Q1.A2000049.h07v05.006.2015136104501.hdf. Skipping it
7: In CPL_gdalbuildvrt(source, destination, options, oo) :
  GDAL Error 4: `Downloads/MOD13Q1.A2000049.h08v05.006.2015136104621.hdf' not recognized as a supported file format.
8: In CPL_gdalbuildvrt(source, destination, options, oo) :
  GDAL Message 1: Can't open Downloads/MOD13Q1.A2000049.h08v05.006.2015136104621.hdf. Skipping it
9: In CPL_gdalbuildvrt(source, destination, options, oo) :
  GDAL Error 4: `Downloads/MOD13Q1.A2000049.h09v05.006.2015136104623.hdf' not recognized as a supported file format.
10: In CPL_gdalbuildvrt(source, destination, options, oo) :
  GDAL Message 1: Can't open Downloads/MOD13Q1.A2000049.h09v05.006.2015136104623.hdf. Skipping it
>

If anyone could help me solve this, I'd be very grateful!

Welcome to Stackoverflow. What does `sf::st_drivers('raster')` say about your HDF4/HDF5 capabilities? It might indicate needing to rebuild gdal in the presence of `libdf` and `libhdf5` if you don't see them. — Chris, Mar 16 '22 at 14:58
It says Name: HDF5, Write:FALSE, Copy:FALSE, is.raster:TRUE , is.vector: FALSE , and vsi: TRUE but nothing on HDF4 — Falloutalley, Mar 17 '22 at 15:10
While not conclusive, and not knowing how to check for presence of `libdf` nor not knowing if the hdf of interest is 4 or 5, gdal, from source, generously builds every driver it can based upon libs/headers discovered in the configure process. [HDF EOS Tools R](https://hdfeos.org/software/r.php) references to mac installs as well, guessing there's a `home brew` session in your near future. — Chris, Mar 17 '22 at 16:08

Christoph · Accepted Answer · 2022-03-18T16:13:30.783

With your sample data below it works.

library(sf)
# Linking to GEOS 3.9.1, GDAL 3.2.1, PROJ 7.2.1; sf_use_s2() is TRUE
library(leaflet)

path <- "C:/temp/"
setwd(path)
filename <- "MOJA_boundary.shp"

y4 <- st_read(dsn = filename)
# Reading layer `MOJA_boundary' from data source `C:\temp\MOJA_boundary.shp' using driver `ESRI Shapefile'
# Simple feature collection with 1 feature and 19 fields
# Geometry type: MULTIPOLYGON
# Dimension:     XY
# Bounding box:  xmin: -116.165 ymin: 34.71693 xmax: -114.9492 ymax: 35.59077
# Geodetic CRS:  NAD83
class(y4)
# [1] "sf"         "data.frame"
st_geometry_type(y4)
# [1] MULTIPOLYGON
# 18 Levels: GEOMETRY POINT LINESTRING POLYGON MULTIPOINT MULTILINESTRING MULTIPOLYGON ... TRIANGLE
st_crs(y4)
# Coordinate Reference System:
#   User input: NAD83 
# wkt:
#   GEOGCRS["NAD83",
#           DATUM["North American Datum 1983",
#                 ELLIPSOID["GRS 1980",6378137,298.257222101,
#                           LENGTHUNIT["metre",1]]],
#           PRIMEM["Greenwich",0,
#                  ANGLEUNIT["degree",0.0174532925199433]],
#           CS[ellipsoidal,2],
#           AXIS["latitude",north,
#                ORDER[1],
#                ANGLEUNIT["degree",0.0174532925199433]],
#           AXIS["longitude",east,
#                ORDER[2],
#                ANGLEUNIT["degree",0.0174532925199433]],
#           ID["EPSG",4269]]
st_bbox(y4)
# xmin       ymin       xmax       ymax 
# -116.16503   34.71693 -114.94916   35.59077 

y6 <- st_transform(y4, "+proj=longlat +datum=WGS84 +no_defs")
st_crs(y6)
# Coordinate Reference System:
#   User input: +proj=longlat +datum=WGS84 +no_defs 
# wkt:
#   GEOGCRS["unknown",
#           DATUM["World Geodetic System 1984",
#                 ELLIPSOID["WGS 84",6378137,298.257223563,
#                           LENGTHUNIT["metre",1]],
#                 ID["EPSG",6326]],
#           PRIMEM["Greenwich",0,
#                  ANGLEUNIT["degree",0.0174532925199433],
#                  ID["EPSG",8901]],
#           CS[ellipsoidal,2],
#           AXIS["longitude",east,
#                ORDER[1],
#                ANGLEUNIT["degree",0.0174532925199433,
#                          ID["EPSG",9122]]],
#           AXIS["latitude",north,
#                ORDER[2],
#                ANGLEUNIT["degree",0.0174532925199433,
#                          ID["EPSG",9122]]]]
plot(y6)
# Warnmeldung:
#   plotting the first 9 out of 20 attributes; use max.plot = 19 to plot all 

y6$NAME
# NULL
y6$geometry
# Geometry set for 1 feature 
# Geometry type: MULTIPOLYGON
# Dimension:     XY
# Bounding box:  xmin: -116.165 ymin: 34.71693 xmax: -114.9492 ymax: 35.59077
# CRS:           +proj=longlat +datum=WGS84 +no_defs
# MULTIPOLYGON (((-114.9529 35.15707, -114.9533 3...
plot(y6$geometry)
y6_2 <- st_transform(y6$geometry, "+init=epsg:4326")
# Warnmeldung:
#   In CPL_crs_from_input(x) :
#   GDAL Message 1: +init=epsg:XXXX syntax is deprecated. It might return a CRS with a non-EPSG compliant axis order.
plot(y6_2)
leaflet() %>% addTiles() %>% 
  addPolygons(data = y6)

y6_2.xy <- st_zm(y6_2) 
y6_2.sp <- sf:::as_Spatial(y6_2.xy) 

leaflet() %>% addTiles() %>% 
  addPolygons(data = y6_2.sp@polygons[[1]])

Hi there, thanks so much for your help! It works up until line '''plot(y6_TG12)''', when it produces an error: Error in plot_sf(x, ...) : NA value(s) in bounding box. Trying to plot empty geometries? — Falloutalley, Mar 17 '22 at 15:21
I also have an alternative shapefile of just Mojave desert boundaries I could use instead of Californiastate boundaries, available at https://catalog.data.gov/dataset/mojave-national-preserve-tract-and-boundary-data/resource/ae577a6e-4948-405c-8ccc-364264676115 — Falloutalley, Mar 17 '22 at 15:25
@Falloutalley See my edit. It works... It turned out, that "your" names are NULL in the example case. But you should be able to adapt this... — Christoph, Mar 18 '22 at 16:14

MODIStsp in R on MacOS: Error in gdal_utils("buildvrt", gdalfile, output.vrt, opts) : gdal_utils buildvrt: an error occured

1 Answers1