library(raster) library(sp) library(rgdal) library(st) library(sf) library(ncdf4) library(gtools) #### create summed up positive and negative anomalies to highlight regional variability #### read the anomaly rasters r = list.files(pattern='.tif') r r <- mixedsort(r) r <- stack(r) r # create negative and positive anomaly rasters and save them to files r_pos <- reclassify(r, cbind(-Inf, 0, 0), right=FALSE) plot(r_pos) r_pos nlayers(r_pos) ND <- unstack(r_pos) outputnames <- paste(seq_along(ND), "positive_anomalies.tif",sep="") for(i in seq_along(ND)){writeRaster(ND[[i]], file=outputnames[i])} r_neg <- reclassify(r, cbind(0, Inf, 0), right=FALSE) plot(r_neg) r_neg nlayers(r_neg) ND <- unstack(r_neg) outputnames <- paste(seq_along(ND), "negative_anomalies.tif",sep="") for(i in seq_along(ND)){writeRaster(ND[[i]], file=outputnames[i])} # load the negative and positive files in a list and create raster sums files_pos = list.files(pattern = "positive") # Get only the .tif files files_pos <- mixedsort(files_pos) files_pos POS <- stack(files_pos) POSsum <- calc(POS, sum) plot(POSsum) files_neg = list.files(pattern = "negative") # Get only the .tif files files_neg <- mixedsort(files_neg) files_neg NEG <- stack(files_neg) NEGsum <- calc(NEG, sum) plot(NEGsum) # write to rasters writeRaster(POSsum, filename="positive_anomaly.tiff", format="GTiff", overwrite=TRUE) writeRaster(NEGsum, filename="negative_anomaly.tiff", format="GTiff", overwrite=TRUE) # Loop through files and print mean for each raster for (i in files_neg) { j = cellStats(brick(i), stat='mean') print(j) } # Loop through files and print mean for each raster for (i in files_pos) { j = cellStats(brick(i), stat='mean') print(j) }