# This set of commands processes the Rosan et al. data in CDO/NCO then generates the plot in Fig. 3a of the manuscript in R.
# ROSAN et al. 2022 data file version to netCDF using CDO operators


#Pixel value 1 = decrease in BA & increase in fragmentation
#Pixel value 2 = increase in BA & increase in fragmentation
#Pixel value 3 = decrease in BA & decrease in fragmentation
#Pixel value 4 = increase in BA & decrease in fragmentation
#NA pixels didn’t present significant trend on one or both (BA, fragmentation), so were excluded from analysis. 


cdo setrtoc,0.5,1.5,-100 Cerrado_frag_fire_combination.nc temp1.nc #decrease in BA & increase in fragmentation
cdo setrtoc,2.5,3.5,-10000 temp1.nc temp2.nc #decrease in BA & decrease in fragmentation
cdo setrtoc,1.5,2.5,10000 temp2.nc temp3.nc #increase in BA & increase in fragmentation
cdo setrtoc,3.5,4.5,100 temp3.nc bla1.nc # increase in BA & decrease in fragmentation
rm -rf temp*
cdo setrtoc,0.5,1.5,-100 Amazon_frag_fire_combination.nc temp1.nc
cdo setrtoc,2.5,3.5,-10000 temp1.nc temp2.nc
cdo setrtoc,1.5,2.5,10000 temp2.nc temp3.nc
cdo setrtoc,3.5,4.5,100 temp3.nc bla2.nc
rm -rf temp*
cdo remapcon,Areas.nc bla1.nc temp1.nc
cdo remapcon,Areas.nc bla2.nc temp2.nc

cdo setrtoc,-10,10,0 temp1.nc tempx.nc
cdo setrtoc,10,1000,1 tempx.nc temp3.nc
cdo setrtoc,1000,20000,10 temp3.nc temp4.nc
cdo setrtoc,-1000,-10,-1 temp4.nc temp5.nc
cdo setrtoc,-20000,-1000,-10 temp5.nc tempa.nc
cdo setctomiss,0 tempa.nc bla3.nc

cdo setrtoc,-10,10,0 temp2.nc tempy.nc
cdo setrtoc,10,1000,1  tempy.nc temp6.nc
cdo setrtoc,1000,20000,10 temp6.nc temp7.nc
cdo setrtoc,-1000,-10,-1 temp7.nc temp8.nc
cdo setrtoc,-20000,-1000,-10 temp8.nc temp9.nc
cdo setctomiss,0 temp9.nc bla4.nc

rm -rf temp*
#ncap2 -s 'where (Areas>0) Areas=1;' Areas.nc ones.nc
#cdo griddes Areas.nc > grid
cdo setgrid,grid bla3.nc temp1.nc
cdo setgrid,grid bla4.nc temp2.nc
cdo setmissval,9.96921e+36 temp1.nc temp3.nc
cdo setmissval,9.96921e+36 temp2.nc temp4.nc

cdo mul ones.nc temp3.nc temp5.nc
cdo mul ones.nc temp4.nc temp6.nc
cdo chname,Areas,Cerrado temp5.nc temp7.nc
cdo chname,Areas,Amazon temp6.nc temp8.nc

cdo merge temp7.nc temp8.nc Frag_Brazil_remap.nc

cdo sellonlatbox,-90,-30,5,-25 Frag_Brazil_remap.nc Frag_Brazil_remap_box.nc
cdo merge AMAZON_dFRAC.nc Frag_Brazil_remap_box.nc Frag_Brazil_dFRAC_remap_box.nc





################## R Script
##################  Script to build Brazil fragmentation fire trend map figures
################## against simulation output


#Clean the dataset loaded
rm(list = ls());
#load the library needed
library(ncdf4);
library(maps)
library(ggplot2)
library(raster) # package for raster manipulation
# theme_set(theme_minimal() +
# theme(legend.position = "top"))


library(plyr)
library(dplyr)
library(eurostat)
library(ggmap)
library(tidyr)
library(geosphere)
library(ggplot2)
library(sp)
library(rgdal)
library(rgeos)
library(maps)
library(maptools)
library(countrycode)
library(rgdal) # package for geospatial analysis
library(lattice)
library(RColorBrewer)
library(maptools)
library(ggpubr)


#Read the file
filein1="Frag_Brazil_dFRAC_remap_box.nc"
nc <- nc_open(filein1)

print(paste("The file has",nc$nvars,"variables"))
print(nc)

lon <- ncvar_get(nc,"lon")
nlon <- dim(lon)
head(lon)
lat <- ncvar_get(nc,"lat")
nlat <- dim(lat)
head(lat)
print(c(nlon,nlat))

data(wrld_simpl)
# get time
time <- ncvar_get(nc,"time_counter")
time



