[R-sig-Geo] MODIS package and gdalUtils
Vladimir Ruslan Wingate
vladimir.wingate at unibas.ch
Thu Apr 13 15:23:54 CEST 2017
Hi all,
I get the following error from running to following code using the MODIS package and gdalUtils. Thanks.
names(s) <- as.Date(modis.hdf[[1]], "A%
Error in `names<-`(`*tmp*`, value = c(NA_real_, NA_real_)) :
incorrect number of layer names
Here is the whole code:
library(MODIS)
library(gdalUtils)
MODISoptions(gdalPath="C:/OSGeo4W64/bin")
#install.packages("MODIS", repos ="http://R-Forge.R-project.org")
setwd("J:/chapter 11")
getwd()
gdal_setInstallation(search_path = "C:/OSGeo4W64/bin", rescan = TRUE,ignore.full_scan = TRUE, verbose = FALSE)
MODISoptions(localArcPath = getwd(), outDirPath =getwd())
MODIS:::checkTools('GDAL')
MODIS:::checkTools('MRT')
viname <- "ndvi"
product <- "MOD13Q1"
ofilename <- paste0(product, "_", viname, "_brick.grd")
ofilename
pth <- paste0(getwd(), "/raster_data/", product)
pth
fileout <- paste(pth, "/", ofilename, sep="")
fileout
if (!file.exists(fileout)) {
if (!file.exists(pth)) {
print("the outfolder does not exist and will be created")
print(pth)
dir.create(pth, recursive = TRUE)
}
}
tileH <- 19
tileV <- 10
begin <- "2015.01.01"
end <- "2015.01.20"
modis.hdf <- getHdf(product = product, begin=begin, end=end, tileH=tileH, tileV=tileV, checkIntegrity=TRUE)
modis.hdf
for (i in 1:length(modis.hdf[[1]])) {
ifl <- unlist(strsplit(unlist(strsplit(modis.hdf[[1]][i], c("[/]"))) [6], "[.]")) [1]#5
ifl
print(ifl)
fn <- paste("cvi_", ifl, ".tif", sep="")
print(fn)
if (is.na(ifl) | file.exists(fn)) {
print("file exists or is not available on the server")
} else {
sds <- get_subdatasets(modis.hdf[[1]][i])
}
}
tmp <- rasterTmpFile()
extension(tmp) <- "tif"
library(gdalUtils)
gdal_translate(sds[1], dst_dataset=tmp)
ndvi <- raster(tmp)/10000
writeRaster(ndvi, filename=paste0("ndvi_", ifl, ".tif", sep=""), dataType="INT2U", format="GTiff", overwrite=TRUE)
tmp2 <- rasterTmpFile()
extension(tmp2) <- "tif"
gdal_translate(sds[12], dst_dataset=tmp2)
rel <- crop(x=raster(tmp2), y=raster(tmp2),
filename=paste("rel_", ifl, ".tif", sep=""),
dataType="INT2U", format="GTiff")
f_cleanvi <- function(vi, rel) {
i <- (rel <= 1)
res <- matrix(NA, length(i), 1)
if (sum(i, na.rm=TRUE) > 0) {
i <- which(i)
res[i] <- vi[i]
}
res
}
clvi <- overlay(ndvi, rel, fun=f_cleanvi, filename=fn, dataType="INT2U", overwrite=TRUE)
clvi
rm(ndvi, rel, clvi, tmp, tmp2)
ofl <- list.files(pattern=glob2rx("rel*.tif"))
class(ofl)
ofl
modis.hdf[[1]]
s <- do.call("brick", lapply(ofl, raster))
names(s) <- as.Date(modis.hdf[[1]], "A%Y%j")
[[alternative HTML version deleted]]
More information about the R-sig-Geo
mailing list