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")
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.