scatterplot_biomes.R

library(rgdal)
library(raster)
library(ggplot2)

####directories#####
dir = "D:/Doutorado/Sanduiche/research/maxent-test/driving_and_obs_overlap/AllConFire_2010_2019/"
dir2 = "D:/Doutorado/Sanduiche/research/maxent-variables/2012-2021/"
dir3 = "D:/Doutorado/Malhas/"

dir_save = "D:/Doutorado/Sanduiche/research/scatterplots/"

##### reading files
biomes = readOGR(paste0(dir3, "biomas_wgs84.shp"))

raster_fire = brick(paste0(dir, "GFED4.1s_Burned_Fraction.nc"))
raster_ppt = brick(paste0(dir, "precip.nc"))
raster_soilM = brick(paste0(dir, "soilM.nc"))
raster_rhumid = brick(paste0(dir, "rhumid.nc"))

raster_tas = brick(paste0(dir, "tas.nc"))
raster_crop = brick(paste0(dir, "crop.nc"))
raster_vpd = brick(paste0(dir, "vpd.nc"))
raster_forest = brick(paste0(dir2, "Forest.nc"))
raster_savanna = brick(paste0(dir2, "Savanna.nc"))
raster_grass = brick(paste0(dir2, "Grassland.nc"))


####### subsetting for 2012
fire_2012 = raster_fire[[25:36]]
ppt_2012 =  raster_ppt[[25:36]]
soil_2012 = raster_soilM[[25:36]]
rhumid_2012 = raster_rhumid[[25:36]]

tas_2012 = raster_tas[[25:36]]
crop_2012 = raster_crop[[25:36]]
vpd_2012 = raster_vpd[[25:36]]

forest_2012 = raster_forest[[1]]
savanna_2012 = raster_savanna[[1]]
grassland_2012 = raster_grassland[[1]]




# Create an empty list to store the results
result_list_fire = vector("list", nlayers(fire_2012))

# Iterate over each layer of the raster stack
for (i in 1:nlayers(fire_2012)) {
        # Extract the current layer
        layer = fire_2012[[i]]
        
        # Perform the operation for the current layer
        pol = rasterToPolygons(layer)
        nearestbiome = over(pol, biomes, byid=TRUE)
        pol$Biome = nearestbiome$Bioma
        result_list_fire[[i]] = pol
        
}

# Create an empty list to store the results
result_list_ppt = vector("list", nlayers(ppt_2012))

# Iterate over each layer of the raster stack
for (i in 1:nlayers(ppt_2012)) {
        # Extract the current layer
        layer = ppt_2012[[i]]
        
        # Perform the operation for the current layer
        pol = rasterToPolygons(layer)
        nearestbiome = over(pol, biomes, byid=TRUE)
        pol$Biome = nearestbiome$Bioma
        result_list_ppt[[i]] = pol
        
}



#change name of shapefiles list - fire
new_column_names = c("Fire", "Biomes")  # Replace with your desired column names

# Iterate through the list of shapefiles
for (i in seq_along(result_list)) {
        colnames(result_list[[i]]@data) = new_column_names
}


##merge them using rbind - fire

rbind_fire = do.call(rbind, result_list)


#change name of shapefiles list - ppt
new_column_names_ppt = c("PPT", "Biomes")  # Replace with your desired column names

# Iterate through the list of shapefiles
for (i in seq_along(result_list_ppt)) {
        colnames(result_list_ppt[[i]]@data) = new_column_names_ppt
}


##merge them using rbind - ppt

rbind_ppt = do.call(rbind, result_list_ppt)


#plot 1 - coloring classes
ggplot(data = data.frame(x = rbind_ppt$PPT, y = rbind_fire$Fire, Class = rbind_fire$Biomes), 
       aes(x = x, y = y, color = Class)) +
        geom_point(size = 1) + xlab("Precipitation") +
        ylab("Burned area fraction") 

#plot 2 - separated plots


# Create a data frame with the combined data
combined_data = data.frame(PPT = rbind_ppt$PPT,
                            Fire = rbind_fire$Fire,
                            Class = rbind_fire$Biomes)

# Rename the class names
combined_data$Class = factor(combined_data$Class, levels = c("Amazônia", "Cerrado", "Caatinga","Mata Atlântica", "Pantanal", "Pampa" ), 
                             labels = c("Amazonia", "Cerrado", "Caatinga","Atlantic Forest", "Pantanal", "Pampa"))

