MacOS 上 R 中的 MODIStsp：gdal_utils(“buildvrt”、gdalfile、output.vrt、opts) 中出现错误：gdal_utils buildvrt：发生错误

Question

我的项目是关于加利福尼亚州的羊，我试图获取该地区随时间变化的植被测量值，看看疾病流行率和植被增加（以 NDVI 测量）之间是否存在任何相关性。为此，我找到了一个介绍如何使用 MODIStsp 的网站： https://flograttarola.com/ post/modis-downloads/。我有疾病存在的坐标，但我需要 NDVI 数据。

我正在尝试获取加利福尼亚州的 NDVI 数据（从此处下载的 shapefile：https:// data.ca.gov/dataset/ca-geographic-boundaries）使用 MODIStsp 包导入 R 并输入以下条件：

#并非所有这些包都是必需的，但我已经包含了我尝试过的包。必要的代码有星号（**）突出显示：

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

但我一直遇到同样的错误：

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

如果有人能帮助我解决这个问题，我将非常感激！

Answer 1

使用下面的示例数据就可以了。

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