# FIRST CONVERT DATA TO DATA FRAME
Cerrado <- ncvar_get(nc, "Cerrado")
#Cerrado[Cerrado == 0] <- NA
fillvalue <- ncatt_get(nc,"Cerrado", "_FillValue")
Cerrado[Cerrado==fillvalue$value] <- NA
library(maptools)
grid <- expand.grid(lon=lon, lat=lat) # CREATES DATA FRAME
cutpts <- c(-10,-5,0,5,10)
temp1<-levelplot(Cerrado ~ lon * lat, data=grid, at=cutpts, cuts=9, pretty=T, 
  col.regions=(rev(brewer.pal(10,"RdBu"))))
temp1
# create dataframe -- reshape data
# matrix (nlon*nlat rows by 2 cols) of lons and lats
lonlat <- as.matrix(expand.grid(lon,lat))
dim(lonlat)
# vector of `tmp` values
Cerrado_tmp <- as.vector(Cerrado)
length(Cerrado_tmp)
# create dataframe and add names
df_Cerrado <- data.frame(cbind(lonlat,Cerrado_tmp))
names(df_Cerrado) <- c("lon","lat","Cerrado")
head(na.omit(df_Cerrado), 10)
#

Amazon <- ncvar_get(nc, "Amazon")
#Amazon[Amazon == 0] <- NA
fillvalue <- ncatt_get(nc,"Amazon", "_FillValue")
Amazon[Amazon==fillvalue$value] <- NA
library(maptools)
grid <- expand.grid(lon=lon, lat=lat) # CREATES DATA FRAME
cutpts <- c(-10,-5,0,5,10)
temp2<-levelplot(Amazon ~ lon * lat, data=grid, at=cutpts, cuts=9, pretty=T, 
  col.regions=(rev(brewer.pal(10,"RdBu"))))
temp2
# create dataframe -- reshape data
# matrix (nlon*nlat rows by 2 cols) of lons and lats
lonlat <- as.matrix(expand.grid(lon,lat))
dim(lonlat)
# vector of `tmp` values
Amazon_tmp <- as.vector(Amazon)
length(Amazon_tmp)
# create dataframe and add names
df_Amazon <- data.frame(cbind(lonlat,Amazon_tmp))
names(df_Amazon) <- c("lon","lat","Amazon")
head(na.omit(df_Amazon), 10)

#



#
D_AREA_BURNT <- ncvar_get(nc, "D_AREA_BURNT")
D_AREA_BURNT[D_AREA_BURNT == 0] <- NA
# D_AREA_BURNT[D_AREA_BURNT >0] <- 1
# D_AREA_BURNT[D_AREA_BURNT <0] <- -1

fillvalue <- ncatt_get(nc,"D_AREA_BURNT", "_FillValue")
D_AREA_BURNT[D_AREA_BURNT==fillvalue$value] <- NA
library(maptools)
grid <- expand.grid(lon=lon, lat=lat) # CREATES DATA FRAME
cutpts <- c(-0.5,-0.4,-0.3,-0.2,-0.1,0,0.1,0.2,0.3,0.4,0.5)
temp5<-levelplot(D_AREA_BURNT ~ lon * lat, data=grid, at=cutpts, cuts=12, pretty=T,  
  col.regions=(rev(brewer.pal(10,"RdBu"))))
temp5
# create dataframe -- reshape data
# matrix (nlon*nlat rows by 2 cols) of lons and lats
lonlat <- as.matrix(expand.grid(lon,lat))
dim(lonlat)
# vector of `tmp` values
D_AREA_BURNT_tmp <- as.vector(D_AREA_BURNT)
length(D_AREA_BURNT_tmp)
# create dataframe and add names
df_D_AREA_BURNT <- data.frame(cbind(lonlat,D_AREA_BURNT_tmp))
names(df_D_AREA_BURNT) <- c("lon","lat","D_AREA_BURNT")
head(na.omit(df_D_AREA_BURNT), 10)
#


df_map <- data.frame(df_Cerrado, df_Amazon,df_D_AREA_BURNT)

library(plyr)
library(dplyr)
library(eurostat)
library(ggmap)
library(tidyr)
library(geosphere)
library(ggplot2)
library(sp)
library(rgdal)
library(rgeos)
library(maps)
library(maptools)
library(countrycode)
library(rgdal) # package for geospatial analysis
library(lattice)
library(RColorBrewer)
library(maptools)
library(ggpubr)



data("world.cities")
head(world.cities)
df <- world.cities
data(wrld_simpl)
plot(wrld_simpl)

data(wrld_simpl)
shape <- wrld_simpl
shape_df <- fortify(shape, region = 'ISO3')
#  test1 + scale_fill_distiller(palette = "Spectral", na.value="gray90", breaks = c(-1,0,1.0,1,1), limits = c(-2,2), oob = squish)+