# Create the scatter plot with facetting
scatter_plot <- ggplot(combined_data, aes(x = PPT, y = Fire, color = Class)) +
        geom_point(size = .7) +
        xlab("Precipitation") +
        ylab("Burned area fraction") +
        ylim(0, 0.10) +
        facet_wrap(~ Class, ncol = 3) +
        guides(color = FALSE)  # Remove the legend

# Display the scatter plot
print(scatter_plot)
ggsave(paste0(dir_save, "scatter_ppt.png"), plot = scatter_plot, width = 10, height = 6, dpi = 300)

#################
#################Forest
#################

# Create an empty list to store the results
result_list_forest = vector("list", nlayers(forest_2012))

# Iterate over each layer of the raster stack
for (i in 1:nlayers(forest_2012)) {
        # Extract the current layer
        layer = forest_2012[[i]]
        
        # Perform the operation for the current layer
        pol = rasterToPolygons(layer)
        nearestbiome = over(pol, biomes, byid=TRUE)
        pol$Biome = nearestbiome$Bioma
        result_list_forest[[i]] = pol
        
}

# Function to change column names
change_column_names <- function(result_list, new_column_names) {
        for (i in seq_along(result_list)) {
                colnames(result_list[[i]]@data) <- new_column_names
        }
        return(result_list)
}

# Change column names for fire and forest shapefiles
new_column_names_fire = c("Fire", "Biomes")
result_list_fire = change_column_names(result_list, new_column_names_fire)

new_column_names_forest = c("Forest", "Biomes")
result_list_forest = change_column_names(result_list_forest, new_column_names_forest)

# Merge shapefiles
rbind_fire <- do.call(rbind, result_list_fire)
rbind_forest <- do.call(rbind, result_list_forest)

#plot 2 - separated plots

# Create a data frame with the combined data
combined_data_for = data.frame(Forest = rbind_forest$Forest,
                           Fire = rbind_fire$Fire,
                           Class = rbind_fire$Biomes)

# Rename the class names
combined_data_for$Class = factor(combined_data$Class, levels = c("Amazônia", "Cerrado", "Caatinga","Mata Atlântica", "Pantanal", "Pampa" ), 
                             labels = c("Amazonia", "Cerrado", "Caatinga","Atlantic Forest", "Pantanal", "Pampa"))

# Create the scatter plot with facetting
scatter_plot_forest = ggplot(combined_data_for, aes(x = Forest, y = Fire, color = Class)) +
        geom_point(size = .7) +
        xlab("Forest cover (%)") +
        ylab("Burned area fraction") +
        ylim(0, 0.10) +
        facet_wrap(~ Class, ncol = 3) +
        guides(color = FALSE)  # Remove the legend

# Display the scatter plot
print(scatter_plot_forest)
ggsave(paste0(dir_save, "scatter_forest.png"), plot = scatter_plot_forest, width = 10, height = 6, dpi = 300)


#################Savanna

#################Grassland

#################
#################SoilM
#################

# Create an empty list to store the results
result_list_soil = vector("list", nlayers(soil_2012))

# Iterate over each layer of the raster stack
for (i in 1:nlayers(soil_2012)) {
        # Extract the current layer
        layer = soil_2012[[i]]
        
        # Perform the operation for the current layer
        pol = rasterToPolygons(layer)
        nearestbiome = over(pol, biomes, byid=TRUE)
        pol$Biome = nearestbiome$Bioma
        result_list_soil[[i]] = pol
        
}

# Function to change column names
change_column_names <- function(result_list, new_column_names) {
        for (i in seq_along(result_list)) {
                colnames(result_list[[i]]@data) <- new_column_names
        }
        return(result_list)
}

# Change column names for fire and forest shapefiles
new_column_names_fire = c("Fire", "Biomes")
result_list_fire = change_column_names(result_list_fire, new_column_names_fire)

new_column_names_soil = c("SoilM", "Biomes")
result_list_soil = change_column_names(result_list_soil, new_column_names_soil)

# Merge shapefiles
rbind_fire <- do.call(rbind, result_list_fire)
rbind_soil <- do.call(rbind, result_list_soil)

#plot 2 - separated plots

# Create a data frame with the combined data
combined_data_soil = data.frame(Soil = rbind_soil$SoilM,
                               Fire = rbind_fire$Fire,
                               Class = rbind_fire$Biomes)

