-
David Minor authoreda709bcb3
biomass_support.R 11.29 KiB
mosaic_list <- function(x, fun, datatype, format, options, overwrite,
tolerance=0.05, filename="") {
mosaic_args <- x
if (!missing(fun)) mosaic_args$fun <- fun
if (!missing(tolerance)) mosaic_args$tolerance <- tolerance
if (!missing(datatype)) mosaic_args$datatype <- datatype
if (!missing(format)) mosaic_args$format <- format
if (!missing(options)) mosaic_args$options <- options
if (!missing(overwrite)) mosaic_args$overwrite <- overwrite
mosaic_args$filename <- filename
do.call(mosaic, mosaic_args)
}
# Read in CCI biomass within a given area of interest. Multiple years can be provided as a vector
# Combines AGBD and SD into as layers in one raster mosaic
readCCI <- function(data_folder, year="2020", version=3, aoi, verbose=FALSE){
# 2010, 2017, 2018 available for CCI Biomass v2.0
# 2020 available for CCI Biomass v3.0
year <- as.character(year)
if(version ==2){
#v2.0 file format
file_root <- paste0('_ESACCI-BIOMASS-L4-AGB-MERGED-100m-', year, '-fv2.0.tif')
file_root_sd <- paste0('_ESACCI-BIOMASS-L4-AGB_SD-MERGED-100m-', year, '-fv2.0.tif')
}
if(version==3){
#v3.0 file format
file_root <- paste0('_ESACCI-BIOMASS-L4-AGB-MERGED-100m-', year, '-fv3.0.tif')
file_root_sd <- paste0('_ESACCI-BIOMASS-L4-AGB_SD-MERGED-100m-', year, '-fv3.0.tif')
}
filename=""
tile_stacks <- c()
for (n in 1:length(tiles)) {
#for (n in c(4,5,7,8)) { #Temporary restriction to tiles while bugs are worked out with rasters
#for (n in 2) { #Temporary restriction to tiles while bugs are worked out with rasters
tile <- tiles[n]
min_x <- bbox(tile)[1, 1]
max_y <- bbox(tile)[2, 2]
if (min_x < 0) {
min_x <- paste0('W', sprintf('%03i', abs(min_x)))
} else {
min_x <- paste0('E', sprintf('%03i', min_x))
}
if (max_y < 0) {
max_y <- paste0('S', sprintf('%02i', abs(max_y)))
} else {
max_y <- paste0('N', sprintf('%02i', max_y))
}
file_prefix <- paste0(max_y, min_x)
band_names2 <- c(paste0('agbd_', year), paste0('agbd_sd_', year))
filenames2 <- c(file.path(data_folder, year, paste0(file_prefix, file_root)), file.path(data_folder, year, paste0(file_prefix, file_root_sd)))
if(verbose){
print(filenames2)
}
if (!file.exists(filenames2[1])){
print(paste0(filenames2, " does not exist"))
next()
}
tile_stack <- crop(stack(filenames2), aoi, datatype="INT1U", format='GTiff', options="COMPRESS=LZW")
names(tile_stack) <- band_names2
tile_stacks <- c(tile_stacks, list(tile_stack))
}
#If there are multiple tiles, create mosaic. Otherwise mosaic only includes 1 tile
if (length(tile_stacks) > 1) {
tile_mosaic <- mosaic_list(tile_stacks, fun='mean', filename=filename,
datatype='INT1U', format='GTiff',
options='COMPRESS=LZW')
} else {
tile_mosaic <- tile_stacks[[1]]
if (filename != '') {
tile_mosaic <- writeRaster(tile_mosaic, filename=filename,
datatype="INT1U", format="GTiff",
options="COMPRESS=LZW")
}
}
names(tile_mosaic) <- band_names2
NAvalue(tile_mosaic) <- -1
cci <- tile_mosaic
return(cci)
}
readGFC <- function(data_folder, image=c('treecover2000', 'treecover2020', 'lossyear', 'gain'), aoi, verbose=FALSE){
dataset <- 'GFC-2020-v1.8'
file_root <- paste0('Hansen_', dataset, '_')
filename=""
tile_stacks <- c()
for (n in 1:length(tiles)) {
#for (n in 2) {
tile <- tiles[n]
min_x <- bbox(tile)[1, 1]
max_y <- bbox(tile)[2, 2]
if (min_x < 0) {
min_x <- paste0(sprintf('%03i', abs(min_x)), 'W')
} else {
min_x <- paste0(sprintf('%03i', min_x), 'E')
}
if (max_y < 0) {
max_y <- paste0(sprintf('%02i', abs(max_y)), 'S')
} else {
max_y <- paste0(sprintf('%02i', max_y), 'N')
}
file_suffix <- paste0('_', max_y, '_', min_x, '.tif')
image_names2 <- image
band_names2 <- image_names2
filenames2 <- file.path(data_folder, paste0(file_root, image_names2, file_suffix))
if(verbose){
print(filenames2)
}
if (!file.exists(filenames2[1])){
print(paste0(filenames2, " does not exist"))
next()
}
#tile_stack <- try(crop(stack(filenames2), aoi, datatype='INT1U', format='GTiff', options="COMPRESS=LZW"))
tile_stack <- crop(stack(filenames2), aoi, datatype='INT1U', format='GTiff', options="COMPRESS=LZW")
names(tile_stack) <- band_names2 tile_stacks <- c(tile_stacks, list(tile_stack))
}
if (length(tile_stacks) > 1) {
tile_mosaic <- mosaic_list(tile_stacks, fun='mean', filename=filename,
datatype='INT1U', format='GTiff',
options='COMPRESS=LZW')
} else {
tile_mosaic <- tile_stacks[[1]]
if (filename != '') {
tile_mosaic <- writeRaster(tile_mosaic, filename=filename,
datatype="INT1U", format="GTiff",
options="COMPRESS=LZW")
}
}
names(tile_mosaic) <- band_names2
NAvalue(tile_mosaic) <- -1
gfc <- tile_mosaic
return(gfc)
}
recode_gfc <- function(treecover2000, treecover2020, lossyear, gain) {
forest_threshold <- 30
forest2000 <- treecover2000 > forest_threshold
forest2020 <- treecover2020 > forest_threshold
# Note forest2000 is a binary variable
lossyear_recode <- lossyear * forest2000
gain_recode <- gain & (!forest2000)
lossgain <- gain & (lossyear != 0)
thresholded <- matrix(c(forest2000, forest2020, lossyear_recode, gain_recode, lossgain),
nrow=length(forest2000),
ncol=5)
return(thresholded)
}
recode_gfc2020 <- function(treecover2020) {
forest_threshold <- 30
forest2020 <- treecover2020 > forest_threshold
thresholded <- matrix(c(forest2020),
nrow=length(forest2020),
ncol=1)
return(thresholded)
}
agbd_table <- function(year=2020, agb_thresholded, agb_masked){
scale_factor=0.0001
agb_boundpoly <- as(extent(agb_thresholded), 'SpatialPolygons')
proj4string(agb_boundpoly) <- proj4string(agb_thresholded)
agb_boundpoly_wgs84 <- spTransform(agb_boundpoly, CRS('+init=epsg:4326'))
aoi_wgs84 <- spTransform(aoi, CRS('+init=epsg:4326'))
if (!gIntersects(agb_boundpoly_wgs84, aoi_wgs84)) {
stop('aoi does not intersect supplied Biomass extract')
}
if ((((xmin(agb_thresholded) >=-180) & (xmax(agb_thresholded) <=180)) || ((xmin(agb_thresholded) >=0) & (xmax(agb_thresholded) <=360))) &&
(ymin(agb_thresholded) >=-90) & (ymax(agb_thresholded) <= 90)) {
# Use the included calc_pixel_area function to calculate the area of
# one cell in each line of the raster, allowing for accurate areal
# estimates in square meters even when imagery is in
# WGS84.
message('Data appears to be in latitude/longitude. Calculating cell areas on a sphere.')
spherical_areas <- TRUE # Calculate the area of a single pixel in each line of the image (to
# avoid repeating this calculation later on)
x <- cci
xleft <- xmin(x)
xright <- xmin(x) + xres(x)
ylower <- seq(from=(ymax(x) - yres(x)), by=(-yres(x)), length.out=nrow(x))
yupper <- seq(from=ymax(x), by=(-yres(x)), length.out=nrow(x))
poly_areas <- function(xl, xr, yl, yu) {
areaPolygon(matrix(c(xl, yl,
xr, yl,
xr, yu,
xl, yu), ncol=2, byrow=TRUE))
}
pixel_areas <- mapply(poly_areas, xleft, xright, ylower, yupper)
} else {
spherical_areas <- FALSE
pixel_areas <- xres(agb_thresholded) * yres(agb_thresholded)
}
aoi <- spTransform(aoi, CRS(proj4string(agb_thresholded)))
if (!('label' %in% names(aoi))) {
aoi$label <- paste('AOI', seq(1:nrow(aoi@data)))
}
uniq_aoi_labels <- unique(aoi$label)
years <- as.numeric(year)
agb_table <- data.frame(year=rep(years, length(uniq_aoi_labels)),
aoi=rep(uniq_aoi_labels, each=length(years)))
agb_table$agb <- 0
agb_table$agbd <- 0
agb_table$area <- 0
agb_table$n_pixel <- 0
#agb_table$agbd_sd <- 0
n_years <- length(years)
for (n in 1:nrow(aoi)) {
#cci_masked <- mask(cci_thresholded, aoi[n, ], datatype='INT1U', format="GTiff",
# options="COMPRESS=LZW")
# Find the first row in the loss table that applies to this aoi
agb_table_st_row <- match(aoi[n, ]$label, agb_table$aoi)
# No loss in first year
#agb_table$agb[loss_table_st_row] <- NA
# Process by blocks to avoid unnecessary reads from disk
bs <- blockSize(agb_masked)
for (block_num in 1:bs$n) {
bl_st_row <- bs$row[block_num]
bl_nrows <- bs$nrows[block_num]
agb_bl <- getValuesBlock(agb_masked, bl_st_row, bl_nrows)
if (spherical_areas) {
bl_pixel_areas <- rep(pixel_areas[bl_st_row:(bl_st_row + bl_nrows - 1)], each=ncol(agb_thresholded))
} else {
bl_pixel_areas <- pixel_areas
}
# Calculate initial aboveground biomass, note that areas are converted to square
# meters using pixel size, then converted using scale_factor
# (default output units are Mg)
# agbinit_col <- grep(paste0('agb_',year[1]), names(cci))
# agb_table$agb[agb_table_st_row] <- agb_table$agb[agb_table_st_row] + # sum(cci_bl[, agbinit_col] * bl_pixel_areas * scale_factor, na.rm=TRUE)
# agb_table$agbd[agb_table_st_row + i] <- agb_table$agbd[agb_table_st_row + i] +
# mean((cci_bl[, agb_col] == i), na.rm=TRUE)
for (i in 0:(n_years-1)) {
#agbd_col <- grep(paste0('agbd_',year[i+1]), names(cci_thresholded))
#sd_col <- grep(paste0('agbd_sd_',year[i+1]), names(cci_thresholded))
# n + 1 because first row is year 2010
agb_table$agb[agb_table_st_row + i] <- agb_table$agb[agb_table_st_row + i] +
sum(agb_bl * bl_pixel_areas * scale_factor, na.rm=TRUE)
agb_table$area[agb_table_st_row + i] <- agb_table$area[agb_table_st_row + i] +
sum((!is.na(agb_bl))*bl_pixel_areas * scale_factor, na.rm=TRUE)
agb_table$n_pixel[agb_table_st_row + i] <- agb_table$n_pixel[agb_table_st_row + i] +
sum(!is.na(agb_bl))
if(bl_st_row==1){
agb_table$agbd[agb_table_st_row + i] <- mean(agb_bl, na.rm=TRUE)
# agb_table$agbd_sd[agb_table_st_row + i] <- mean(cci_bl[, sd_col], na.rm=TRUE)
} else{
agb_table$agbd[agb_table_st_row + i] <- weighted.mean(c(agb_table$agbd[agb_table_st_row + i],
mean(agb_bl, na.rm=TRUE)), w=c((bl_st_row-1), bl_nrows), na.rm=TRUE)
# agb_table$agbd_sd[agb_table_st_row + i] <- weighted.mean(c(agb_table$agbd_sd[agb_table_st_row + i],
# mean(cci_bl[, sd_col], na.rm=TRUE)), w=c((bl_st_row-1), bl_nrows), na.rm=TRUE)
}
}
}
}
return(agb_table)
}