Quantcast

MODIS package and gdalUtils

Previous Topic Next Topic
 
classic Classic list List threaded Threaded
1 message Options
Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

MODIS package and gdalUtils

Vladimir Ruslan Wingate
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]]

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Loading...