require(raster) require(rgdal) require(plyr) rm(list=ls()) #_____________________________________________________________________________________________________________________________ # Environmental heterogeneity quantification on oceanic islands #ŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻ # Load environmental heterogeneity components setwd(" ") # set the directory where you stored the data islands_poly <- readOGR("islands.shp") # load island shapefile (polygons) temp <-raster("CHELSA_temp_1979-2013_V1_1MeanannualtemperatureVersion1.1.tif") # download http://chelsa-climate.org/downloads/ prec <-raster("CHELSA_prec_1979-2013_V1_1MeanannualprecipitationVersion1.1.tif") # download http://chelsa-climate.org/downloads/ elev <-raster("srtm_250.tif") # download http://srtm.csi.cgiar.org/srtmdata/ heat <-raster("heat.tif") # This was calculated from the elevation data (loaded in line 15) using the Geomorphometry & Gradient Metrics toolbox version 2.0-0 (Evans et al. 2014) in ArcGIS version 10.4. Check https://github.com/jeffreyevans/GradientMetrics # Note: The elevation raster layer we provide in Drayad was preprocessed after download: # 1. Tiles were merged 2. Raster layer was masked for islands only 3. Negative values were eliminated. # Quantification of whole-island metrics v <- extract(prec, islands_poly) # Island range min<- unlist(lapply(v, function(x) if (!is.null(x)) min(x, na.rm=TRUE) else NA )) max<- unlist(lapply(v, function(x) if (!is.null(x)) max(x, na.rm=TRUE) else NA )) PRECrg <- max-min v <- extract(temp, islands_poly) min<- unlist(lapply(v, function(x) if (!is.null(x)) min(x, na.rm=TRUE) else NA )) max<- unlist(lapply(v, function(x) if (!is.null(x)) max(x, na.rm=TRUE) else NA )) TEMPrg <- max-min v <- extract(elev, islands_poly) min<- unlist(lapply(v, function(x) if (!is.null(x)) min(x, na.rm=TRUE) else NA )) max<- unlist(lapply(v, function(x) if (!is.null(x)) max(x, na.rm=TRUE) else NA )) ELEVrg <- max-min v <- extract(heat, islands_poly) min<- unlist(lapply(v, function(x) if (!is.null(x)) min(x, na.rm=TRUE) else NA )) max<- unlist(lapply(v, function(x) if (!is.null(x)) max(x, na.rm=TRUE) else NA )) HLIrg <- max-min v <- extract(prec, islands_poly) # Island standard deviation PRECsd <- unlist(lapply(v, function(x) if (!is.null(x)) sd(x, na.rm=TRUE) else NA )) v <- extract(temp, islands_poly) TEMPsd <- unlist(lapply(v, function(x) if (!is.null(x)) sd(x, na.rm=TRUE) else NA )) v <- extract(elev, islands_poly) ELEVsd <- unlist(lapply(v, function(x) if (!is.null(x)) sd(x, na.rm=TRUE) else NA )) v <- extract(heat, islands_poly) HLIsd <- unlist(lapply(v, function(x) if (!is.null(x)) sd(x, na.rm=TRUE) else NA )) # Quantification moving-window metrics: summarized EH per island by mean values of the heterogeneity rasters # Heterogeneity rasters notes: # All heterogeneity rasters were computed using ArcGIS with the Toolbox for Surface Gradient and Geomorphometric Modeling, version 2.0-0 (Evans et al. 2014) # You can calculate the heterogeneity rasters in R: # Use the `dissection` function from package spatialEco v1.2-1 https://www.rdocumentation.org/packages/spatialEco/versions/1.2-1/topics/dissection # Use the `tri` function (i.e., Terrain Ruggedness Index) from package spatialEco v1.2-1 https://www.rdocumentation.org/packages/spatialEco/versions/1.2-1/topics/tri # Use `focal` function from the raster package and calculate stadard deviation in a selected window size https://www.rdocumentation.org/packages/raster/versions/3.0-12/topics/focal # A window size 9 km² is comprised by 12*12 pixels when using 250 m elevation data and 3*3 pixels when using 1 km² CHELSA data # Load rasters layers provided in Dryad elev_msd<- raster("elev_msd12.tif") elev_dis<- raster("elev_dis12.tif") elev_rou<- raster("elev_rou12.tif") prec_msd <- raster("prec_msd3.tif") prec_dis <- raster("prec_dis3.tif") prec_rou <- raster("prec_rou3.tif") temp_msd <- raster("temp_msd3.tif") temp_dis <- raster("temp_dis3.tif") temp_rou <- raster("temp_rou3.tif") heat_msd <- raster("heat_msd12.tif") heat_dis <- raster("heat_dis12.tif") heat_rou <- raster("heat_rou12.tif") # Once you have the heterogeneity rasters, (you computed them either using a GIS software or R) you can summarized values per island: v <- extract(elev_dis, islands_poly) ELEVdis <- unlist(lapply(v, function(x) if (!is.null(x)) mean(x, na.rm=TRUE) else NA )) v <- extract(elev_msd, islands_poly) ELEVmsd <- unlist(lapply (v, function(x) if (!is.null(x)) mean(x, na.rm=TRUE) else NA )) v <- extract(elev_rou, islands_poly) ELEVrou <- unlist(lapply(v, function(x) if (!is.null(x)) mean(x, na.rm=TRUE) else NA )) v <- extract(prec_dis, islands_poly) PRECdis <- unlist(lapply(v, function(x) if (!is.null(x)) mean(x, na.rm=TRUE) else NA )) v <- extract(prec_msd, islands_poly) PRECmsd <- unlist(lapply(v, function(x) if (!is.null(x)) mean(x, na.rm=TRUE) else NA )) v <- extract(prec_rou, islands_poly) PRECrou <- unlist(lapply(v, function(x) if (!is.null(x)) mean(x, na.rm=TRUE) else NA )) v <- extract(temp_dis, islands_poly) TEMPdis <- unlist(lapply(v, function(x) if (!is.null(x)) mean(x, na.rm=TRUE) else NA )) v <- extract(temp_msd, islands_poly) TEMPmsd <- unlist(lapply(v, function(x) if (!is.null(x)) mean(x, na.rm=TRUE) else NA )) v <- extract(temp_rou, islands_poly) TEMProu <- unlist(lapply(v, function(x) if (!is.null(x)) mean(x, na.rm=TRUE) else NA )) v <- extract(heat_dis, islands_poly) HLIdis <- unlist(lapply(v, function(x) if (!is.null(x)) mean(x, na.rm=TRUE) else NA )) v <- extract(heat_msd, islands_poly) HLImsd <- unlist(lapply(v, function(x) if (!is.null(x)) mean(x, na.rm=TRUE) else NA )) v <- extract(heat_rou, islands_poly) HLIrou <- unlist(lapply(v, function(x) if (!is.null(x)) mean(x, na.rm=TRUE) else NA )) eh_metrics <- data.frame( PRECrg, TEMPrg, ELEVrg, HLIrg, PRECsd, TEMPsd, ELEVsd, HLIsd, PRECdis, PRECmsd, PRECrou, TEMPdis, TEMPmsd, TEMProu, ELEVdis, ELEVmsd, ELEVrou, HLIdis, HLImsd, HLIrou) eh_metrics <- cbind(islands_poly@data[1], eh_metrics)# Adding island names to EH metrics eh_metrics <- data.frame( entt_ID = eh_metrics$entt_ID, ELEVmsd = eh_metrics$ELEVmsd, ELEVsd = eh_metrics$ELEVsd, ELEVrou = eh_metrics$ELEVrou, ELEVdis = eh_metrics$ELEVdis, ELEVrg = eh_metrics$ELEVrg, PRECmsd = eh_metrics$PRECmsd, PRECsd = eh_metrics$PRECsd, PRECrou = eh_metrics$PRECrou, PRECdis = eh_metrics$PRECdis, PRECrg = eh_metrics$PRECrg, TEMPmsd = eh_metrics$TEMPmsd+10, # +10 to avoid cero when log-tranforming. TEMPsd = eh_metrics$TEMPsd+10, TEMProu = eh_metrics$TEMProu+10, TEMPdis = eh_metrics$TEMPdis+10, TEMPrg = eh_metrics$TEMPrg+10, HLImsd = eh_metrics$HLImsd, HLIsd = eh_metrics$HLIsd, HLIrou = eh_metrics$HLIrou, HLIdis = eh_metrics$HLIdis, HLIrg = eh_metrics$HLIrg) #_____________________________________________________________________________________________________________________________ # Test of the effect of Window size and spatial grain on the EH quantification per island #ŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻ # Note 1: Raster layers of elevation at 500 m and 1 km² were obtained from elevation at 250 m. We used in ArcGIS 10.4 the billinear resampling technique from Resample function. # You can also resample a raster in R using `resample` function from raster package in R. # Note 2: We used the 500 m and 1 km² layers and calculated standard deviation in three additional window sizes 3, 25 and 49 km². # This produced heterogeneity rasters with the different spatial grains calculated with three additional window sizes: # Load rasters layers provided in dryad in zipfile: Scale_effect_analysis.zip elev250_msd_3 <- raster("srtm250m_msd7.tif") # Window size 3 km² comprised 7*7 pixels elev250_msd_25 <- raster("srtm250m_msd20.tif") # Window size 25 km² comprised 20*20 pixels elev250_msd_49 <- raster("srtm250m_msd28.tif") # Window size 49 km² comprised 28*28 pixels elev500m_msd_49<- raster("srtm15asec_msd7.tif") # Window size 49 km² comprised 14*14 pixels elev1km_msd_9 <- raster("srtm30asec_msd3.tif") # Window size 9 km² comprised 3*3 pixels elev1km_msd_49 <- raster("srtm30asec_msd7.tif") # Window size 49 km² comprised 7*7 pixels # Summarized EH per island by calculating mean values of the heterogeneity rasters v <- extract(elev250_msd_3, islands_poly) elev250msd_3 <- unlist(lapply(v, function(x) if (!is.null(x)) mean(x, na.rm=TRUE) else NA )) v <- extract(elev250_msd_9, islands_poly) elev250msd_9 <- unlist(lapply(v, function(x) if (!is.null(x)) mean(x, na.rm=TRUE) else NA )) v <- extract(elev250_msd_25, islands_poly) elev250msd_25 <- unlist(lapply(v, function(x) if (!is.null(x)) mean(x, na.rm=TRUE) else NA )) v <- extract(elev250_msd_49, islands_poly) elev250msd_49 <- unlist(lapply(v, function(x) if (!is.null(x)) mean(x, na.rm=TRUE) else NA )) v <- extract(elev500m_msd_49, islands_poly) elev500msd_49 <- unlist(lapply(v, function(x) if (!is.null(x)) mean(x, na.rm=TRUE) else NA )) v <- extract(elev1km_msd_9, islands_poly) elev1msd_9 <- unlist(lapply(v, function(x) if (!is.null(x)) mean(x, na.rm=TRUE) else NA )) v <- extract(elev1km_msd_49, islands_poly) elev1msd_49 <- unlist(lapply(v, function(x) if (!is.null(x)) mean(x, na.rm=TRUE) else NA )) test <- data.frame(elev250msd_3, elev250msd_9, elev250msd_25, elev250msd_49, elev500msd_49, elev1msd_9, elev1msd_49) #_______________________ # Figure 6 Appendix S2 #ŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻ require(GGally) ggpairs(data=test, columns=1:7, mapping= NULL, diag = list(continuous = wrap("densityDiag")), lower = list(continuous = wrap("points", size=0.8)))