Plot GAM values in a worlwide raster file

Question

I am trying to make a SDM (Species Distribution Models), to this I built a GAM model (let's say for species A:

  response<- SDM_data$presence
  predictors <- c("bioclim.bio2","bioclim.bio3", "bioclim.bio8","bioclim.bio9","bioclim.bio13","bioclim.bio14","bioclim.bio15","bioclim.bio18","bioclim.bio19", "Elevation", "Evapotranspiration", "Slope","Roads")
  formula <- as.formula(paste("presence ~", paste0("s(", predictors, ", k = 5)", collapse = "+")))
  
  GAM_model <- gam(formula, data = SDM_data, family = binomial, select = TRUE)

In that context, the "presence" refers to a binary factor structure where 1 represents presence and 0 represents absence (background points). This binary factor structure is combined with a set of numeric raster predictors, e.g., Elevation, Slope, and various bioclimatic variables. Note: The presence/absence data is generated globally. Furthermore, to avoid the overfitting of the model, I have restricted the number of knots to five.

  env_subset <- env[, predictors] #Select only the predictors used where env refer to a dataframe with the name of predictors
  
prediction <- predict(GAM_model, newdata = env_subset, type = "response", na.action = "na.pass")  # predict of my GAM model



#Empty raster file: 
    world_extent <- extent(-180, 179, -56, 84)
    prediction_raster <- raster(world_extent, resolution = 0.08333333)
        prediction_matrix <- matrix(prediction, nrow = nrow(prediction_raster), ncol = ncol(prediction_raster), byrow = TRUE)

 #Convert to matrix

#Note Warming message: In matrix(prediction, nrow = nrow(prediction_raster), ncol = ncol(prediction_raster), : data length [175565] is not a sub-multiple or multiple of the number of rows [1680]

prediction_raster <- raster(prediction_matrix, xmn = extent(prediction_raster)[1], xmx = extent(prediction_raster)[2],
                            ymn = extent(prediction_raster)[3], ymx = extent(prediction_raster)[4],
                            crs = projection(prediction_raster))


values(prediction_raster) <- prediction_matrix
plot(prediction_raster) #Plot

However, as you can see in the next images, it appears that the process of placing the prediction back onto the raster was not accurate...

This is how looks the str of my data:

tibble [11,234 x 30] (S3: tbl_df/tbl/data.frame)
 $ cell.id           : num [1:11234] 2250831 2655686 6345649 6098248 2246524 ...
 $ bioclim.bio1      : num [1:11234] 151 177 119 163 150 153 160 148 145 159 ...
 $ bioclim.bio2      : num [1:11234] 77 59 97 111 78 138 114 95 102 137 ...
 $ bioclim.bio3      : num [1:11234] 45 47 47 58 45 57 54 48 48 55 ...
 $ bioclim.bio4      : num [1:11234] 3291 2425 3660 2689 3326 ...
 $ bioclim.bio5      : num [1:11234] 235 245 228 254 233 273 272 258 263 263 ...
 $ bioclim.bio6      : num [1:11234] 66 120 23 63 63 32 64 61 53 18 ...
 $ bioclim.bio7      : num [1:11234] 169 125 205 191 170 241 208 197 210 245 ...
 $ bioclim.bio8      : num [1:11234] 114 161 72 170 112 160 125 123 139 197 ...
 $ bioclim.bio9      : num [1:11234] 189 201 150 135 188 186 202 193 194 104 ...
 $ bioclim.bio10     : num [1:11234] 191 210 165 197 190 197 203 193 194 200 ...
 $ bioclim.bio11     : num [1:11234] 106 149 72 129 104 106 117 102 96 104 ...
 $ bioclim.bio12     : num [1:11234] 981 612 1253 560 1007 ...
 $ bioclim.bio13     : num [1:11234] 137 101 124 55 140 53 139 76 81 139 ...
 $ bioclim.bio14     : num [1:11234] 11 3 83 36 11 28 22 44 49 7 ...
 $ bioclim.bio15     : num [1:11234] 53 70 12 11 53 17 60 17 16 73 ...
 $ bioclim.bio16     : num [1:11234] 396 293 351 152 406 133 403 213 227 364 ...
 $ bioclim.bio17     : num [1:11234] 64 19 265 123 66 93 77 136 154 22 ...
 $ bioclim.bio18     : num [1:11234] 74 30 278 145 77 103 82 136 154 354 ...
 $ bioclim.bio19     : num [1:11234] 388 245 351 125 397 120 395 188 195 22 ...
 $ Elevation         : num [1:11234] 9.91 194.43 412.11 321.96 23.72 ...
 $ Slope             : num [1:11234] 0.259 4.464 2.02 5.697 0.636 ...
 $ Evapotranspiration: num [1:11234] 10186 12878 0 10291 9359 ...
 $ Rugosity          : num [1:11234] 3.6 52.56 24.25 64.6 8.42 ...
 $ Roads             : num [1:11234] 0 8.25 3.67 11.08 1.33 ...
 $ Roads1            : num [1:11234] 1 8.25 3.67 11.08 1.33 ...
 $ species           : chr [1:11234] "Acacia mearnsii" "Acacia mearnsii" "Acacia mearnsii" "Acacia mearnsii" ...
 $ species_ID        : num [1:11234] 1 1 1 1 1 1 1 1 1 1 ...
 $ Order             : chr [1:11234] "Plants" "Plants" "Plants" "Plants" ...
 $ presence          : Factor w/ 2 levels "0","1": 2 2 2 2 2 2 2 2 2 2 ...

All comments and suggestions are more than welcome to understand how to place back the prediction of the GAM model in to a raster file

Plot GAM values in a worlwide raster file

Answers (1)

Related Questions