LandCover <- function(Dir = ".", Band) { ########## Define land cover classes for each lc band. LC_CLASS <- list( Land_Cover_Type_1 = c("Water" = 0, "Evergreen Needleleaf forest" = 1, "Evergreen Broadleaf forest" = 2, "Deciduous Needleleaf forest" = 3, "Deciduous Broadleaf forest" = 4, "Mixed forest" = 5, "Closed shrublands" = 6, "Open shrublands" = 7, "Woody savannas" = 8, "Savannas" = 9, "Grasslands" = 10, "Permanent wetlands" = 11, "Croplands" = 12, "Urban & built-up" = 13, "Cropland/Natural vegetation mosaic" = 14, "Snow & ice" = 15, "Barren/Sparsely vegetated" = 16, "Unclassified" = 254, "NoDataFill" = 255), Land_Cover_Type_2 = c("Water" = 0, "Evergreen Needleleaf forest" = 1, "Evergreen Broadleaf forest" = 2, "Deciduous Needleleaf forest" = 3, "Deciduous Broadleaf forest" = 4, "Mixed forest" = 5, "Closed shrublands" = 6, "Open shrublands" = 7, "Woody savannas" = 8, "Savannas" = 9, "Grasslands" = 10, "Croplands" = 12, "Urban & built-up" = 13, "Barren/Sparsely vegetated" = 16, "Unclassified" = 254, "NoDataFill" = 255), Land_Cover_Type_3 = c("Water" = 0, "Grasses/Cereal crops" = 1, "Shrubs" = 2, "Broadleaf crops" = 3, "Savanna" = 4, "Evergreen Broadleaf forest" = 5, "Deciduous Broadleaf forest" = 6, "Evergreen Needleleaf forest" = 7, "Deciduous Needleleaf forest" = 8, "Non-vegetated" = 9, "Urban" = 10, "Unclassified" = 254, "NoDataFill" = 255), Land_Cover_Type_4 = c("Water" = 0, "Evergreen Needleleaf forest" = 1, "Evergreen Broadleaf forest" = 2, "Deciduous Needleleaf forest" = 3, "Deciduous Broadleaf forest" = 4, "Annual Broadleaf vegetation" = 5, "Annual grass vegetation" = 6, "Non-vegetated land" = 7, "Urban" = 8, "Unclassified" = 254, "NoDataFill" = 255), Land_Cover_Type_5 = c("Water" = 0, "Evergreen Needleleaf forest" = 1, "Evergreen Broadleaf forest" = 2, "Deciduous Needleleaf forest" = 3, "Deciduous Broadleaf forest" = 4, "Shrub" = 5, "Grass" = 6, "Cereal crop" = 7, "Broadleaf crop" = 8, "Urban & built-up" = 9, "Snow & ice" = 10, "Barren/Sparsely vegetated" = 11, "Unclassified" = 254, "NoDataFill" = 255) ) NUM_METADATA_COLS <- 10 ########## if(!file.exists(Dir)) stop("Character string input for Dir argument does not resemble an existing file path.") file.list <- list.files(path = Dir, pattern = "MCD12Q1.*asc$") if(length(file.list) == 0) stop("Found no MODIS Land Cover ASCII files in Dir.") if(!any(GetBands("MCD12Q1") == Band)) stop("LandCover is for land cover data. Band specified is not for this product.") lc.type.set <- LC_CLASS[[which(names(LC_CLASS) == Band)]] NoDataFill <- unname(lc.type.set["NoDataFill"]) ValidRange <- unname(lc.type.set) lc.summary <- list(NA) for(i in 1:length(file.list)){ cat("Processing file ", i, " of ", length(file.list), "...\n", sep="") lc.subset <- read.csv(paste(Dir, "/", file.list[i], sep = ""), header = FALSE, as.is = TRUE) names(lc.subset) <- c("nrow", "ncol", "xll", "yll", "pixelsize", "row.id", "land.product.code", "MODIS.acq.date", "where", "MODIS.proc.date", 1:(ncol(lc.subset) - NUM_METADATA_COLS)) where.long <- regexpr("Lon", lc.subset$where[1]) where.samp <- regexpr("Samp", lc.subset$where[1]) where.land <- regexpr("Land", lc.subset$row.id) lat <- as.numeric(substr(lc.subset$where[1], 4, where.long - 1)) long <- as.numeric(substr(lc.subset$where[1], where.long + 3, where.samp - 1)) band.codes <- substr(lc.subset$row.id, where.land, nchar(lc.subset$row.id)) ifelse(any(grepl(Band, lc.subset$row.id)), which.are.band <- which(band.codes == Band), stop("Cannot find which rows in LoadDat are band data. Make sure the only ascii files in the directory are those downloaded from MODISSubsets.")) lc.tiles <- as.matrix(lc.subset[which.are.band,(NUM_METADATA_COLS+1):ncol(lc.subset)], nrow = length(which.are.band), ncol = length((NUM_METADATA_COLS+1):ncol(lc.subset))) if(!all(lc.tiles %in% ValidRange)) stop("Some values fall outside the valid range for the data band specified.") # Screen pixels in lc.tiles: pixels = NoDataFill, or whose corresponding pixel in qc.tiles < QualityThreshold. lc.tiles <- matrix(ifelse(lc.tiles != NoDataFill, lc.tiles, NA), nrow = length(which.are.band)) # Extract year and day from the metadata and make POSIXlt dates (YYYY-MM-DD), ready for time-series analysis. year <- as.numeric(substr(lc.subset$MODIS.acq.date, 2, 5)) day <- as.numeric(substr(lc.subset$MODIS.acq.date, 6, 8)) lc.subset$date <- strptime(paste(year, "-", day, sep = ""), "%Y-%j") # Initialise objects to store landscape summaries lc.mode.class <- rep(NA, nrow(lc.tiles)) lc.richness <- rep(NA, nrow(lc.tiles)) simp.even <- rep(NA, nrow(lc.tiles)) simp.d <- rep(NA, nrow(lc.tiles)) no.fill <- rep(NA, nrow(lc.tiles)) poor.quality <- rep(NA, nrow(lc.tiles)) for(x in 1:nrow(lc.tiles)){ # Calculate mode - most frequent lc class lc.freq <- table(lc.tiles[x, ]) lc.freq <- lc.freq / ncol(lc.tiles) lc.freq <- sum(lc.freq^2) # Calculate Simpson's D diversity index simp.d[x] <- 1 / lc.freq lc.mode <- which.max(table(lc.tiles[x, ])) lc.mode.class[x] <- names(which(lc.type.set == lc.mode)) # Calculate landscape richness lc.richness[x] <- length(table(lc.tiles[x, ])) # Calculate Simpson's measure of evenness simp.even[x] <- simp.d[x] / lc.richness[x] no.fill[x] <- paste(round((sum(lc.subset[x,(NUM_METADATA_COLS+1):ncol(lc.subset)] == NoDataFill) / length(lc.tiles[x, ])) * 100, 2), "% (", sum(lc.subset[x,(NUM_METADATA_COLS+1):ncol(lc.subset)] == NoDataFill), "/", length(lc.tiles[x, ]), ")", sep = "") } # End of loop that summaries tiles at each time-step, for the ith ASCII file. # Compile summaries into a table. lc.summary[[i]] <- data.frame(lat = lat, long = long, date = lc.subset$date[which(band.codes == Band)], modis.band = Band, most.common = lc.mode.class, richness = lc.richness, simpsons.d = simp.d, simpsons.evenness = simp.even, no.data.fill = no.fill) } # End of loop that reiterates for each ascii file. # Write output summary file by appending summary data from all files, producing one file of summary output. lc.summary <- do.call("rbind", lc.summary) write.table(lc.summary, file = paste(Dir, "/", "MODIS_Land_Cover_Summary ", Sys.Date(), ".csv", sep = ""), sep = ",", row.names = FALSE) cat("Done! Check the 'MODIS Land Cover Summary' output file.\n") }