-
R语言：利用biomod2进行生态位建模
在这里主要是分享一个不错的代码，喜欢的可以慢慢研究。我看了一遍，觉得里面有很多有意思的东西，供大家学习和参考。
利用PCA轴总结的70个环境变量，利用biomod2进行生态位建模：
#----------------------------------------------------------# # NICHE MODELLING WITH BIOMOD2 USING ####### # 70 ENVIRONMENTAL VARIABLES (10-km RESOLUTION) ####### # SUMMARIZED IN PCA AXES ####### #-------------------------------------------------------# # Contact: Pedro V. Eisenlohr (pedro.eisenlohr@unemat.br) #------------------------------------------------- Acknowledgments ------------------------------------------------------------#### ### Dr. Guarino Colli's team of Universidade de BrasÃlia. ######################################################################### ### Dr. Diogo Souza Bezerra Rocha (Instituto de Pesquisas Jardim Botânico/RJ). #################################################### ### Drª Marinez Ferreira de Siqueira (Instituto de Pesquisas Jardim Botânico/RJ). ################################################# ### My students of Ecology Lab, mainly J.C. Pires-de-Oliveira. #################################################################### #----------------------------------------- ---------------------------------------------------------------------------------------# #------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------# ### Environmental data source (70 variables): ### Temperature and precipitation (19 variables): CHELSA (http://chelsa-climate.org/). #Bio1 = Annual Mean Temperature #Bio2 = Mean Diurnal Range #Bio3 = Isothermality #Bio4 = Temperature Seasonality #Bio5 = Max Temperature of Warmest Month #Bio6 = Min Temperature of Coldest Month #Bio7 = Temperature Annual Range #Bio8 = Mean Temperature of Wettest Quarter #Bio9 = Mean Temperature of Driest Quarter #Bio10 = Mean Temperature of Warmest Quarter #Bio11 = Mean Temperature of Coldest Quarter #Bio12 = Annual Precipitation #Bio13 = Precipitation of Wettest Month #Bio14 = Precipitation of Driest Month #Bio15 = Precipitation Seasonality #Bio16 = Precipitation of Wettest Quarter #Bio17 = Precipitation of Driest Quarter #Bio18 = Precipitation of Warmest Quarter #Bio19 = Precipitation of Coldest Quarter ### Solar radiation (3 variables), water vapor pressure (3 variables) and wind speed (3 variables): WorldClim 2.0 (http://worldclim.org/version2). #Solar Radiation: Maximum, Minimum and Mean #Water Vapor Pressure: Maximum, Minimum and Mean #Wind Speed: Maximum, Minimum and Mean ### Cloud Cover (3 variables): CRU-TS v3.10.01 Historic Climate Database for GIS (http://www.cgiar-csi.org/data/uea-cru-ts-v3-10-01-historic-climate-database). #Cloud Cover: Maximum, Minimum and Mean ### Enhanced Vegetation Index (3 variables): http://www.earthenv.org/. #Coefficient of variation of EVI = Normalized dispersion of EVI #Range of EVI #Standard deviation of EVI ### Forest Coverage (1 variable): http://www.fao.org/soils-portal/soil-survey/soil-maps-and-databases/harmonized-world-soil-database-v12/en/. #Forest land, calibrated to FRA2000 land statistics ### Grassland/Scrub/Woodland Coverage (1 variable): http://www.fao.org/soils-portal/soil-survey/soil-maps-and-databases/harmonized-world-soil-database-v12/en/. ### Water Bodies Coverage (1 variable): http://www.fao.org/soils-portal/soil-survey/soil-maps-and-databases/harmonized-world-soil-database-v12/en/. ### Elevation (1 variable): CGIAR-CSI (2006): NASA Shuttle Radar Topographic Mission (SRTM) (http://srtm.csi.cgiar.org/). ### Slope (1 variable) and Aspect (1 variable): obtained from Elevation. #Topographic variables obtained by applying 'terrain' function of 'raster' package. ### Topographic Wetness Index (1 variable): ENVIREM - ENVIronmental Rasters for Ecological Modeling (http://envirem.github.io/#varTable). ### Global Relief Model (1 variable): UNEP - http://geodata.grid.unep.ch/results.php #Global Relief Model of Earth's surface that integrates land topography and ocean bathymetry. ### Terrain Roughness Index (1 variable): ENVIREM - ENVIronmental Rasters for Ecological Modeling (http://envirem.github.io/#varTable). ### Potential Evapotranspiration - PET (6 variables) and Aridity Index (1 variable): Global Aridity and PET Database (http://www.cgiar-csi.org/data/global-aridity-and-pet-database) # and ENVIREM - ENVIronmental Rasters for Ecological Modeling (http://envirem.github.io/#varTable). #Annual Potential Evapotranspiration. #Mean Monthly PET of Coldest Quarter. #Mean Monthly PET of Driest Quarter. #PET Seasonality: monthly variability in potential evapotranspiration. #Mean Monthly PET of Warmest Quarter. #Mean Monthly PET of Wettest quarter #Global Annual Aridity Index. ### AET (1 variable) and Soil Water Stress (3 variables): Global High-Resolution Soil-Water Balance (http://www.cgiar-csi.org/data/global-high-resolution-soil-water-balance#download). #Mean Annual Actual Evapotranspiration. #Soil Water Stress: Maximum, Minimum and Mean. ### Relative Humidity (6 variables): Climond (https://www.climond.org/RawClimateData.aspx). #Relative Humidity at 9 am: Maximum, Minimum and Mean. #Relative Humidity at 3 pm: Maximum, Minimum and Mean. ### Soil Variables (10 variables): Soil grids (https://soilgrids.org) #BulkDensity = Bulk density (fine earth) in kg/cubic–meter #Clay = Clay content (0–2 micro meter) mass fraction in % #Coarse = Coarse fragments volumetric in % #Sand = Sand content (50–2000 micro meter) mass fraction in % #Silt = Silt content (2–50 micro meter) mass fraction in % #BDRLOG = Predicted probability of occurrence (0–100%) of R horizon #BDRICM = Depth to bedrock (R horizon) up to 200 cm #CARBON = Soil organic carbon content (fine earth fraction) in g per kg #pH_H20 = Soil pH x 10 in H2O #CEC = Cation exchange capacity of soil in cmolc/kg #------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------# #----------------------------# ## SET WORKING DIRECTORY #### #--------------------------# # Each user should adjust this! setwd(choose.dir()) getwd() list.files() # Among the listed files, there must be one called # "Environmental layers" and another called "Shapefiles". #---------------------------------------------# ## INSTALL AND LOAD THE REQUIRED PACKAGES #### #-------------------------------------------# #install.packages("biomod2", dep=T) #install.packages("colorRamps", dep=T) #install.packages("dismo", dep=T) #install.packages("dplyr", dep=T) #install.packages("maps", dep=T) #install.packages("maptools", dep=T) #install.packages("plotKML", dep=T) #install.packages("raster", dep=T) #install.packages("rgdal", dep=T) #install.packages("RStoolbox", dep=T) #install.packages("foreach", dep=T) #install.packages("doParallel", dep=T) library(biomod2) library(colorRamps) library(dismo) library(dplyr) library(maps) library(maptools) library(plotKML) library(raster) library(rgdal) library(RStoolbox) library(foreach) library(doParallel) library(virtualspecies) library(filesstrings) # Creating output folder # if (dir.exists("outputs") == F) { dir.create("outputs") } # Parallel processing # # cores <- detectCores()/2 # Assigning 50% of the cores for modeling #getDoParWorkers() #cl <- parallel::makeCluster(cores, outfile =paste0("./outputs/", "Log.log")) #registerDoParallel(cl) #getDoParWorkers() #--------------------------------------------------------------------------------------------# ### IF YOU HAVE ALREADY DOWNLOADED AND TREATED ALL LAYERS, YOU SHOULD SKIP THE STEPS BELOW #### #------------------------------------------------------------------------------------------# #---------------------------------------------------------------------# # Loading CHELSA layers (Temperature and Precipitation - 1979-2013) #### #---------------------------------------------------------------------# # First, load a 10-km resolution mask to resample: #bio.wc <- list.files("./Environmental layers/WorldClim 2.0", full.names=TRUE) #bio.wc <- stack(bio.wc) #bio.wc #res(bio.wc) # Crop mask layers #neotrop <- readOGR("./Shapefiles/ShapeNeo/neotropic.shp") #bio.wc <- mask(crop(bio.wc,neotrop),neotrop) #bio.wc #res(bio.wc) # Resampling CHELSA layers #bioclim <- list.files("./Environmental layers/CHELSA", full.names=TRUE, pattern=".grd") #bioclim <- stack(bioclim) #bioclim <- mask(crop(bioclim,neotrop),neotrop) #names(bioclim) #res(bioclim) #bioclim <-resample(bioclim, bio.wc) #res(bioclim) #plot(bioclim[[1]]) #names(bioclim) #bio1<-(bioclim[[1]]) #writeRaster(bio1, "bio01") #bio10<-(bioclim[[2]]) #writeRaster(bio10,"bio10") #bio11<-(bioclim[[3]]) #writeRaster(bio11,"bio11") #bio12<-(bioclim[[4]]) #writeRaster(bio12,"bio12") #bio13<-(bioclim[[5]]) #writeRaster(bio13,"bio13") #bio14<-(bioclim[[6]]) #writeRaster(bio14,"bio14") #bio15<-(bioclim[[7]]) #writeRaster(bio15,"bio15") #bio16<-(bioclim[[8]]) #writeRaster(bio16,"bio16") #bio17<-(bioclim[[9]]) #writeRaster(bio17,"bio17") #bio18<-(bioclim[[10]]) #writeRaster(bio18,"bio18") #bio19<-(bioclim[[11]]) #writeRaster(bio19,"bio19") #bio2<-(bioclim[[12]]) #writeRaster(bio2,"bio2") #bio3<-(bioclim[[13]]) #writeRaster(bio3,"bio3") #bio4<-(bioclim[[14]]) #writeRaster(bio4,"bio4") #bio5<-(bioclim[[15]]) #writeRaster(bio5,"bio5") #bio6<-(bioclim[[16]]) #writeRaster(bio6,"bio6") #bio7<-(bioclim[[17]]) #writeRaster(bio7,"bio7") #bio8<-(bioclim[[18]]) #writeRaster(bio8,"bio8") #bio9<-(bioclim[[19]]) #writeRaster(bio9,"bio9") #----------------------------------------# #----------------------------------------# ### Compiling other rasters to stack #### #--------------------------------------# #Solar Radiation: #solar.radiation <- list.files("./Environmental layers/Solar Radiation", pattern=".tif", full.names=TRUE) #solar.radiation <- stack(solar.radiation) #solar.radiation.mean <- mean(solar.radiation) #solar.radiation.max <- max(solar.radiation) #solar.radiation.min <- min(solar.radiation) #solar.radiation.mean <- mask(crop(solar.radiation.mean, neotrop),neotrop) #writeRaster(solar.radiation.mean,"SolarRadiationMean") #solar.radiation.max <- mask(crop(solar.radiation.max, neotrop),neotrop) #writeRaster(solar.radiation.max,"SolarRadiationMax") #solar.radiation.min <- mask(crop(solar.radiation.min, neotrop),neotrop) #writeRaster(solar.radiation.min,"SolarRadiationMin") #res(solar.radiation.mean) #plot(solar.radiation.mean) #res(solar.radiation.max) #plot(solar.radiation.max) #res(solar.radiation.min) #plot(solar.radiation.min) #Water Vapor Pressure: #water.vapor.pressure <- list.files("./Environmental layers/Water Vapor Pressure", pattern=".tif", full.names=TRUE) #water.vapor.pressure <-stack(water.vapor.pressure) #water.vapor.pressure.mean <-mean(water.vapor.pressure) #water.vapor.pressure.max <-max(water.vapor.pressure) #water.vapor.pressure.min <-min(water.vapor.pressure) #water.vapor.pressure.mean <- mask(crop(water.vapor.pressure.mean, neotrop),neotrop) #writeRaster(water.vapor.pressure.mean,"WaterVaporPressureMean") #water.vapor.pressure.max <- mask(crop(water.vapor.pressure.max, neotrop),neotrop) #writeRaster(water.vapor.pressure.max,"WaterVaporPressureMax") #water.vapor.pressure.min <- mask(crop(water.vapor.pressure.min, neotrop),neotrop) #writeRaster(water.vapor.pressure.min,"WaterVaporPressureMin") #res(water.vapor.pressure.mean) #plot(water.vapor.pressure.mean) #res(water.vapor.pressure.max) #plot(water.vapor.pressure.max) #res(water.vapor.pressure.min) #plot(water.vapor.pressure.min) #Wind Speed: #wind.speed <- list.files("./Environmental layers/Wind Speed", pattern=".tif", full.names=TRUE) #wind.speed <- stack(wind.speed) #wind.speed.mean <-mean(wind.speed) #wind.speed.max <-max(wind.speed) #wind.speed.min <-min(wind.speed) #wind.speed.mean <-mask(crop(wind.speed.mean, neotrop),neotrop) #writeRaster(wind.speed.mean, "WindSpeedMean") #wind.speed.max <-mask(crop(wind.speed.max, neotrop),neotrop) #writeRaster(wind.speed.max, "WindSpeedMax") #wind.speed.min <-mask(crop(wind.speed.min, neotrop),neotrop) #writeRaster(wind.speed.min, "WindSpeedMin") #res(wind.speed.mean) #plot(wind.speed.mean) #res(wind.speed.max) #plot(wind.speed.max) #res(wind.speed.min) #plot(wind.speed.min) #Cloud Cover: #cloud.cover<-list.files("./Environmental layers/Cloud Cover",pattern=".asc", full.names=TRUE) #cloud.cover<-stack(cloud.cover) #cloud.cover.mean<-mean(cloud.cover) #cloud.cover.max<-max(cloud.cover) #cloud.cover.min<-min(cloud.cover) #cloud.cover.mean<-mask(crop(cloud.cover.mean, neotrop),neotrop) #cloud.cover.mean<-resample(cloud.cover.mean,bioclim) #writeRaster(cloud.cover.mean,"CloudCoverMean") #cloud.cover.max<-mask(crop(cloud.cover.max, neotrop),neotrop) #cloud.cover.max<-resample(cloud.cover.max,bioclim) #writeRaster(cloud.cover.max,"CloudCoverMax") #cloud.cover.min<-mask(crop(cloud.cover.min, neotrop),neotrop) #cloud.cover.min<-resample(cloud.cover.min,bioclim) #writeRaster(cloud.cover.min,"CloudCoverMin") #res(cloud.cover.mean) #plot(cloud.cover.mean) #res(cloud.cover.max) #plot(cloud.cover.max) #res(cloud.cover.min) #plot(cloud.cover.min) #Enhanced Vegetation Index - Coeficient of Variation: #EVI.cv <- list.files("./Environmental layers/Enhanced Vegetation Index_cv",pattern=".tif", full.names=TRUE) #EVI.cv <- stack(EVI.cv) #EVI.cv <- mask(crop(EVI.cv,neotrop),neotrop) #EVI.cv.10km <- resample(EVI.cv,bioclim) #writeRaster(EVI.cv.10km, "EVIcv10km") #res(EVI.cv.10km) #plot(EVI.cv.10km) #Enhanced Vegetation Index - Range: #EVI.rng <- list.files("./Environmental layers/Enhanced Vegetation Index_range",pattern=".tif", full.names=TRUE) #EVI.rng <- stack(EVI.rng) #EVI.rng <- mask(crop(EVI.rng,neotrop),neotrop) #EVI.rng.10km <- resample(EVI.rng,bioclim) #writeRaster(EVI.rng.10km, "EVIrng10km") #res(EVI.rng.10km) #plot(EVI.rng.10km) #Enhanced Vegetation Index - Standard Deviation: #EVI.std <- list.files("./Environmental layers/Enhanced Vegetation Index_std",pattern=".tif", full.names=TRUE) #EVI.std <- stack(EVI.std) #EVI.std <- mask(crop(EVI.std,neotrop),neotrop) #EVI.std.10km <- resample(EVI.std,bioclim) #writeRaster(EVI.std.10km, "EVIstd10km") #res(EVI.std.10km) #plot(EVI.std.10km) #Forest Coverage: #FOR.cov <- list.files("./Environmental layers/Vegetation coverage/Forest Coverage",pattern=".asc", full.names=TRUE) #FOR.cov <- stack(FOR.cov) #FOR.cov <- mask(crop(FOR.cov,neotrop),neotrop) #writeRaster(FOR.cov, "FORcov") #res(FOR.cov) #plot(FOR.cov) #Grassland/Scrub/Woodland Coverage: #GRASS.cov <- list.files("./Environmental layers/Vegetation coverage/Grassland Coverage",pattern=".asc", full.names=TRUE) #GRASS.cov <- stack(GRASS.cov) #GRASS.cov <- mask(crop(GRASS.cov,neotrop),neotrop) #writeRaster(GRASS.cov, "GRASScov") #res(GRASS.cov) #plot(GRASS.cov) #Water Bodies: #WATB.cov <- list.files("./Environmental layers/Vegetation coverage/Water Bodies",pattern=".asc", full.names=TRUE) #WATB.cov <- stack(WATB.cov) #WATB.cov <- mask(crop(WATB.cov,neotrop),neotrop) #writeRaster(WATB.cov, "WATBcov") #res(WATB.cov) #plot(WATB.cov) #Elevation: #elevation <-list.files("./Environmental layers/Elevation",pattern=".asc", full.names=TRUE) #elevation <-stack(elevation) #elevation <-mask(crop(elevation, neotrop),neotrop) #elevation.10km <-resample(elevation,bioclim) #writeRaster(elevation.10km,"Elevation10km") #res(elevation.10km) #plot(elevation.10km) # Global Relief Model: #relief <- list.files("./Environmental layers/Global Relief Model", pattern="tif", full.names=TRUE) #relief <- stack(relief) #relief <- mask(crop(relief,neotrop),neotrop) #relief.10km <- resample(relief, bioclim) #writeRaster(relief.10km, "relief10km") #res(relief.10km) #plot(relief.10km) #Slope and Aspect: #slope <- terrain(elevation.10km, opt="slope") #writeRaster(slope,"Slope") #res(slope) #plot(slope) #aspect <- terrain(elevation.10km, opt="aspect") #writeRaster(aspect,"Aspect") #res(aspect) #plot(aspect) #Terrain Roughness Index: #roughness <-list.files("./Environmental layers/Terrain Roughness Index",pattern=".tif", full.names=TRUE) #roughness <- stack(roughness) #roughness <-mask(crop(roughness, neotrop),neotrop) #roughness.10km <-resample(roughness,bioclim) #writeRaster(roughness.10km,"Roughness10km") #res(roughness.10km) #plot(roughness.10km) #Topographic Wetness Index: #topowet <-list.files("./Environmental layers/Topographic Wetness Index",pattern=".tif", full.names=TRUE) #topowet <- stack(topowet) #topowet <-mask(crop(topowet, neotrop),neotrop) #topowet.10km <-resample(topowet,bioclim) #writeRaster(topowet.10km,"TopoWet10km") #res(topowet.10km) #plot(topowet.10km) #Potential Evapotranspiration - PET: ### Annual PET: #PET.1km <- raster("./Environmental layers/Potential Evapotranspiration/Global PET - Annual/PET_he_annual/pet_he_yr/w001001.adf") #PET.1km <- mask(crop(PET.1km,neotrop),neotrop) #PET.10km <- resample(PET.1km,bioclim) #writeRaster(PET.10km, "PET10km") #res(PET.10km) #plot(PET.10km) ### PET Coldest Quarter: #PET.cq <- list.files("./Environmental layers/Potential Evapotranspiration/PET Coldest Quarter",pattern=".tif", full.names=TRUE) #PET.cq <- stack(PET.cq) #PET.cq <-mask(crop(PET.cq, neotrop),neotrop) #PET.cq <-resample(PET.cq,bioclim) #writeRaster(PET.cq,"PETcq") #res(PET.cq) #plot(PET.cq) ### PET Driest Quarter: #PET.dq <- list.files("./Environmental layers/Potential Evapotranspiration/PET Driest Quarter",pattern=".tif", full.names=TRUE) #PET.dq <- stack(PET.dq) #PET.dq <-mask(crop(PET.dq, neotrop),neotrop) #PET.dq <-resample(PET.dq,bioclim) #writeRaster(PET.dq,"PETdq") #res(PET.dq) #plot(PET.dq) ### PET Warmest Quarter: #PET.wq <- list.files("./Environmental layers/Potential Evapotranspiration/PET Warmest Quarter",pattern=".tif", full.names=TRUE) #PET.wq <- stack(PET.wq) #PET.wq <- mask(crop(PET.wq, neotrop),neotrop) #PET.wq <- resample(PET.wq,bioclim) #writeRaster(PET.wq,"PETwq") #res(PET.wq) #plot(PET.wq) ### PET Wettest Quarter: #PET.wetq <- list.files("./Environmental layers/Potential Evapotranspiration/PET Wettest Quarter",pattern=".tif", full.names=TRUE) #PET.wetq <- stack(PET.wetq) #PET.wetq <- mask(crop(PET.wetq, neotrop),neotrop) #PET.wetq <- resample(PET.wetq,bioclim) #writeRaster(PET.wetq,"PETwetq") #res(PET.wetq) #plot(PET.wetq) ### PET Seasonality: #PET.seas <- list.files("./Environmental layers/Potential Evapotranspiration/PET Seasonality",pattern=".tif", full.names=TRUE) #PET.seas <- stack(PET.seas) #PET.seas <- mask(crop(PET.seas, neotrop),neotrop) #PET.seas <- resample(PET.seas,bioclim) #writeRaster(PET.seas,"PETseas") #res(PET.seas) #plot(PET.seas) #Aridity Index: #Aridity.1km <- raster("./Environmental layers/Global Aridity and PET database/Global Aridity - Annual/AI_annual/ai_yr/w001001.adf") #Aridity.1km <- mask(crop(Aridity.1km,neotrop),neotrop) #Aridity.10km <- resample(Aridity.1km,bioclim) #writeRaster(Aridity.10km, "Aridity10km") #res(Aridity.10km) #plot(Aridity.10km) #Actual Evapotranspiration: #AET.1km <- raster("./Environmental layers/Global Soil Water Balance and AET/Mean Annual AET/AET_YR/aet_yr/w001001.adf") #AET.1km <- mask(crop(AET.1km,neotrop),neotrop) #AET.10km <- resample(AET.1km,bioclim) #writeRaster(AET.10km, "AET10km") #res(AET.10km) #plot(AET.10km) #Soil Water Stress: #SWS.jan <-raster("./Environmental layers/Global Soil Water Balance and AET/Monthly Soil Water Stress/swc_fr/swc_fr_1/w001001.adf") #SWS.feb <-raster("./Environmental layers/Global Soil Water Balance and AET/Monthly Soil Water Stress/swc_fr/swc_fr_2/w001001.adf") #SWS.mar <-raster("./Environmental layers/Global Soil Water Balance and AET/Monthly Soil Water Stress/swc_fr/swc_fr_3/w001001.adf") #SWS.apr <-raster("./Environmental layers/Global Soil Water Balance and AET/Monthly Soil Water Stress/swc_fr/swc_fr_4/w001001.adf") #SWS.may <-raster("./Environmental layers/Global Soil Water Balance and AET/Monthly Soil Water Stress/swc_fr/swc_fr_5/w001001.adf") #SWS.jun <-raster("./Environmental layers/Global Soil Water Balance and AET/Monthly Soil Water Stress/swc_fr/swc_fr_6/w001001.adf") #SWS.jul <-raster("./Environmental layers/Global Soil Water Balance and AET/Monthly Soil Water Stress/swc_fr/swc_fr_7/w001001.adf") #SWS.aug <-raster("./Environmental layers/Global Soil Water Balance and AET/Monthly Soil Water Stress/swc_fr/swc_fr_8/w001001.adf") #SWS.sep <-raster("./Environmental layers/Global Soil Water Balance and AET/Monthly Soil Water Stress/swc_fr/swc_fr_9/w001001.adf") #SWS.oct <-raster("./Environmental layers/Global Soil Water Balance and AET/Monthly Soil Water Stress/swc_fr/swc_fr_10/w001001.adf") #SWS.nov <-raster("./Environmental layers/Global Soil Water Balance and AET/Monthly Soil Water Stress/swc_fr/swc_fr_11/w001001.adf") #SWS.dec <-raster("./Environmental layers/Global Soil Water Balance and AET/Monthly Soil Water Stress/swc_fr/swc_fr_12/w001001.adf") #SWS.stack <-stack(SWS.jan,SWS.feb,SWS.mar,SWS.apr,SWS.may,SWS.jun,SWS.jul, # SWS.aug,SWS.sep,SWS.oct,SWS.nov,SWS.dec) #SWS.mean.1km <-mean(SWS.stack) #SWS.mean.1km <-mask(crop(SWS.mean.1km,neotrop),neotrop) #SWS.mean.10km <-resample(SWS.mean.1km, bioclim) #writeRaster(SWS.mean.10km,"SWSmean10km") #res(SWS.mean.10km) #plot(SWS.mean.10km) #SWS.max.1km <-max(SWS.stack) #SWS.max.1km <-mask(crop(SWS.max.1km,neotrop),neotrop) #SWS.max.10km <-resample(SWS.max.1km, bioclim) #writeRaster(SWS.max.10km,"SWSmax10km") #res(SWS.max.10km) #plot(SWS.max.10km) #SWS.min.1km <-min(SWS.stack) #SWS.min.1km <-mask(crop(SWS.min.1km,neotrop),neotrop) #SWS.min.10km <-resample(SWS.min.1km, bioclim) #writeRaster(SWS.min.10km,"SWSmin10km") #res(SWS.min.10km) #plot(SWS.min.10km) #Relative Humidity at 3pm: #Humidity.3pm.jan <-raster("./Environmental layers/Relative Humidity at 3 pm/CM10_1975H_Raw_ESRI_RHpm_V1.2/CM10_1975H_Raw_ESRI_RHpm_V1.2/rhpm01/w001001.adf") #Humidity.3pm.feb <-raster("./Environmental layers/Relative Humidity at 3 pm/CM10_1975H_Raw_ESRI_RHpm_V1.2/CM10_1975H_Raw_ESRI_RHpm_V1.2/rhpm02/w001001.adf") #Humidity.3pm.mar <-raster("./Environmental layers/Relative Humidity at 3 pm/CM10_1975H_Raw_ESRI_RHpm_V1.2/CM10_1975H_Raw_ESRI_RHpm_V1.2/rhpm03/w001001.adf") #Humidity.3pm.apr <-raster("./Environmental layers/Relative Humidity at 3 pm/CM10_1975H_Raw_ESRI_RHpm_V1.2/CM10_1975H_Raw_ESRI_RHpm_V1.2/rhpm04/w001001.adf") #Humidity.3pm.may <-raster("./Environmental layers/Relative Humidity at 3 pm/CM10_1975H_Raw_ESRI_RHpm_V1.2/CM10_1975H_Raw_ESRI_RHpm_V1.2/rhpm05/w001001.adf") #Humidity.3pm.jun <-raster("./Environmental layers/Relative Humidity at 3 pm/CM10_1975H_Raw_ESRI_RHpm_V1.2/CM10_1975H_Raw_ESRI_RHpm_V1.2/rhpm06/w001001.adf") #Humidity.3pm.jul <-raster("./Environmental layers/Relative Humidity at 3 pm/CM10_1975H_Raw_ESRI_RHpm_V1.2/CM10_1975H_Raw_ESRI_RHpm_V1.2/rhpm07/w001001.adf") #Humidity.3pm.aug <-raster("./Environmental layers/Relative Humidity at 3 pm/CM10_1975H_Raw_ESRI_RHpm_V1.2/CM10_1975H_Raw_ESRI_RHpm_V1.2/rhpm08/w001001.adf") #Humidity.3pm.sep <-raster("./Environmental layers/Relative Humidity at 3 pm/CM10_1975H_Raw_ESRI_RHpm_V1.2/CM10_1975H_Raw_ESRI_RHpm_V1.2/rhpm09/w001001.adf") #Humidity.3pm.oct <-raster("./Environmental layers/Relative Humidity at 3 pm/CM10_1975H_Raw_ESRI_RHpm_V1.2/CM10_1975H_Raw_ESRI_RHpm_V1.2/rhpm10/w001001.adf") #Humidity.3pm.nov <-raster("./Environmental layers/Relative Humidity at 3 pm/CM10_1975H_Raw_ESRI_RHpm_V1.2/CM10_1975H_Raw_ESRI_RHpm_V1.2/rhpm11/w001001.adf") #Humidity.3pm.dec <-raster("./Environmental layers/Relative Humidity at 3 pm/CM10_1975H_Raw_ESRI_RHpm_V1.2/CM10_1975H_Raw_ESRI_RHpm_V1.2/rhpm12/w001001.adf") #Humidity.3pm.stack <-stack(Humidity.3pm.jan, Humidity.3pm.feb, Humidity.3pm.mar, Humidity.3pm.apr, Humidity.3pm.may, Humidity.3pm.jun, Humidity.3pm.jul, # Humidity.3pm.aug, Humidity.3pm.sep, Humidity.3pm.oct, Humidity.3pm.nov, Humidity.3pm.dec) #Humidity.3pm.mean.20km <-mean(Humidity.3pm.stack) #Humidity.3pm.mean.20km <-mask(crop(Humidity.3pm.mean.20km,neotrop),neotrop) #Humidity.3pm.mean.10km <-resample(Humidity.3pm.mean.20km, bioclim) #writeRaster(Humidity.3pm.mean.10km,"Humidity3pmMean10km") #res(Humidity.3pm.mean.10km) #plot(Humidity.3pm.mean.10km) #Humidity.3pm.max.20km <-max(Humidity.3pm.stack) #Humidity.3pm.max.20km <-mask(crop(Humidity.3pm.max.20km,neotrop),neotrop) #Humidity.3pm.max.10km <-resample(Humidity.3pm.max.20km, bioclim) #writeRaster(Humidity.3pm.max.10km,"Humidity3pmMax10km") #res(Humidity.3pm.max.10km) #plot(Humidity.3pm.max.10km) #Humidity.3pm.min.20km <-min(Humidity.3pm.stack) #Humidity.3pm.min.20km <-mask(crop(Humidity.3pm.min.20km,neotrop),neotrop) #Humidity.3pm.min.10km <-resample(Humidity.3pm.min.20km, bioclim) #writeRaster(Humidity.3pm.min.10km,"Humidity3pmMin10km") #res(Humidity.3pm.min.10km) #plot(Humidity.3pm.min.10km) #Relative Humidity at 9am: #Humidity.9am.jan <-raster("./Environmental layers/Relative Humidity at 9 am/CM10_1975H_Raw_ESRI_RHam_V1.2/CM10_1975H_Raw_ESRI_RHam_V1.2/rham01/w001001.adf") #Humidity.9am.feb <-raster("./Environmental layers/Relative Humidity at 9 am/CM10_1975H_Raw_ESRI_RHam_V1.2/CM10_1975H_Raw_ESRI_RHam_V1.2/rham02/w001001.adf") #Humidity.9am.mar <-raster("./Environmental layers/Relative Humidity at 9 am/CM10_1975H_Raw_ESRI_RHam_V1.2/CM10_1975H_Raw_ESRI_RHam_V1.2/rham03/w001001.adf") #Humidity.9am.apr <-raster("./Environmental layers/Relative Humidity at 9 am/CM10_1975H_Raw_ESRI_RHam_V1.2/CM10_1975H_Raw_ESRI_RHam_V1.2/rham04/w001001.adf") #Humidity.9am.may <-raster("./Environmental layers/Relative Humidity at 9 am/CM10_1975H_Raw_ESRI_RHam_V1.2/CM10_1975H_Raw_ESRI_RHam_V1.2/rham05/w001001.adf") #Humidity.9am.jun <-raster("./Environmental layers/Relative Humidity at 9 am/CM10_1975H_Raw_ESRI_RHam_V1.2/CM10_1975H_Raw_ESRI_RHam_V1.2/rham06/w001001.adf") #Humidity.9am.jul <-raster("./Environmental layers/Relative Humidity at 9 am/CM10_1975H_Raw_ESRI_RHam_V1.2/CM10_1975H_Raw_ESRI_RHam_V1.2/rham07/w001001.adf") #Humidity.9am.aug <-raster("./Environmental layers/Relative Humidity at 9 am/CM10_1975H_Raw_ESRI_RHam_V1.2/CM10_1975H_Raw_ESRI_RHam_V1.2/rham08/w001001.adf") #Humidity.9am.sep <-raster("./Environmental layers/Relative Humidity at 9 am/CM10_1975H_Raw_ESRI_RHam_V1.2/CM10_1975H_Raw_ESRI_RHam_V1.2/rham09/w001001.adf") #Humidity.9am.oct <-raster("./Environmental layers/Relative Humidity at 9 am/CM10_1975H_Raw_ESRI_RHam_V1.2/CM10_1975H_Raw_ESRI_RHam_V1.2/rham10/w001001.adf") #Humidity.9am.nov <-raster("./Environmental layers/Relative Humidity at 9 am/CM10_1975H_Raw_ESRI_RHam_V1.2/CM10_1975H_Raw_ESRI_RHam_V1.2/rham11/w001001.adf") #Humidity.9am.dec <-raster("./Environmental layers/Relative Humidity at 9 am/CM10_1975H_Raw_ESRI_RHam_V1.2/CM10_1975H_Raw_ESRI_RHam_V1.2/rham12/w001001.adf") #Humidity.9am.stack <-stack(Humidity.9am.jan, Humidity.9am.feb, Humidity.9am.mar, Humidity.9am.apr, Humidity.9am.may, Humidity.9am.jun, Humidity.9am.jul, # Humidity.9am.aug, Humidity.9am.sep, Humidity.9am.oct, Humidity.9am.nov, Humidity.9am.dec) #Humidity.9am.mean.20km <-mean(Humidity.9am.stack) #Humidity.9am.mean.20km <-mask(crop(Humidity.9am.mean.20km,neotrop),neotrop) #Humidity.9am.mean.10km <-resample(Humidity.9am.mean.20km, bioclim) #writeRaster(Humidity.9am.mean.10km,"Humidity9amMean10km") #res(Humidity.9am.mean.10km) #plot(Humidity.9am.mean.10km) #Humidity.9am.max.20km <-max(Humidity.9am.stack) #Humidity.9am.max.20km <-mask(crop(Humidity.9am.max.20km,neotrop),neotrop) #Humidity.9am.max.10km <-resample(Humidity.9am.max.20km, bioclim) #writeRaster(Humidity.9am.max.10km,"Humidity9amMax10km") #res(Humidity.9am.max.10km) #plot(Humidity.9am.max.10km) #Humidity.9am.min.20km <-min(Humidity.9am.stack) #Humidity.9am.min.20km <-mask(crop(Humidity.9am.min.20km,neotrop),neotrop) #Humidity.9am.min.10km <-resample(Humidity.9am.min.20km, bioclim) #writeRaster(Humidity.9am.min.10km,"Humidity9amMin10km") #res(Humidity.9am.min.10km) #plot(Humidity.9am.min.10km) ### Soil Grids: # Bulk Density #BulkDensity.0 <- raster("./Environmental layers/Soil Grids/Bulk Density/BLDFIE_M_sl1_250m.tif") #BulkDensity.5 <- raster("./Environmental layers/Soil Grids/Bulk Density/BLDFIE_M_sl2_250m.tif") #BulkDensity.15 <- raster("./Environmental layers/Soil Grids/Bulk Density/BLDFIE_M_sl3_250m.tif") #BulkDensity.30 <- raster("./Environmental layers/Soil Grids/Bulk Density/BLDFIE_M_sl4_250m.tif") #BulkDensity <- stack(BulkDensity.0, BulkDensity.5, BulkDensity.15, BulkDensity.30) #BulkDensity <- mean(BulkDensity) #BulkDensity <- mask(crop(BulkDensity,neotrop),neotrop) #BulkDensity <- resample(BulkDensity,bioclim) #writeRaster(BulkDensity, "BulkDensity.grd") #res(BulkDensity) # Clay Content #Clay.0 <- raster("./Environmental layers/Soil Grids/Clay Content/CLYPPT_M_sl1_250m.tif") #Clay.5 <- raster("./Environmental layers/Soil Grids/Clay Content/CLYPPT_M_sl2_250m.tif") #Clay.15 <- raster("./Environmental layers/Soil Grids/Clay Content/CLYPPT_M_sl3_250m.tif") #Clay.30 <- raster("./Environmental layers/Soil Grids/Clay Content/CLYPPT_M_sl4_250m.tif") #Clay <- stack(Clay.0,Clay.5,Clay.15,Clay.30) #Clay <- mean(Clay) #Clay <- mask(crop(Clay,neotrop),neotrop) #Clay <- resample(Clay, bioclim) #writeRaster(Clay, "Clay.grd") #res(Clay) # Coarse Fragments #Coarse.0 <- raster("./Environmental layers/Soil Grids/Coarse Fragments/CRFVOL_M_sl1_250m.tif") #Coarse.5 <- raster("./Environmental layers/Soil Grids/Coarse Fragments/CRFVOL_M_sl2_250m.tif") #Coarse.15 <- raster("./Environmental layers/Soil Grids/Coarse Fragments/CRFVOL_M_sl3_250m.tif") #Coarse.30 <- raster("./Environmental layers/Soil Grids/Coarse Fragments/CRFVOL_M_sl4_250m.tif") #Coarse <- stack(Coarse.0,Coarse.5,Coarse.15,Coarse.30) #Coarse <- mean(Coarse) #Coarse <- mask(crop(Coarse,neotrop),neotrop) #Coarse <- resample(Coarse, bioclim) #writeRaster(Coarse, "Coarse.grd") #res(Coarse) # Sand Content #Sand.0 <- raster("./Environmental layers/Soil Grids/Sand Content/SNDPPT_M_sl1_250m.tif") #Sand.5 <- raster("./Environmental layers/Soil Grids/Sand Content/SNDPPT_M_sl2_250m.tif") #Sand.15 <- raster("./Environmental layers/Soil Grids/Sand Content/SNDPPT_M_sl3_250m.tif") #Sand.30 <- raster("./Environmental layers/Soil Grids/Sand Content/SNDPPT_M_sl4_250m.tif") #Sand <- stack(Sand.0,Sand.5,Sand.15,Sand.30) #Sand <- mean(Sand) #Sand <- mask(crop(Sand,neotrop),neotrop) #Sand <- resample(Sand, bioclim) #writeRaster(Sand, "Sand.grd") #res(Sand) # Silt Content #Silt.0 <- raster("./Environmental layers/Soil Grids/Silt Content/SLTPPT_M_sl1_250m.tif") #Silt.5 <- raster("./Environmental layers/Soil Grids/Silt Content/SLTPPT_M_sl2_250m.tif") #Silt.15 <- raster("./Environmental layers/Soil Grids/Silt Content/SLTPPT_M_sl3_250m.tif") #Silt.30 <- raster("./Environmental layers/Soil Grids/Silt Content/SLTPPT_M_sl4_250m.tif") #Silt <- stack(Silt.0,Silt.5,Silt.15,Silt.30) #Silt <- mean(Silt) #Silt <- mask(crop(Silt,neotrop),neotrop) #Silt <- resample(Silt, bioclim) #writeRaster(Silt, "Silt.grd") #res(Silt) # Predicted Probability of Occurrence of R horizon #BDRLOG <- raster("./Environmental layers/Soil Grids/BDRLOG/BDRLOG_M_250m.tif") #BDRLOG <- stack(BDRLOG) #BDRLOG <- mask(crop(BDRLOG,neotrop),neotrop) #BDRLOG <- resample(BDRLOG, bioclim) #writeRaster(BDRLOG, "BDRLOG.grd") #res(BDRLOG) # Depth to bedrock up to 200m #BDRICM <- raster("./Environmental layers/Soil Grids/Depth to Bedrock/BDRICM_M_250m.tif") #BDRICM <- stack(BDRICM) #BDRICM <- mask(crop(BDRICM,neotrop),neotrop) #BDRICM <- resample(BDRICM, bioclim) #writeRaster(BDRICM, "BDRICM.grd") #res(BDRICM) # Soil organic carbon stock #CARBON.0 <- raster("./Environmental layers/Soil Grids/Carbon stock/OCSTHA_M_sd1_250m.tif") #CARBON.5 <- raster("./Environmental layers/Soil Grids/Carbon stock/OCSTHA_M_sd2_250m.tif") #CARBON.15 <- raster("./Environmental layers/Soil Grids/Carbon stock/OCSTHA_M_sd3_250m.tif") #CARBON.30 <- raster("./Environmental layers/Soil Grids/Carbon stock/OCSTHA_M_sd4_250m.tif") #CARBON <- stack(CARBON.0, CARBON.5, CARBON.15, CARBON.30) #CARBON <- mean (CARBON) #CARBON <- mask(crop(CARBON,neotrop),neotrop) #CARBON <- resample(CARBON, bioclim) #writeRaster(CARBON, "CARBON.grd") #res(CARBON) # pH in H20 #pH_w.0 <- raster("./Environmental layers/Soil Grids/PHIHOX/PHIHOX_M_sl1_250m.tif") #pH_w.5 <- raster("./Environmental layers/Soil Grids/PHIHOX/PHIHOX_M_sl2_250m.tif") #pH_w.15 <- raster("./Environmental layers/Soil Grids/PHIHOX/PHIHOX_M_sl3_250m.tif") #pH_w.30 <- raster("./Environmental layers/Soil Grids/PHIHOX/PHIHOX_M_sl4_250m.tif") #pH_w <- stack(pH_w.0,pH_w.5,pH_w.15,pH_w.30) #pH_w <- mean (pH_w) #pH_w <- mask(crop(pH_w,neotrop),neotrop) #pH_w <- resample(pH_w, bioclim) #writeRaster(pH_w, "pH_w.grd") #res(pH_w) # pH in KCl #pH_k.0 <- raster("./Environmental layers/Soil Grids/PHIKCL/PHIKCL_M_sl1_250m.tif") #pH_k.5 <- raster("./Environmental layers/Soil Grids/PHIKCL/PHIKCL_M_sl2_250m.tif") #pH_k.15 <- raster("./Environmental layers/Soil Grids/PHIKCL/PHIKCL_M_sl3_250m.tif") #pH_k.30 <- raster("./Environmental layers/Soil Grids/PHIKCL/PHIKCL_M_sl4_250m.tif") #pH_k <- stack(pH_k.0,pH_k.5,pH_k.15,pH_k.30) #pH_k <- mean (pH_k) #pH_k <- mask(crop(pH_k,neotrop),neotrop) #pH_k <- resample(pH_k, bioclim) #writeRaster(pH_k, "pH_k.grd", overwrite=TRUE) #res(pH_k) #ORCDRC #ORCDRC.0 <- raster("./Environmental layers/Soil Grids/ORCDRC/ORCDRC_M_sl1_250m.tif") #ORCDRC.5 <- raster("./Environmental layers/Soil Grids/ORCDRC/ORCDRC_M_sl2_250m.tif") #ORCDRC.15 <- raster("./Environmental layers/Soil Grids/ORCDRC/ORCDRC_M_sl3_250m.tif") #ORCDRC.30 <- raster("./Environmental layers/Soil Grids/ORCDRC/ORCDRC_M_sl4_250m.tif") #ORC <- stack(ORCDRC.0,ORCDRC.5,ORCDRC.15,ORCDRC.30) #ORC <- mean (ORC) #ORC <- mask(crop(ORC,neotrop),neotrop) #ORC <- resample(ORC, bioclim) #writeRaster(ORC, "ORC.grd") #res(ORC) # CEC #CEC.0 <- raster("./Environmental layers/Soil Grids/CECSOL/CECSOL_M_sl1_250m.tif") #CEC.5 <- raster("./Environmental layers/Soil Grids/CECSOL/CECSOL_M_sl2_250m.tif") #CEC.15 <- raster("./Environmental layers/Soil Grids/CECSOL/CECSOL_M_sl3_250m.tif") #CEC.30 <- raster("./Environmental layers/Soil Grids/CECSOL/CECSOL_M_sl4_250m.tif") #CEC <- stack(CEC.0,CEC.5,CEC.15,CEC.30) #CEC <- mean (CEC) #CEC <- mask(crop(CEC,neotrop),neotrop) #CEC <- resample(CEC, bioclim) #writeRaster(CEC, "CEC.grd") #res(CEC) #--------------------------------------------------------------------------------------------# ### IF YOU HAVE ALREADY DOWNLOAD AND TREATED ALL LAYERS, YOU SHOULD CONTINUE FROM HERE ###### #------------------------------------------------------------------------------------------# #-----------------------------------# # Loading environmental layers ##### #-----------------------------------# bioclim <- list.files("./Environmental layers/CHELSA", pattern="grd", full.names=TRUE) bioclim <- stack(bioclim) solar.radiation.mean <-raster("./Environmental layers/Solar Radiation/SolarRadiationMean.grd") names(solar.radiation.mean) = "Solar Rad_Mean" solar.radiation.max <-raster("./Environmental layers/Solar Radiation/SolarRadiationMax.grd") names(solar.radiation.max) = "Solar Rad_Max" solar.radiation.min <-raster("./Environmental layers/Solar Radiation/SolarRadiationMin.grd") names(solar.radiation.min) = "Solar Rad_Min" water.vapor.pressure.mean<-raster("./Environmental layers/Water Vapor Pressure/WaterVaporPressureMean.grd") names(water.vapor.pressure.mean) = "Water Vapor Press_Mean" water.vapor.pressure.max <-raster("./Environmental layers/Water Vapor Pressure/WaterVaporPressureMax.grd") names(water.vapor.pressure.max) = "Water Vapor Press_Max" water.vapor.pressure.min <-raster("./Environmental layers/Water Vapor Pressure/WaterVaporPressureMin.grd") names(water.vapor.pressure.min) = "Water Vapor Press_Min" wind.speed.mean <-raster("./Environmental layers/Wind Speed/WindSpeedMean.grd") names(wind.speed.mean) = "Wind Speed_Mean" wind.speed.max <-raster("./Environmental layers/Wind Speed/WindSpeedMax.grd") names(wind.speed.max) = "Wind Speed_Max" wind.speed.min <-raster("./Environmental layers/Wind Speed/WindSpeedMin.grd") names(wind.speed.min) = "Wind Speed_Min" cloud.cover.mean <-raster("./Environmental layers/Cloud Cover/CloudCoverMean.grd") names(cloud.cover.mean) = "Cloud Cover_Mean" cloud.cover.max <- raster("./Environmental layers/Cloud Cover/CloudCoverMax.grd") names(cloud.cover.max) = "Cloud Cover_Max" cloud.cover.min <- raster("./Environmental layers/Cloud Cover/CloudCoverMin.grd") names(cloud.cover.min) = "Cloud Cover_Min" EVI.cv.10km <- raster("./Environmental layers/Enhanced Vegetation Index_cv/EVIcv10km.grd") names(EVI.cv.10km) = "EVI_cv" EVI.rng.10km <- raster("./Environmental layers/Enhanced Vegetation Index_rng/EVIrng10km.grd") names(EVI.rng.10km) = "EVI_rng" EVI.std.10km <- raster("./Environmental layers/Enhanced Vegetation Index_std/EVIstd10km.grd") names(EVI.std.10km) = "EVI_std" FOR.cov <- raster("./Environmental layers/Vegetation coverage/Forest coverage/FORcov.grd") names(FOR.cov) = "FOREST_cov" GRASS.cov <- raster("./Environmental layers/Vegetation coverage/Grassland coverage/GRASScov.grd") names(GRASS.cov) = "GRASS_cov" WATB.cov <- raster("./Environmental layers/Vegetation coverage/Water Bodies/WATBcov.grd") names(WATB.cov) = "WATBODIES_cov" elevation.10km <- raster("./Environmental layers/Elevation/Elevation10km.grd") names(elevation.10km) = "Elevation" slope <-raster("./Environmental layers/Slope/Slope.grd") names(slope) = "Slope" aspect <-raster("./Environmental layers/Aspect/Aspect.grd") names(aspect) = "Aspect" roughness.10km <- raster("./Environmental layers/Terrain Roughness Index/Roughness10km.grd") names(roughness.10km) = "Roughness" topowet.10km <- raster("./Environmental layers/Topographic Wetness Index/TopoWet10km.grd") names(topowet.10km) = "TopoWet" PET.10km <- raster("./Environmental layers/Potential Evapotranspiration/Global PET - Annual/PET10km.grd") names(PET.10km) = "Annual PET" PET.cq <- raster("./Environmental layers/Potential Evapotranspiration/PET Coldest Quarter/PETcq.grd") names(PET.cq) = "PET_ColdQuart" PET.dq <- raster("./Environmental layers/Potential Evapotranspiration/PET Driest Quarter/PETdq.grd") names(PET.dq) = "PET_DriQuart" PET.wq <- raster("./Environmental layers/Potential Evapotranspiration/PET Warmest Quarter/PETwq.grd") names(PET.wq) = "PET_WarmQuart" PET.wetq <-raster("./Environmental layers/Potential Evapotranspiration/PET Wettest Quarter/PETwetq.grd") names(PET.wetq) = "PET_WetQuart" PET.seas <-raster("./Environmental layers/Potential Evapotranspiration/PET Seasonality/PETseas.grd") names(PET.seas) = "PET_Seas" Aridity.10km <-raster("./Environmental layers/Global Aridity/Global Aridity - Annual/Aridity10km") names(Aridity.10km) = "Aridity" AET.10km <-raster("./Environmental layers/Actual Evapotranspiration/Mean Annual AET/AET10km.grd") names(AET.10km) = "AET" SWS.mean.10km <-raster("./Environmental layers/Soil Water Stress/Monthly Soil Water Stress/SWSmean10km.grd") names(SWS.mean.10km) = "SWS_mean" SWS.max.10km <-raster("./Environmental layers/Soil Water Stress/Monthly Soil Water Stress/SWSmax10km.grd") names(SWS.max.10km) = "SWS_max" SWS.min.10km <-raster("./Environmental layers/Soil Water Stress/Monthly Soil Water Stress/SWSmin10km.grd") names(SWS.min.10km) = "SWS_min" relief.10km <-raster("./Environmental layers/Global Relief Model/relief10km.grd") names(relief.10km) = "Relief" Humidity.3pm.mean.10km <-raster("./Environmental layers/Relative Humidity 3pm/Humidity3pmMean10km.grd") names(Humidity.3pm.mean.10km) = "Humidity3pm_mean" Humidity.3pm.min.10km <-raster("./Environmental layers/Relative Humidity 3pm/Humidity3pmMin10km.grd") names(Humidity.3pm.min.10km) = "Humidity3pm_min" Humidity.3pm.max.10km <-raster("./Environmental layers/Relative Humidity 3pm/Humidity3pmMax10km.grd") names(Humidity.3pm.max.10km) = "Humidity3pm_max" Humidity.9am.mean.10km <-raster("./Environmental layers/Relative Humidity 9am/Humidity9amMean10km.grd") names(Humidity.9am.mean.10km) = "Humidity9am_mean" Humidity.9am.max.10km <-raster("./Environmental layers/Relative Humidity 9am/Humidity9amMax10km.grd") names(Humidity.9am.max.10km) = "Humidity9am_max" Humidity.9am.min.10km <-raster("./Environmental layers/Relative Humidity 9am/Humidity9amMin10km.grd") names(Humidity.9am.min.10km) = "Humidity9am_min" BulkDensity <- raster("./Environmental layers/Soil Grids/Bulk Density/BulkDensity.grd") names(BulkDensity) = "BulkDensity" Clay <- raster("./Environmental layers/Soil Grids/Clay Content/Clay.grd") names(Clay) = "Clay" Coarse <- raster("./Environmental layers/Soil Grids/Coarse Fragments/Coarse.grd") names(Coarse) = "Coarse" Sand <- raster("./Environmental layers/Soil Grids/Sand Content/Sand.grd") names(Sand) = "Sand" Silt <- raster("./Environmental layers/Soil Grids/Silt Content/Silt.grd") names(Silt) = "Silt" BDRLOG <- raster("./Environmental layers/Soil Grids/BDRLOG/BDRLOG.grd") names(BDRLOG) = "BDRLOG" BDRICM <- raster("./Environmental layers/Soil Grids/Depth to Bedrock/BDRICM.grd") names(BDRICM) = "BDRICM" CARBON <- raster("./Environmental layers/Soil Grids/Carbon stock/CARBON.grd") names(CARBON) = "CARBON" pH_H20 <- raster("./Environmental layers/Soil Grids/PHIHOX/pH_w.grd") names(pH_H20) = "pH_H20" CEC <- raster("./Environmental layers/Soil Grids/CECSOL/CEC.grd") names(CEC) = "CEC" #------------------------------------------------------------------------# ############### Stacking all environmental layers ####################### #----------------------------------------------------------------------# # If you wish to use the layers from WorldClim 2.0 instead of the layers # from CHELSA, you should replace bioclim by bio.wc below. bio.crop <- stack(bioclim, solar.radiation.mean, solar.radiation.max, solar.radiation.min, water.vapor.pressure.mean, water.vapor.pressure.max, water.vapor.pressure.min, wind.speed.mean, wind.speed.max, wind.speed.min, cloud.cover.mean, cloud.cover.max, cloud.cover.min, EVI.cv.10km, EVI.rng.10km, EVI.std.10km, FOR.cov, GRASS.cov, WATB.cov, elevation.10km, relief.10km, slope, aspect, roughness.10km, topowet.10km, PET.10km, PET.cq, PET.dq, PET.wq, PET.wetq, PET.seas, Aridity.10km, AET.10km, SWS.mean.10km, SWS.min.10km, SWS.max.10km, Humidity.3pm.mean.10km, Humidity.3pm.min.10km, Humidity.3pm.max.10km, Humidity.9am.mean.10km, Humidity.9am.max.10km, Humidity.9am.min.10km, BulkDensity, Clay, Coarse, Sand, Silt, BDRLOG, BDRICM, CARBON, pH_H20,CEC) bio.crop res(bio.crop) ##0.083 = aprox. 10km #----------------------------------------------------------------# ##################### PCA ####################################### #--------------------------------------------------------------# #install.packages("FactoMineR") #library(FactoMineR) #bio.crop.df<-as.data.frame(bio.crop) #PCA<-PCA(bio.crop.df) memory.limit(1000000) env.selected1 <- rasterPCA(bio.crop, nComp=13,scores = TRUE, cor=TRUE, spca = TRUE, bylayer=TRUE, filename="PCA.grd", overwrite=TRUE) # Here I selected the first 13 components because they account for more than 90% # of the total variance considering the 70 predictors of this routine for the # entire Neotropical Region (10-km resolution). #env.selected1$model$loadings #write.table(env.selected1$model$loadings, 'cont.csv', sep = ',') summary(env.selected1$model) #to verify the explanation of each PCA component env.selected <-stack(env.selected1$map) env.selected res(env.selected) plot(env.selected) names(env.selected) #---------------------------------------# ### Loading species occurrence data #### #-------------------------------------# #The species matrix should be exactly as demonstrated below: #sp lon lat #Genera.species1 -000.00 -000.00 #Genera.species1 -000.00 -000.00 #Genera.species1 -000.00 -000.00 #Don't forget the '.' between genera and species' epithet #The same name for the same species #negative coordinates for South Hemisphere #positive coordinates for North Hemisphere spp<-read.table(file.choose(),header=T,sep=",") dim(spp) View(spp) #If you would like to obtain values of the 70 environmental predictors #for each of your occurrence records: spp1<-spp[,-1] View(spp1) ext<-extract(bio.crop,spp1) ext<-cbind(spp,ext) View(ext) write.table(ext,"Variables for each site.csv") # Visualizing species occurrence records on a map # data(wrld_simpl) plot(wrld_simpl, xlim=c(-85, -35), ylim=c(-55, 15), col="lightgray", axes=TRUE) points(spp$lon, spp$lat, col="black", bg="red", pch=21, cex=1.0, lwd=1.0) # Formating occurences data table(spp$sp) #The second code (after '$') needs to match the code entered in the matrix spp especies <- unique(spp$sp) #ditto especies # Creating objects for models calibration models1<-c("CTA","RF", "GBM") models2<-c("MAXENT.Phillips", "GLM", "GAM", "MARS","ANN", "FDA") n.runs = 2 # number of RUNs (use at least 10) n.algo1 = length(models1)# number of algorithms n.algo2 = length(models2) #numero de algorithms n.conj.pa2 = 2 # set of pseudo-absences (use at least 10) env.selected = bio.crop especie = especies[1] # To model without a loop, remove the '#' of this line and add it to the 'for', 'foreach' and '.packages' #-------------------------# #beginning of the loop#### #-----------------------# # for(especie in especies[1:length(especies)]){ # foreach(especie = especies, # For parallel looping (Multiple Species) # .packages = c("raster", "biomod2", 'sp', "sdmvspecies", "filesstrings")) %dopar% { # ini1 = Sys.time() # criando tabela para uma especie occs <- spp[spp$sp == especie, c("lon", "lat")] # nome = strsplit(as.vector(especie), " ") # especie = paste(nome[[1]][1], nome[[1]][2], sep = ".") # Selecionado pontos espacialmente únicos # mask <- env.selected[[1]] {(cell <- cellFromXY(mask, occs[, 1:2])) # get the cell number for each point (x<-(cbind(occs[, 1:2], cell))) #dup <- duplicated(cbind(occs[, 1:2], cell)) (dup2 <- duplicated(cbind(cell))) xv<-data.frame(x,dup2) xv[xv=="TRUE"]<-NA (xv<-na.omit(xv)) xv<-xv[,1:2] occs =xv # select the records that are not duplicated } occs #pontos espacialmente únicos dim(occs) #-----------------------------------------------# # GENERATING OTHER REQUIRED OBJECTS FOR SDM #### #---------------------------------------------# # Convert dataset to SpatialPointsDataFrame (only presences) myRespXY <- occs[, c("lon", "lat")] #Caso dê algum erro aqui, veja como você intitulou as colunas da sua matriz. # Creating occurrence data object occurrence.resp <- rep(1, length(myRespXY$lon)) #------------------------------------------# # FIT SPECIES DISTRIBUTION MODELS - SDMS #### #----------------------------------------# try({ coord1 = occs sp::coordinates(coord1) <- ~ lon + lat raster::crs(coord1) <- raster::crs(env.selected) dist.mean <- mean(sp::spDists( x = coord1, longlat = T, segments = FALSE )) dist.min = 5 dist.min <- min(sp::spDists(x = coord1, longlat = T, segments = F)) dist.min = 5 write.table( c(dist.min, dist.mean), paste0('./outputs/', especie,"_", ".csv"), row.names = F, sep = "," ) }) dim(occs) PA.number <- length(occs[, 1]) PA.number #número de pontos de ocorrência espacialmente únicos diretorio = paste0("Occurrence.", especie) ##### FORMATING DATA ##### # Preparando para CTA, GBM e RF: sppBiomodData.PA.equal <- BIOMOD_FormatingData( resp.var = occurrence.resp, expl.var = env.selected, resp.xy = myRespXY, resp.name = diretorio, PA.nb.rep = n.conj.pa2, #numero de datasets de pseudoausencias PA.nb.absences = PA.number, #= numero de pseudoausencias = numero de pontos espacialmente unicos PA.strategy = "disk", # PA.sre.quant = 0.10, PA.dist.min = dist.min * 1000, PA.dist.max = dist.mean * 1000, na.rm = TRUE ) sppBiomodData.PA.equal #Preparando para os demais algoritmos: sppBiomodData.PA.10000 <- BIOMOD_FormatingData( resp.var = occurrence.resp, expl.var = env.selected, resp.xy = myRespXY, resp.name = diretorio, PA.nb.rep = n.conj.pa2, PA.nb.absences = 1000, PA.strategy = "disk", # PA.sre.quant = 0.10, PA.dist.min = dist.min * 1000, PA.dist.max = dist.mean * 1000, na.rm = TRUE ) sppBiomodData.PA.10000 #Alocar o Maxent no diretorio correto (certifique-se que o java esteja instalado e atualizado) #MaxEnt .jar jar <- paste0(system.file(package = "dismo"), "/java/maxent.jar") if (file.exists(jar) != T) { url = "http://biodiversityinformatics.amnh.org/open_source/maxent/maxent.php?op=download" download.file(url, dest = "maxent.zip", mode = "wb") unzip("maxent.zip", files = "maxent.jar", exdir = system.file("java", package = "dismo")) unlink("maxent.zip") warning("Maxent foi colocado no diret?rio") } system.file("java", package = "dismo") myBiomodOption <- BIOMOD_ModelingOptions(MAXENT.Phillips = list(path_to_maxent.jar = jar)) # save.image() #---------------# # Modeling #### #-------------# # Com partição treino x teste: sppModelOut.PA.equal <- BIOMOD_Modeling( sppBiomodData.PA.equal, models =models1, models.options = NULL, NbRunEval = n.runs, #número de repeticoes para cada algoritmo DataSplit = 70,#percentagem de pts para treino. Prevalence = 0.5, VarImport = 0,#caso queira avaliar a importancia das variaveis, mudar para 10 ou 100 permutacoes models.eval.meth = c("TSS", "ROC"), SaveObj = TRUE, rescal.all.models = TRUE, do.full.models = FALSE, modeling.id = "spp_presente" ) # import.var.equal<-data.frame(sppModelOut.PA.equal@variables.importances@val) # names(import.var.equal)<-rep(c('GBM','CTA','RF'),n.runs + n.conj.pa2) # import.var.equal # write.table(import.var.equal, # paste0("./outputs/", especie, "_", "Var.import.PA.equal.csv"), sep = ',') sppModelOut.PA.10000 <- BIOMOD_Modeling( sppBiomodData.PA.10000, models = models2, models.options = myBiomodOption, NbRunEval = n.runs, #número de repetições para cada algoritmo DataSplit = 70, #percentagem de pts para treino. Prevalence = 0.5, VarImport = 0, #caso queira avaliar a importancia das variaveis, mudar para 10 ou 100 permutacoes models.eval.meth = c("TSS", "ROC"), SaveObj = TRUE, rescal.all.models = TRUE, do.full.models = FALSE, modeling.id = "spp_presente" ) # import.var.1000<-data.frame(sppModelOut.PA.10000@variables.importances@val) # names(import.var.1000)<-rep(c("MAXENT.Phillips", "GLM", "GAM", "ANN", "FDA", "MARS"),n.runs + n.conj.pa2) # import.var.1000 # write.table(import.var.1000, # paste0("./outputs/", especie, "_", "Var.import.PA.1000.csv"), sep = ',') #---------------------------------# # EVALUATE MODELS USING BIOMOD2 ## #-------------------------------# # Sobre as metricas avaliativas, # ver http://www.cawcr.gov.au/projects/verification/#Methods_for_dichotomous_forecasts ##### Evaluation of Models #### sppModelEval.PA.equal <- get_evaluations(sppModelOut.PA.equal)#GBM, CTA e RF sppModelEval.PA.equal write.table( sppModelEval.PA.equal, paste0("./outputs/", especie, "_", "EvaluationsAll_1.csv") ) sppModelEval.PA.10000 <- get_evaluations(sppModelOut.PA.10000) #Os demais. sppModelEval.PA.10000 write.table( sppModelEval.PA.10000, paste0("./outputs/", especie, "_", "EvaluationsAll_2.csv") ) # Sumarizando as métricas avaliativas sdm.models1 <-models1 sdm.models1 eval.methods1 <- c("TSS", "ROC") #2 evaluation methods eval.methods1 ##### Eval.1 #### means.i1 <- numeric(0) for (i in 1:n.algo1) { m1 <- sppModelEval.PA.equal[paste(eval.methods1[1]), "Testing.data", paste(sdm.models1[i]), ,] means.i1 = c(means.i1, m1) } summary.eval.equal <- data.frame(rep(sdm.models1, each = n.runs*n.conj.pa2), rep(1:n.conj.pa2, each = n.runs), rep(1:n.runs, n.algo1), means.i1) names(summary.eval.equal) <- c("Model", "PA","Run", "TSS") summary.eval.equal write.table( summary.eval.equal, paste0("./outputs/", especie, "_", "Models1_Evaluation.csv") ) #----------------------------------------------------------------------------------------# means.i1 <- numeric(0) for (i in 1:n.algo1) { m1 <- sppModelEval.PA.equal[paste(eval.methods1[2]), "Sensitivity", paste(sdm.models1[i]), ,] means.i1 = c(means.i1, m1) } summary.eval.equal.1 <- data.frame(means.i1) summary.eval.equal.1 (test1<-cbind(summary.eval.equal,summary.eval.equal.1)) names(test1)<-c("Model", "PA","Run","TSS","Se") test1 #----------------------------------------------------------------------------------------# means.i1.1 <- numeric(0) means.j1.1 <- numeric(2) for (i in 1:n.algo1){ for (j in 1:2){ means.j1.1[j] <- mean(sppModelEval.PA.equal[paste(eval.methods1[j]),"Testing.data",paste(sdm.models1[i]),,]) } means.i1.1 <- c(means.i1.1, means.j1.1) } summary.eval.equal.mean <- data.frame(rep(sdm.models1,each=j), rep(eval.methods1,i), means.i1.1) names(summary.eval.equal.mean) <- c("Model", "Method", "Mean") summary.eval.equal.mean write.table(summary.eval.equal.mean, paste0("./outputs/", especie, "_", "Models1_Evaluation_Mean.csv")) sd.i1 <- numeric(0) sd.j1 <- numeric(2) for (i in 1:n.algo1) { for (j in 1:2) { sd.j1[j] <- sd(sppModelEval.PA.equal[paste(eval.methods1[j]), "Testing.data", paste(sdm.models1[i]), ,]) } sd.i1 <- c(sd.i1, sd.j1) } summary.eval.equal.sd <- data.frame(rep(sdm.models1, each = 2), rep(eval.methods1, n.algo1), sd.i1) names(summary.eval.equal.sd) <- c("Model", "Method", "SD") summary.eval.equal.sd write.table( summary.eval.equal.sd, paste0("./outputs/", especie, "_", "Models1_Evaluation_SD.csv") ) sdm.models2 <-models2 #7 models sdm.models2 eval.methods2 <- c("TSS", "ROC") #2 evaluation methods eval.methods2 ##### Eval.2 #### means.i2 <- numeric(0) for (i2 in 1:n.algo2) { m2 <- sppModelEval.PA.10000[paste(eval.methods2[1]), "Testing.data", paste(sdm.models2[i2]), ,] means.i2 = c(means.i2, m2) } summary.eval.10000 <- data.frame(rep(sdm.models2, each = n.runs*n.conj.pa2), rep(1:n.conj.pa2, each = n.runs), rep(1:n.runs, n.algo2), means.i2) names(summary.eval.10000) <- c("Model", "PA","Run", "TSS") summary.eval.10000 write.table( summary.eval.10000, paste0("./outputs/", especie, "_", "Models2_Evaluation.csv") ) #----------------------------------------------------------------------------------------# means.i21 <- numeric(0) for (i21 in 1:n.algo2) { m21 <- sppModelEval.PA.10000[paste(eval.methods2[2]), "Sensitivity", paste(sdm.models2[i21]), ,] means.i21 = c(means.i21, m21) } summary.eval.10000.1 <- data.frame(means.i21) summary.eval.10000.1 (test2<-cbind(summary.eval.10000,summary.eval.10000.1)) names(test2)<-c("Model", "PA","Run","TSS","Se") test2 #----------------------------------------------------------------------------------------# means.i2.2 <- numeric(0) means.j2.2 <- numeric(2) for (i in 1:n.algo2){ for (j in 1:2){ means.j2.2[j] <- mean(sppModelEval.PA.10000[paste(eval.methods2[j]),"Testing.data",paste(sdm.models2[i]),,], na.rm = T) } means.i2.2 <- c(means.i2.2, means.j2.2) } summary.eval.10000.mean <- data.frame(rep(sdm.models2,each=j), rep(eval.methods2,i), means.i2.2) names(summary.eval.10000.mean) <- c("Model", "Method", "Mean") summary.eval.10000.mean write.table(summary.eval.10000.mean, paste0("./outputs/", especie, "_", "Models2_Evaluation_Mean.csv")) sd.i2 <- numeric(0) sd.j2 <- numeric(2) for (i in 1:n.algo2) { for (j in 1:2) { sd.j2[j] <- sd(sppModelEval.PA.10000[paste(eval.methods2[j]), "Testing.data", paste(sdm.models2[i]), ,]) } sd.i2 <- c(sd.i2, sd.j2) } summary.eval.10000.sd <- data.frame(rep(sdm.models2, each = 2), rep(eval.methods2, n.algo2), sd.i2) names(summary.eval.10000.sd) <- c("Model", "Method", "SD") summary.eval.10000.sd write.table( summary.eval.10000.sd, paste0("./outputs/", especie, "_", "Models2_Evaluation_SD.csv") ) #-----------------------------# # BUILDING OF PROJECTIONS #### #---------------------------# spp.projections_1 <- BIOMOD_Projection( modeling.output = sppModelOut.PA.equal, new.env = env.selected, proj.name = "Cur1_presente", selected.models = "all", #binary.meth = "ROC", output.format = ".grd" ) spp.projections_2 <- BIOMOD_Projection( modeling.output = sppModelOut.PA.10000, new.env = env.selected, proj.name = "Cur2_presente", selected.models = "all", #binary.meth = "ROC", output.format = ".grd" ) # save.image() ### Definir diretório onde está o arquivo proj_Cur1_presente_Occurrence.grd projections_1 <- stack( paste0( "./", diretorio, "/proj_Cur1_presente/proj_Cur1_presente_Occurrence.", especie, ".grd" ) ) names(projections_1) summary.eval.equal_1<-test1 x1<-length(na.omit(summary.eval.equal_1$TSS)) summary.eval.equal_1 <-na.omit(summary.eval.equal_1) summary.eval.equal_1 = summary.eval.equal_1[order(summary.eval.equal_1$Run),] summary.eval.equal_1 = summary.eval.equal_1[order(summary.eval.equal_1$PA),] summary.eval.equal_1$ID = 1:x1 sel = summary.eval.equal_1[summary.eval.equal_1[, "TSS"] > 0.400,] sel <- na.omit(sel) projections.1 = (subset(projections_1, sel[, "ID"])) proj.select1 <- names(projections.1) ### Definir diretório onde está o arquivo proj_Cur2_presente_Occurrence.grd projections_2 <- stack( paste0( "./", diretorio, "/proj_Cur2_presente/proj_Cur2_presente_Occurrence.", especie, ".grd" ) ) names(projections_2) summary.eval.10000_1<-test2 x2<-length(na.omit(summary.eval.10000_1$TSS)) summary.eval.10000_1 <-na.omit(summary.eval.10000_1) summary.eval.10000_1 = summary.eval.10000_1[order(summary.eval.10000_1$Run),] summary.eval.10000_1 = summary.eval.10000_1[order(summary.eval.10000_1$PA),] summary.eval.10000_1$ID = 1:x2 sel2 = summary.eval.10000_1[summary.eval.10000_1[, "TSS"] > 0.400,] sel2 <- na.omit(sel2) projections.2 = (subset(projections_2, sel2[, "ID"])) proj.select2 <- names(projections.2) #-----------------------------------------------# # Mean of the models by algorithm (Present) #### #---------------------------------------------# projections.all1 <- stack(projections.1) projections.all2 <- stack(projections.2) #--------------------------------# # Ensemble - Current Climate #### #------------------------------# all.pres<-stack(projections.1, projections.2) # Regression RG<-c("GLM", "GAM", "FDA", "MARS") fam.reg<-stack() for (l in 1:length(RG)) { fam.reg<- stack(fam.reg, subset(all.pres, grep(RG[l], names(all.pres)))) } fam.reg fam.reg.m<-mean(fam.reg) writeRaster( fam.reg.m, filename = paste0("./outputs/", especie, "_", "Regression - Current Climate.tif"), format = "GTiff", overwrite = TRUE ) # Machine Learning MC<-c("MAXENT.Phillips", "RF", "ANN","GBM", "CTA") fam.mac<-stack() for (l in 1:length(MC)) { fam.mac<- stack(fam.mac, subset(all.pres, grep(MC[l], names(all.pres)))) } fam.mac fam.mac.m<-(mean(fam.mac)) writeRaster( fam.mac.m, filename = paste0("./outputs/", especie, "_", "Machine - Current Climate.tif"), format = "GTiff", overwrite = TRUE ) # All # try({ projections.all.mean <- mean(fam.reg.m,fam.mac.m) / 1000 writeRaster( projections.all.mean, filename = paste0("./outputs/", especie, "_", "Ensemble - Current Climate.tif"), format = "GTiff", overwrite = TRUE ) # }) #--------------------------# # Scores ROC Threshold #### #------------------------# scores_ROC_equal<-subset(sel, select = c(Model, Se)) scores_ROC_equal[scores_ROC_equal=='-Inf']<-NA scores_ROC_equal[scores_ROC_equal=='Inf']<-NA scores_ROC_equal<-na.omit(scores_ROC_equal) write.table(scores_ROC_equal, paste0("./outputs/",especie, "_", "scores_equal_.csv")) ## Evaluation Scores of the Projections with PA.10000 scores_ROC_10000<-subset(sel2, select = c(Model, Se)) scores_ROC_10000[scores_ROC_10000=='-Inf']<-NA scores_ROC_10000[scores_ROC_10000=='Inf']<-NA scores_ROC_10000<-na.omit(scores_ROC_10000) write.table(scores_ROC_10000, paste0("./outputs/",especie, "_", "scores_10000_.csv")) #Scores mean t<-rbind(scores_ROC_equal, scores_ROC_10000) (score.1<-mean(sel$Se)) (score.2<-mean(sel2$Se)) (score.all<-(mean(cbind(score.1,score.2)/100))) # write.table(th_mean, paste0("./outputs/",especie, "_", "scores_mean.csv")) # Regression fam.reg.d<-NULL for (l in 1:length(RG)) { fam.reg.d<- rbind(fam.reg.d, subset(t, Model== RG[l], select = c(Model, Se))) } fam.reg.d.m<-mean(fam.reg.d$Se) # Machine Learning fam.mac.d<-NULL for (l in 1:length(MC)) { fam.mac.d<- rbind(fam.mac.d, subset(t, Model== MC[l], select = c(Model, Se))) } fam.mac.d.m<-mean(fam.mac.d$Se) # score mean (s.m<-mean(fam.reg.d.m,fam.mac.d.m)/100) #-------------------------------------------------------# # Binary models by each algorithm (Current Climate) #### #-----------------------------------------------------# {th<- function(x,y){ if("RasterLayer" %in% class(x)){ v<-as.data.frame(x, h=T,xy=F) v[v=='0']<-NA v.l<-na.omit(v) (vlen<-length(v.l)) n<-raster::ncell(x) (PR<-vlen/n) # PR }else{ cat("x need be raste layer object") } if("numeric" %in% class(y)){ (Se<-y) #Sencitivity 0 to 1 (VDl <- Se-PR) }else stop( # VDI cat("y need be numeric object")) PA <- convertToPA( x, PA.method = "probability", prob.method = "logistic", beta = VDl, alpha = -0.05, plot = T ) } } #---------------------# # Ensenble Binary #### #-------------------# Convert.p<-th(projections.all.mean,s.m) projections.binary.all <- Convert.p$pa.raster writeRaster( projections.binary.all, filename = paste0("./outputs/", especie, "_","Ensemble Binary - Current Climate.tif"), format = "GTiff", overwrite = TRUE ) #--------------------# # Move the files #### #------------------# #install.packages("filesstrings") results<-list.files( "./outputs/",paste0(especie, "_"), full.names = TRUE ) file.move((list.files( "./outputs/",paste0(especie, "_"), full.names = TRUE )), (paste0("./outputs/", especie)), overwrite = TRUE) #--------------------# # Time Computing #### #------------------# sink("./outputs/tempo.txt", append = T) print(especie) print(Sys.time() - ini1) sink() } #END- 1
- 2
- 3
- 4
- 5
- 6
- 7
- 8
- 9
- 10
- 11
- 12
- 13
- 14
- 15
- 16
- 17
- 18
- 19
- 20
- 21
- 22
- 23
- 24
- 25
- 26
- 27
- 28
- 29
- 30
- 31
- 32
- 33
- 34
- 35
- 36
- 37
- 38
- 39
- 40
- 41
- 42
- 43
- 44
- 45
- 46
- 47
- 48
- 49
- 50
- 51
- 52
- 53
- 54
- 55
- 56
- 57
- 58
- 59
- 60
- 61
- 62
- 63
- 64
- 65
- 66
- 67
- 68
- 69
- 70
- 71
- 72
- 73
- 74
- 75
- 76
- 77
- 78
- 79
- 80
- 81
- 82
- 83
- 84
- 85
- 86
- 87
- 88
- 89
- 90
- 91
- 92
- 93
- 94
- 95
- 96
- 97
- 98
- 99
- 100
- 101
- 102
- 103
- 104
- 105
- 106
- 107
- 108
- 109
- 110
- 111
- 112
- 113
- 114
- 115
- 116
- 117
- 118
- 119
- 120
- 121
- 122
- 123
- 124
- 125
- 126
- 127
- 128
- 129
- 130
- 131
- 132
- 133
- 134
- 135
- 136
- 137
- 138
- 139
- 140
- 141
- 142
- 143
- 144
- 145
- 146
- 147
- 148
- 149
- 150
- 151
- 152
- 153
- 154
- 155
- 156
- 157
- 158
- 159
- 160
- 161
- 162
- 163
- 164
- 165
- 166
- 167
- 168
- 169
- 170
- 171
- 172
- 173
- 174
- 175
- 176
- 177
- 178
- 179
- 180
- 181
- 182
- 183
- 184
- 185
- 186
- 187
- 188
- 189
- 190
- 191
- 192
- 193
- 194
- 195
- 196
- 197
- 198
- 199
- 200
- 201
- 202
- 203
- 204
- 205
- 206
- 207
- 208
- 209
- 210
- 211
- 212
- 213
- 214
- 215
- 216
- 217
- 218
- 219
- 220
- 221
- 222
- 223
- 224
- 225
- 226
- 227
- 228
- 229
- 230
- 231
- 232
- 233
- 234
- 235
- 236
- 237
- 238
- 239
- 240
- 241
- 242
- 243
- 244
- 245
- 246
- 247
- 248
- 249
- 250
- 251
- 252
- 253
- 254
- 255
- 256
- 257
- 258
- 259
- 260
- 261
- 262
- 263
- 264
- 265
- 266
- 267
- 268
- 269
- 270
- 271
- 272
- 273
- 274
- 275
- 276
- 277
- 278
- 279
- 280
- 281
- 282
- 283
- 284
- 285
- 286
- 287
- 288
- 289
- 290
- 291
- 292
- 293
- 294
- 295
- 296
- 297
- 298
- 299
- 300
- 301
- 302
- 303
- 304
- 305
- 306
- 307
- 308
- 309
- 310
- 311
- 312
- 313
- 314
- 315
- 316
- 317
- 318
- 319
- 320
- 321
- 322
- 323
- 324
- 325
- 326
- 327
- 328
- 329
- 330
- 331
- 332
- 333
- 334
- 335
- 336
- 337
- 338
- 339
- 340
- 341
- 342
- 343
- 344
- 345
- 346
- 347
- 348
- 349
- 350
- 351
- 352
- 353
- 354
- 355
- 356
- 357
- 358
- 359
- 360
- 361
- 362
- 363
- 364
- 365
- 366
- 367
- 368
- 369
- 370
- 371
- 372
- 373
- 374
- 375
- 376
- 377
- 378
- 379
- 380
- 381
- 382
- 383
- 384
- 385
- 386
- 387
- 388
- 389
- 390
- 391
- 392
- 393
- 394
- 395
- 396
- 397
- 398
- 399
- 400
- 401
- 402
- 403
- 404
- 405
- 406
- 407
- 408
- 409
- 410
- 411
- 412
- 413
- 414
- 415
- 416
- 417
- 418
- 419
- 420
- 421
- 422
- 423
- 424
- 425
- 426
- 427
- 428
- 429
- 430
- 431
- 432
- 433
- 434
- 435
- 436
- 437
- 438
- 439
- 440
- 441
- 442
- 443
- 444
- 445
- 446
- 447
- 448
- 449
- 450
- 451
- 452
- 453
- 454
- 455
- 456
- 457
- 458
- 459
- 460
- 461
- 462
- 463
- 464
- 465
- 466
- 467
- 468
- 469
- 470
- 471
- 472
- 473
- 474
- 475
- 476
- 477
- 478
- 479
- 480
- 481
- 482
- 483
- 484
- 485
- 486
- 487
- 488
- 489
- 490
- 491
- 492
- 493
- 494
- 495
- 496
- 497
- 498
- 499
- 500
- 501
- 502
- 503
- 504
- 505
- 506
- 507
- 508
- 509
- 510
- 511
- 512
- 513
- 514
- 515
- 516
- 517
- 518
- 519
- 520
- 521
- 522
- 523
- 524
- 525
- 526
- 527
- 528
- 529
- 530
- 531
- 532
- 533
- 534
- 535
- 536
- 537
- 538
- 539
- 540
- 541
- 542
- 543
- 544
- 545
- 546
- 547
- 548
- 549
- 550
- 551
- 552
- 553
- 554
- 555
- 556
- 557
- 558
- 559
- 560
- 561
- 562
- 563
- 564
- 565
- 566
- 567
- 568
- 569
- 570
- 571
- 572
- 573
- 574
- 575
- 576
- 577
- 578
- 579
- 580
- 581
- 582
- 583
- 584
- 585
- 586
- 587
- 588
- 589
- 590
- 591
- 592
- 593
- 594
- 595
- 596
- 597
- 598
- 599
- 600
- 601
- 602
- 603
- 604
- 605
- 606
- 607
- 608
- 609
- 610
- 611
- 612
- 613
- 614
- 615
- 616
- 617
- 618
- 619
- 620
- 621
- 622
- 623
- 624
- 625
- 626
- 627
- 628
- 629
- 630
- 631
- 632
- 633
- 634
- 635
- 636
- 637
- 638
- 639
- 640
- 641
- 642
- 643
- 644
- 645
- 646
- 647
- 648
- 649
- 650
- 651
- 652
- 653
- 654
- 655
- 656
- 657
- 658
- 659
- 660
- 661
- 662
- 663
- 664
- 665
- 666
- 667
- 668
- 669
- 670
- 671
- 672
- 673
- 674
- 675
- 676
- 677
- 678
- 679
- 680
- 681
- 682
- 683
- 684
- 685
- 686
- 687
- 688
- 689
- 690
- 691
- 692
- 693
- 694
- 695
- 696
- 697
- 698
- 699
- 700
- 701
- 702
- 703
- 704
- 705
- 706
- 707
- 708
- 709
- 710
- 711
- 712
- 713
- 714
- 715
- 716
- 717
- 718
- 719
- 720
- 721
- 722
- 723
- 724
- 725
- 726
- 727
- 728
- 729
- 730
- 731
- 732
- 733
- 734
- 735
- 736
- 737
- 738
- 739
- 740
- 741
- 742
- 743
- 744
- 745
- 746
- 747
- 748
- 749
- 750
- 751
- 752
- 753
- 754
- 755
- 756
- 757
- 758
- 759
- 760
- 761
- 762
- 763
- 764
- 765
- 766
- 767
- 768
- 769
- 770
- 771
- 772
- 773
- 774
- 775
- 776
- 777
- 778
- 779
- 780
- 781
- 782
- 783
- 784
- 785
- 786
- 787
- 788
- 789
- 790
- 791
- 792
- 793
- 794
- 795
- 796
- 797
- 798
- 799
- 800
- 801
- 802
- 803
- 804
- 805
- 806
- 807
- 808
- 809
- 810
- 811
- 812
- 813
- 814
- 815
- 816
- 817
- 818
- 819
- 820
- 821
- 822
- 823
- 824
- 825
- 826
- 827
- 828
- 829
- 830
- 831
- 832
- 833
- 834
- 835
- 836
- 837
- 838
- 839
- 840
- 841
- 842
- 843
- 844
- 845
- 846
- 847
- 848
- 849
- 850
- 851
- 852
- 853
- 854
- 855
- 856
- 857
- 858
- 859
- 860
- 861
- 862
- 863
- 864
- 865
- 866
- 867
- 868
- 869
- 870
- 871
- 872
- 873
- 874
- 875
- 876
- 877
- 878
- 879
- 880
- 881
- 882
- 883
- 884
- 885
- 886
- 887
- 888
- 889
- 890
- 891
- 892
- 893
- 894
- 895
- 896
- 897
- 898
- 899
- 900
- 901
- 902
- 903
- 904
- 905
- 906
- 907
- 908
- 909
- 910
- 911
- 912
- 913
- 914
- 915
- 916
- 917
- 918
- 919
- 920
- 921
- 922
- 923
- 924
- 925
- 926
- 927
- 928
- 929
- 930
- 931
- 932
- 933
- 934
- 935
- 936
- 937
- 938
- 939
- 940
- 941
- 942
- 943
- 944
- 945
- 946
- 947
- 948
- 949
- 950
- 951
- 952
- 953
- 954
- 955
- 956
- 957
- 958
- 959
- 960
- 961
- 962
- 963
- 964
- 965
- 966
- 967
- 968
- 969
- 970
- 971
- 972
- 973
- 974
- 975
- 976
- 977
- 978
- 979
- 980
- 981
- 982
- 983
- 984
- 985
- 986
- 987
- 988
- 989
- 990
- 991
- 992
- 993
- 994
- 995
- 996
- 997
- 998
- 999
- 1000
- 1001
- 1002
- 1003
- 1004
- 1005
- 1006
- 1007
- 1008
- 1009
- 1010
- 1011
- 1012
- 1013
- 1014
- 1015
- 1016
- 1017
- 1018
- 1019
- 1020
- 1021
- 1022
- 1023
- 1024
- 1025
- 1026
- 1027
- 1028
- 1029
- 1030
- 1031
- 1032
- 1033
- 1034
- 1035
- 1036
- 1037
- 1038
- 1039
- 1040
- 1041
- 1042
- 1043
- 1044
- 1045
- 1046
- 1047
- 1048
- 1049
- 1050
- 1051
- 1052
- 1053
- 1054
- 1055
- 1056
- 1057
- 1058
- 1059
- 1060
- 1061
- 1062
- 1063
- 1064
- 1065
- 1066
- 1067
- 1068
- 1069
- 1070
- 1071
- 1072
- 1073
- 1074
- 1075
- 1076
- 1077
- 1078
- 1079
- 1080
- 1081
- 1082
- 1083
- 1084
- 1085
- 1086
- 1087
- 1088
- 1089
- 1090
- 1091
- 1092
- 1093
- 1094
- 1095
- 1096
- 1097
- 1098
- 1099
- 1100
- 1101
- 1102
- 1103
- 1104
- 1105
- 1106
- 1107
- 1108
- 1109
- 1110
- 1111
- 1112
- 1113
- 1114
- 1115
- 1116
- 1117
- 1118
- 1119
- 1120
- 1121
- 1122
- 1123
- 1124
- 1125
- 1126
- 1127
- 1128
- 1129
- 1130
- 1131
- 1132
- 1133
- 1134
- 1135
- 1136
- 1137
- 1138
- 1139
- 1140
- 1141
- 1142
- 1143
- 1144
- 1145
- 1146
- 1147
- 1148
- 1149
- 1150
- 1151
- 1152
- 1153
- 1154
- 1155
- 1156
- 1157
- 1158
- 1159
- 1160
- 1161
- 1162
- 1163
- 1164
- 1165
- 1166
- 1167
- 1168
- 1169
- 1170
- 1171
- 1172
- 1173
- 1174
- 1175
- 1176
- 1177
- 1178
- 1179
- 1180
- 1181
- 1182
- 1183
- 1184
- 1185
- 1186
- 1187
- 1188
- 1189
- 1190
- 1191
- 1192
- 1193
- 1194
- 1195
- 1196
- 1197
- 1198
- 1199
- 1200
- 1201
- 1202
- 1203
- 1204
- 1205
- 1206
- 1207
- 1208
- 1209
- 1210
- 1211
- 1212
- 1213
- 1214
- 1215
- 1216
- 1217
- 1218
- 1219
- 1220
- 1221
- 1222
- 1223
- 1224
- 1225
- 1226
- 1227
- 1228
- 1229
- 1230
- 1231
- 1232
- 1233
- 1234
- 1235
- 1236
- 1237
- 1238
- 1239
- 1240
- 1241
- 1242
- 1243
- 1244
- 1245
- 1246
- 1247
- 1248
- 1249
- 1250
- 1251
- 1252
- 1253
- 1254
- 1255
- 1256
- 1257
- 1258
- 1259
- 1260
- 1261
- 1262
- 1263
- 1264
- 1265
- 1266
- 1267
- 1268
- 1269
- 1270
- 1271
- 1272
- 1273
- 1274
- 1275
- 1276
- 1277
- 1278
- 1279
- 1280
- 1281
- 1282
- 1283
- 1284
- 1285
- 1286
- 1287
- 1288
- 1289
- 1290
- 1291
- 1292
- 1293
- 1294
- 1295
- 1296
- 1297
- 1298
- 1299
- 1300
- 1301
- 1302
- 1303
- 1304
- 1305
- 1306
- 1307
- 1308
- 1309
- 1310
- 1311
- 1312
- 1313
- 1314
- 1315
- 1316
- 1317
- 1318
- 1319
- 1320
- 1321
- 1322
- 1323
- 1324
- 1325
- 1326
- 1327
- 1328
- 1329
- 1330
- 1331
- 1332
- 1333
- 1334
- 1335
- 1336
- 1337
- 1338
- 1339
- 1340
- 1341
- 1342
- 1343
- 1344
- 1345
- 1346
- 1347
- 1348
- 1349
- 1350
- 1351
- 1352
- 1353
- 1354
- 1355
- 1356
- 1357
- 1358
- 1359
- 1360
- 1361
- 1362
- 1363
- 1364
- 1365
- 1366
- 1367
- 1368
- 1369
- 1370
- 1371
- 1372
- 1373
- 1374
- 1375
- 1376
- 1377
- 1378
- 1379
- 1380
- 1381
- 1382
- 1383
- 1384
- 1385
- 1386
- 1387
- 1388
- 1389
- 1390
- 1391
- 1392
- 1393
- 1394
- 1395
- 1396
- 1397
- 1398
- 1399
- 1400
- 1401
- 1402
- 1403
- 1404
- 1405
- 1406
- 1407
- 1408
- 1409
- 1410
- 1411
- 1412
- 1413
- 1414
- 1415
- 1416
- 1417
- 1418
- 1419
- 1420
- 1421
- 1422
- 1423
- 1424
- 1425
- 1426
- 1427
- 1428
- 1429
- 1430
- 1431
- 1432
- 1433
- 1434
- 1435
- 1436
- 1437
- 1438
- 1439
- 1440
- 1441
- 1442
- 1443
- 1444
- 1445
- 1446
- 1447
- 1448
- 1449
- 1450
- 1451
- 1452
- 1453
- 1454
- 1455
- 1456
- 1457
- 1458
- 1459
- 1460
- 1461
- 1462
- 1463
- 1464
- 1465
- 1466
- 1467
- 1468
- 1469
- 1470
- 1471
- 1472
- 1473
- 1474
- 1475
- 1476
- 1477
- 1478
- 1479
- 1480
- 1481
- 1482
- 1483
- 1484
- 1485
- 1486
- 1487
- 1488
- 1489
- 1490
- 1491
- 1492
- 1493
- 1494
- 1495
- 1496
- 1497
- 1498
- 1499
- 1500
- 1501
- 1502
- 1503
- 1504
- 1505
- 1506
- 1507
- 1508
- 1509
- 1510
- 1511
- 1512
- 1513
- 1514
- 1515
- 1516
- 1517
- 1518
- 1519
- 1520
- 1521
- 1522
- 1523
- 1524
- 1525
- 1526
- 1527
- 1528
- 1529
- 1530
- 1531
- 1532
- 1533
- 1534
- 1535
- 1536
- 1537
- 1538
- 1539
- 1540
- 1541
- 1542
- 1543
- 1544
- 1545
- 1546
- 1547
- 1548
- 1549
- 1550
- 1551
- 1552
- 1553
- 1554
- 1555
- 1556
- 1557
- 1558
- 1559
- 1560
- 1561
- 1562
- 1563
- 1564
- 1565
- 1566
- 1567
- 1568
- 1569
- 1570
- 1571
- 1572
- 1573
- 1574
- 1575
- 1576
- 1577
- 1578
- 1579
- 1580
- 1581
- 1582
- 1583
- 1584
- 1585
- 1586
- 1587
- 1588
- 1589
- 1590
- 1591
- 1592
- 1593
- 1594
- 1595
- 1596
- 1597
- 1598
- 1599
- 1600
-
相关阅读:
AI框架之Spring AI与Spring Cloud Alibaba AI使用讲解
BigDecimal保留两位小数失败问题
javascript 根据数组指定字段值，实现升序/降序
23面向对象案例1
太空网络攻击
基于ssm旅社客房收费管理系统+vue
通过电子文档管理改善患者体验和办公效率
word 导出 excel导出
Visual Studio：Entity设置表之间的关联关系
scrm系统哪款好？快鲸scrm私域流量运营专家
- 原文地址：https://blog.csdn.net/amyniez/article/details/134493988