#theme(panel.grid.major = element_line(colour = gray(0.5), linetype = "dashed", size = 0.5)+

# Northern S. America  
 test1<- geom_raster(data=df_D_AREA_BURNT, aes(x=lon, y=lat, fill = D_AREA_BURNT))
  ggplot() +
test1 + scale_fill_distiller(palette = "RdBu", na.value="gray98", breaks = c(-0.05,0.05), limits = c(-0.1,0.1), oob=scales::squish)+
theme(panel.grid.major = element_line(colour = gray(0.5), linetype = "dashed", size = 0.5))+
geom_point(data = df_map, aes(x = lon, y = lat,color=cut(Amazon,c(-12,-3,-0.5,0,0.5,3,12))),size=1,shape=16)+
scale_color_manual(name = "Amazon",values = c("(-12,-3]" = "#0000FF","(-3,-0.5]" = "#33FFFF","(0.5,3]" = "orange","(3,12]" = "#FF0000"),na.value = NA)+ 
# geom_point(data = df_map, aes(x = lon, y = lat,color=cut(Cerrado,c(-12,-3,-0.5,0,0.5,3,12))),size=1,shape=17)+
# scale_color_manual(name = "Cerrado",values = c("(-12,-3]" = "#0000FF","(-3,-0.5]" = "#33FFFF","(0.5,3]" = "orange","(3,12]" = "#FF3300"),na.value = NA)+ 
geom_polygon(data = shape , aes(x = long, y = lat, group = group), color = 'black', fill=NA) +
 theme(legend.position = "right", legend.key.width = unit(1.5, "cm"))+
  coord_equal(xlim = c(-82,-35), ylim = c(-23.5,3.5))

#####
# THIS VERSION WORKS
# Northern S. America  
 test1<- geom_raster(data=df_D_AREA_BURNT, aes(x=lon, y=lat, fill = D_AREA_BURNT))
  ggplot() +
test1 + scale_fill_distiller(palette = "RdBu", na.value="gray98", breaks = c(-0.05,0.05), limits = c(-0.1,0.1), oob=scales::squish)+
theme(panel.grid.major = element_line(colour = gray(0.5), linetype = "dashed", size = 0.5))+
geom_point(data = df_map, aes(x = lon, y = lat,color=cut(Amazon,c(-12,-3,-0.5,0,0.5,3,12))),size=1,shape=16)+
scale_color_manual(name = "Amazon",values = c("(-12,-3]" = "#0000FF","(-3,-0.5]" = "#33FFFF","(0.5,3]" = "orange","(3,12]" = "#FF0000"),na.value = NA)+ 
geom_point(data = df_map, aes(x = lon, y = lat,color=cut(Cerrado,c(-12,-3,-0.5,0,0.5,3,12))),size=1,shape=17)+
scale_color_manual(name = "Cerrado",values = c("(-12,-3]" = "#0000FF","(-3,-0.5]" = "#33FFFF","(0.5,3]" = "orange","(3,12]" = "#FF3300"),na.value = NA)+ 
geom_polygon(data = shape , aes(x = long, y = lat, group = group), color = 'black', fill=NA) +
 theme(legend.position = "bottom", legend.key.width = unit(1, "cm"))+
  coord_equal(xlim = c(-82,-35), ylim = c(-23.5,3.5))



  
  # Northern S. America  
 test1<- geom_raster(data=df_D_AREA_BURNT, aes(x=lon, y=lat, fill = D_AREA_BURNT))
  ggplot() +
  test1 + scale_fill_distiller(palette = "RdBu", na.value="gray98", breaks = c(-0.05,0,0.05), limits = c(-0.1,0.1), oob=scales::squish)+
theme(panel.grid.major = element_line(colour = gray(0.5), linetype = "dashed", size = 0.5))+
geom_point(data = df_map, aes(x = lon, y = lat,color=cut(Amazon,c(-12,-3,0,3,12))),size=1,shape=16)+
scale_color_manual(name = "Amazon",values = c("(-12,-3]" = "blue","(-3,0)]" = "white","(0,3]" = "orange","(3,12]" = "red"),na.value = NA)+ 
geom_point(data = df_map, aes(x = lon, y = lat,color=cut(Cerrado,c(-12,-3,0,3,12))),size=1,shape=17)+
scale_color_manual(name = "Cerrado",values = c("(-12,-3]" = "blue","(-3,0)]" = "white","(0,3]" = "orange","(3,12]" = "red"),na.value = NA)+ 
geom_polygon(data = shape , aes(x = long, y = lat, group = group), color = 'black', fill=NA) +
 theme(legend.position = "bottom", legend.key.width = unit(0.3, "cm"))+
  coord_equal(xlim = c(-82,-35), ylim = c(-23.5,3.5))

  
  
  
  
  