# Rename the class names
combined_data_soil$Class = factor(combined_data_soil$Class, levels = c("Amazônia", "Cerrado", "Caatinga","Mata Atlântica", "Pantanal", "Pampa" ), 
                                 labels = c("Amazonia", "Cerrado", "Caatinga","Atlantic Forest", "Pantanal", "Pampa"))

# Create the scatter plot with facetting
scatter_plot_soil = ggplot(combined_data_soil, aes(x = Soil, y = Fire, color = Class)) +
        geom_point(size = .7) +
        xlab("Soil Moisture") +
        ylab("Burned area fraction") +
        ylim(0, 0.10) +
        facet_wrap(~ Class, ncol = 3) +
        guides(color = FALSE)  # Remove the legend

# Display the scatter plot
print(scatter_plot_soil)
ggsave(paste0(dir_save, "scatter_soilM.png"), plot = scatter_plot_soil, width = 10, height = 6, dpi = 300)


#################
#################Rhumid
#################

# Create an empty list to store the results
result_list_rhumid = vector("list", nlayers(rhumid_2012))

# Iterate over each layer of the raster stack
for (i in 1:nlayers(rhumid_2012)) {
        # Extract the current layer
        layer = rhumid_2012[[i]]
        
        # Perform the operation for the current layer
        pol = rasterToPolygons(layer)
        nearestbiome = over(pol, biomes, byid=TRUE)
        pol$Biome = nearestbiome$Bioma
        result_list_rhumid[[i]] = pol
        
}


# Change column names for fire and forest shapefiles using function

new_column_names_rhumid = c("rhumid", "Biomes")
result_list_rhumid = change_column_names(result_list_rhumid, new_column_names_rhumid)

# Merge shapefiles
rbind_rhumid <- do.call(rbind, result_list_rhumid)

#plot 2 - separated plots

# Create a data frame with the combined data
combined_data_rhumid = data.frame(rhumid = rbind_rhumid$rhumid,
                                Fire = rbind_fire$Fire,
                                Class = rbind_fire$Biomes)

# Rename the class names
combined_data_rhumid$Class = factor(combined_data_rhumid$Class, levels = c("Amazônia", "Cerrado", "Caatinga","Mata Atlântica", "Pantanal", "Pampa" ), 
                                  labels = c("Amazonia", "Cerrado", "Caatinga","Atlantic Forest", "Pantanal", "Pampa"))

# Create the scatter plot with facetting
scatter_plot_rhumid = ggplot(combined_data_rhumid, aes(x = rhumid, y = Fire, color = Class)) +
        geom_point(size = .7) +
        xlab("Relative Humidity") +
        ylab("Burned area fraction") +
        ylim(0, 0.10) +
        facet_wrap(~ Class, ncol = 3) +
        guides(color = FALSE)  # Remove the legend

# Display the scatter plot
print(scatter_plot_rhumid)
ggsave(paste0(dir_save, "scatter_rhumid.png"), plot = scatter_plot_rhumid, width = 10, height = 6, dpi = 300)


#################
#################tas
#################

# Create an empty list to store the results
result_list_tas = vector("list", nlayers(tas_2012))

# Iterate over each layer of the raster stack
for (i in 1:nlayers(tas_2012)) {
        # Extract the current layer
        layer = tas_2012[[i]]
        
        # Perform the operation for the current layer
        pol = rasterToPolygons(layer)
        nearestbiome = over(pol, biomes, byid=TRUE)
        pol$Biome = nearestbiome$Bioma
        result_list_tas[[i]] = pol
        
}


# Change column names for fire and forest shapefiles using function

new_column_names_tas = c("tas", "Biomes")
result_list_tas = change_column_names(result_list_tas, new_column_names_tas)

# Merge shapefiles
rbind_tas <- do.call(rbind, result_list_tas)

#plot 2 - separated plots

# Create a data frame with the combined data
combined_data_tas = data.frame(tas = rbind_tas$tas,
                                  Fire = rbind_fire$Fire,
                                  Class = rbind_fire$Biomes)

# Rename the class names
combined_data_tas$Class = factor(combined_data_tas$Class, levels = c("Amazônia", "Cerrado", "Caatinga","Mata Atlântica", "Pantanal", "Pampa" ), 
                                    labels = c("Amazonia", "Cerrado", "Caatinga","Atlantic Forest", "Pantanal", "Pampa"))

