It would be more integrated with the previous functions if it incorporates the propagation of the parameter uncertainty.
For that purpose, I suggest repeating the operation for each row of the thermal_suitability_bounds data frame, (one row per bootstrapped simulations).
The output of the function would be a raster with the total months per year with optimal temperatures for the species by averaging the value of all the obtained rasters. It may also provide an option for an uncertainty raster whose values represent the standard deviation based on the boostrapped estimates of thermal bounds, rather than the number of months.
It would be more integrated with the previous functions if it incorporates the propagation of the parameter uncertainty.
For that purpose, I suggest repeating the operation for each row of the thermal_suitability_bounds data frame, (one row per bootstrapped simulations).
The output of the function would be a raster with the total months per year with optimal temperatures for the species by averaging the value of all the obtained rasters. It may also provide an option for an uncertainty raster whose values represent the standard deviation based on the boostrapped estimates of thermal bounds, rather than the number of months.