# Script to reproduce hypoxia models in the paper "Identifying areas prone to coastal hypoxia - the role of topography " # Virtanen EA et al., Biogeosciences. # Hypoxia_data does not include predictor variables, as they have been excluded due to national restrictions # Script here reproduces model for occasional hypoxia (<4.6 mg L -1), and can be modified accordingly for all 10 versions # and for other frequency-thresholded hypoxia models # Parameters shown are from the base model run of occasional hypoxia (<4.6 mg L -1) ######### # OH4_6 # ######### ###BASE### oxygen=read.dbf("data_source/Hypoxia_data.dbf") oxygen1$study_area=as.factor(oxygen1$study_area) str(oxygen1) oxygen=oxygen1 set.seed(4) inTrain <- createDataPartition(y = oxygen$OH4_6, p = .7, list = FALSE) training <- oxygen[ inTrain,] testing <- oxygen[-inTrain,] str(training) str(testing) mOH4_6<- gbm.step(data=training, gbm.x = 2:12, gbm.y = "OH4_6", family = "bernoulli", tree.complexity = 6, verbose=TRUE, learning.rate = 0.001, bag.fraction = 0.7) ############# # Thresholding (training) preds_mOH4_6<- predict.gbm(mOH4_6, training, n.trees=mOH4_6$gbm.call$best.trees, type="response") predResults1 <- data.frame(preds_mOH4_6) dat1=cbind(training, predResults1) dat1$plot_ID<- 1:nrow(dat1) data1=data.frame(dat1$plot_ID, dat1$OH4_6, dat1$preds_mOH4_6) optimal.thresholds(data1) # Thresholding (testing) preds_mOH4_6<- predict.gbm(mOH4_6, testing, n.trees=mOH4_6$gbm.call$best.trees, type="response") predResults1 <- data.frame(preds_mOH4_6) dat1=cbind(testing, predResults1) dat1$plot_ID<- 1:nrow(dat1) data1=data.frame(dat1$plot_ID, dat1$OH4_6, dat1$preds_mOH4_6) optimal.thresholds(data1) ############## # PREDICTORS # ############## ACR=raster("data_source/ACR.tif") SWMd20=raster("data_source/SWMd20.tif") depth20=raster("data_source/depth20.tif") TSI20=raster("data_source/TSI20.tif") BPI40_100=raster("data_source/BPI40_100.tif") VRM=raster("data_source/VRM.tif") BPI16_40=raster("data_source/BPI16_40.tif") BPI2_5=raster("data_source/BPI2_5.tif") BPI6_15=raster("data_source/BPI6_15.tif") BPI10_25=raster("data_source/BPI10_25.tif") study_area=raster("data_source/study_area.tif") # stack rasters envars=stack(ACR,SWMd20,depth20,TSI20,BPI40_100,VRM, BPI16_40, BPI2_5,BPI6_15, BPI10_25,study_area) # Predict (to a raster stack) pOH4_6=predict(envars, mOH4_6, n.trees=mOH4_6$gbm.call$best.trees, type="response", progress='text') # write rasters writeRaster(pOH4_6, file="data_source/pOH4_6.tif",format="GTiff", progress="text") # DONE!!