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!

Answer 1

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]])