rmacosgdalr-modis
MODIStsp in R on MacOS: Error in gdal_utils("buildvrt", gdalfile, output.vrt, opts) : gdal_utils buildvrt: an error occured
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!
Solution
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]])
