MacOS 上 R 中的 MODIStsp：gdal_utils("buildvrt"、gdalfile、output.vrt、opts) 中的错误：gdal_utils buildvrt：发生错误
MODIStsp in R on MacOS: Error in gdal_utils("buildvrt", gdalfile, output.vrt, opts) : gdal_utils buildvrt: an error occured
我的项目是关于加利福尼亚州的绵羊，我正在尝试获取该地区随时间推移的植被测量值，以查看疾病流行与植被增加（以 NDVI 衡量）之间是否存在任何相关性。为此，我找到了一个介绍如何使用 MODIStsp 的网站：https://flograttarola.com/post/modis-downloads/。我有疾病存在的坐标，但我需要 NDVI 数据。
我正在尝试使用 MODIStsp 包将加利福尼亚的 NDVI 数据（从此处下载的 shapefile：https://data.ca.gov/dataset/ca-geographic-boundaries）导入 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
>
如果有人能帮我解决这个问题，我将不胜感激！
下面的示例数据有效。
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]])
我的项目是关于加利福尼亚州的绵羊，我正在尝试获取该地区随时间推移的植被测量值，以查看疾病流行与植被增加（以 NDVI 衡量）之间是否存在任何相关性。为此，我找到了一个介绍如何使用 MODIStsp 的网站：https://flograttarola.com/post/modis-downloads/。我有疾病存在的坐标，但我需要 NDVI 数据。
我正在尝试使用 MODIStsp 包将加利福尼亚的 NDVI 数据（从此处下载的 shapefile：https://data.ca.gov/dataset/ca-geographic-boundaries）导入 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
>
如果有人能帮我解决这个问题，我将不胜感激！
下面的示例数据有效。
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]])