MODISGrid.R 3.8 KB

1234567891011121314151617181920212223242526272829303132333435363738394041424344454647484950515253545556575859606162636465666768697071727374757677
  1. MODISGrid <-
  2. function(Dir = ".", DirName = "MODIS_GRID", SubDir = TRUE, NoDataValues)
  3. {
  4. ## DEFINE
  5. NUM_METADATA_COLS <- 10
  6. if(Dir == '.') cat('Files downloaded will be written to ', file.path(getwd(), DirName), '.\n', sep = '')
  7. if(Dir != '.') cat('Files downloaded will be written to ', file.path(Dir, DirName), '.\n', sep = '')
  8. ## Create directory for storing new grid files.
  9. if(!file.exists(file.path(Dir, DirName))) dir.create(path = file.path(Dir, DirName))
  10. ## Find all MODIS data files in Dir.
  11. file.list <- list.files(path = Dir, pattern = ".asc$")
  12. file.list <- file.list[grepl("___", file.list)]
  13. if(length(file.list) == 0) stop("Could not find any MODIS data files in directory specified.")
  14. ## Check NoDataValues is a list of named vectors.
  15. if(!is.list(NoDataValues)) stop("NoDataValues should be a list of named vectors. See help doc.")
  16. ## Check the number of products in NoDataValues list equals the number found in file.list.
  17. prod.set <- unique(do.call(c, as.vector(lapply(
  18. file.path(Dir, file.list), function(x) read.csv(x, header = FALSE, as.is = TRUE)[ ,7]
  19. ))))
  20. if(any(nchar(prod.set) == 0)) stop("A subset was incompletely downloaded. Check MODISSubsets output and retry the subset.")
  21. if(!all(prod.set %in% names(NoDataValues))) stop("Mismatch between NoDataValues and data products found in files.")
  22. ## Check that NoDataValues value is specified for every data band found in file.list.
  23. band.set <- unique(as.vector(sapply(
  24. lapply(file.path(Dir, file.list), function(x) read.csv(x, header = FALSE, as.is = TRUE)[ ,6]),
  25. function(x) unique(substr(x, (gregexpr(".", x, fixed = TRUE)[[1]][5] + 1), nchar(x)))
  26. )))
  27. if(!all(band.set %in% names(do.call(c, unname(NoDataValues))))){
  28. stop("Mismatch between NoDataValues and data bands found in files.")
  29. }
  30. for(i in 1:length(file.list))
  31. {
  32. cat("Creating new GIS ASCII files from MODIS data file", i, "out of", length(file.list), "\n")
  33. data.file <- read.csv(file.path(Dir, file.list[i]), header = FALSE, as.is = TRUE)
  34. names(data.file) <- c("nrow", "ncol", "xll", "yll", "pixelsize", "row.id", "product.code", "MODIS.acq.date",
  35. "where", "MODIS.proc.date", 1:(ncol(data.file) - NUM_METADATA_COLS))
  36. ## Create directory for this data file if SubDir = TRUE.
  37. sub.dir <- substr(file.list[i], 1, regexpr(".asc$", file.list[i])-1)
  38. if(SubDir & !file.exists(file.path(Dir, DirName, sub.dir))) dir.create(path = file.path(Dir, DirName, sub.dir))
  39. for(n in 1:nrow(data.file))
  40. {
  41. data.band <- substr(data.file$row.id[n],
  42. gregexpr(".", data.file$row.id[n], fixed = TRUE)[[1]][5] + 1,
  43. nchar(data.file$row.id[n]))
  44. data.date <- data.file$MODIS.acq.date[n]
  45. path <- ifelse(SubDir,
  46. file.path(Dir, DirName, sub.dir,
  47. paste0("GRID_", sub.dir, "_", data.band, "_", data.date)),
  48. file.path(Dir, DirName,
  49. paste0("GRID_", sub.dir, "_", data.band, "_", data.date)))
  50. write(c(sprintf("ncols\t\t %i", data.file$ncol[n]),
  51. sprintf("nrows\t\t %i", data.file$nrow[n]),
  52. sprintf("xllcorner\t %.2f", data.file$xll[n]),
  53. sprintf("yllcorner\t %.2f", data.file$yll[n]),
  54. sprintf("cellsize\t %s", as.character(data.file$pixelsize[n])),
  55. sprintf("NODATA_value\t %s", as.character(NoDataValues[[data.file$product.code[n]]][data.band]))),
  56. file = file.path(paste0(path,".asc")))
  57. WritePRJ(Path = file.path(paste0(path,".prj")))
  58. grid.data <- matrix(data.file[n,(NUM_METADATA_COLS+1):ncol(data.file)],
  59. nrow = data.file$nrow[n], ncol = data.file$ncol[n], byrow = TRUE)
  60. write.table(grid.data, file = file.path(paste0(path,".asc")), append = TRUE, col.names = FALSE, row.names = FALSE)
  61. }
  62. }
  63. }