# Create the scatter plot with facetting
scatter_plot_tas = ggplot(combined_data_tas, aes(x = tas, y = Fire, color = Class)) +
        geom_point(size = .7) +
        xlab("Mean Temperature") +
        ylab("Burned area fraction") +
        ylim(0, 0.10) +
        facet_wrap(~ Class, ncol = 3) +
        guides(color = FALSE)  # Remove the legend

# Display the scatter plot
print(scatter_plot_tas)
ggsave(paste0(dir_save, "scatter_tas.png"), plot = scatter_plot_tas, width = 10, height = 6, dpi = 300)

#################
#################crop
#################

# Create an empty list to store the results
result_list_crop = vector("list", nlayers(crop_2012))

# Iterate over each layer of the raster stack
for (i in 1:nlayers(crop_2012)) {
        # Extract the current layer
        layer = crop_2012[[i]]
        
        # Perform the operation for the current layer
        pol = rasterToPolygons(layer)
        nearestbiome = over(pol, biomes, byid=TRUE)
        pol$Biome = nearestbiome$Bioma
        result_list_crop[[i]] = pol
        
}


# Change column names for fire and forest shapefiles using function

new_column_names_crop = c("crop", "Biomes")
result_list_crop = change_column_names(result_list_crop, new_column_names_crop)

# Merge shapefiles
rbind_crop = do.call(rbind, result_list_crop)

#plot 2 - separated plots

# Create a data frame with the combined data
combined_data_crop = data.frame(crop = rbind_crop$crop,
                               Fire = rbind_fire$Fire,
                               Class = rbind_fire$Biomes)

# Rename the class names
combined_data_crop$Class = factor(combined_data_crop$Class, levels = c("Amazônia", "Cerrado", "Caatinga","Mata Atlântica", "Pantanal", "Pampa" ), 
                                 labels = c("Amazonia", "Cerrado", "Caatinga","Atlantic Forest", "Pantanal", "Pampa"))

# Create the scatter plot with facetting
scatter_plot_crop = ggplot(combined_data_crop, aes(x = crop, y = Fire, color = Class)) +
        geom_point(size = .7) +
        xlab("Crop fraction") +
        ylab("Burned area fraction") +
        ylim(0, 0.10) +
        facet_wrap(~ Class, ncol = 3) +
        guides(color = FALSE)  # Remove the legend

# Display the scatter plot
print(scatter_plot_crop)
ggsave(paste0(dir_save, "scatter_crop.png"), plot = scatter_plot_crop, width = 10, height = 6, dpi = 300)


#################
#################vpd
#################

# Create an empty list to store the results
result_list_vpd = vector("list", nlayers(vpd_2012))

# Iterate over each layer of the raster stack
for (i in 1:nlayers(vpd_2012)) {
        # Extract the current layer
        layer = vpd_2012[[i]]
        
        # Perform the operation for the current layer
        pol = rasterToPolygons(layer)
        nearestbiome = over(pol, biomes, byid=TRUE)
        pol$Biome = nearestbiome$Bioma
        result_list_vpd[[i]] = pol
        
}



new_column_names_vpd = c("vpd", "Biomes")
result_list_vpd = change_column_names(result_list_vpd, new_column_names_vpd)

# Merge shapefiles
rbind_vpd <- do.call(rbind, result_list_vpd)

#plot 2 - separated plots

# Create a data frame with the combined data
combined_data_vpd = data.frame(vpd = rbind_vpd$vpd,
                                Fire = rbind_fire$Fire,
                                Class = rbind_fire$Biomes)

# Rename the class names
combined_data_vpd$Class = factor(combined_data_vpd$Class, levels = c("Amazônia", "Cerrado", "Caatinga","Mata Atlântica", "Pantanal", "Pampa" ), 
                                  labels = c("Amazonia", "Cerrado", "Caatinga","Atlantic Forest", "Pantanal", "Pampa"))

# Create the scatter plot with facetting
scatter_plot_vpd = ggplot(combined_data_vpd, aes(x = vpd, y = Fire, color = Class)) +
        geom_point(size = .7) +
        xlab("Vapour Pressure Deficit") +
        ylab("Burned area fraction") +
        ylim(0, 0.10) +
        facet_wrap(~ Class, ncol = 3) +
        guides(color = FALSE)  # Remove the legend

# Display the scatter plot
print(scatter_plot_vpd)
ggsave(paste0(dir_save, "scatter_vpd.png"), plot = scatter_plot_vpd, width = 10, height = 6, dpi = 300)