diff --git a/DEM2POINT.R b/DEM2POINT.R index 6aaca5320f826f5dfd8544918f3b7bcc61881220..041a51426738d27f2e54529ab256fd455eb7f119 100644 --- a/DEM2POINT.R +++ b/DEM2POINT.R @@ -1,110 +1,117 @@ -# 1. PREPARATION #------------------------------------- -# DO ONLY ONCE - -# package to load many pckgs (to show how2install pckg) -install.packages("pacman") # install from web to PC -library(pacman) # load(attach) from PC's library folder - -pacman::p_load(raster, maptools, sp, dismo, rasterVis, rgdal) -# notice: when '#' the line is "neutral = Not run" - -# set path to DEM-raster -rasterpath<-"C:/Users/CHORO_04/Onedrive/PRACTICALS/AREALKUNDE/SDM" -# set path to species dots shapefile - -# 2. GET THE DATA #----------------------------------------- -# DEM raster -setwd(rasterpath) -elev <- raster("CHELSA-DEM.tif") #from the defined path - -# Create Study area "SWA" -ext <- extent(elev) # get extent of DEM raster -SWA <- as(ext, 'SpatialPolygons')#convert ext to a SpatialPolygon -proj4string(SWA) <- CRS("+init=epsg:4326") #define KBS WGS84=4326 - -# Species occurrences -setwd("C:/Users/CHORO_04/OneDrive/AREALKUNDE/PRACTICALS/SPEC01") -# the resp. path to species dots -# you only need to change *THIS PATH* and *ELEVATION* range for -# all next species - -dots <- maptools::readShapePoints("DOTS.shp") -proj4string(dots) <- CRS("+init=epsg:4326") #define KBS WGS84=4326 - -# 3. OUT GBIF DOTS #---------------------------------- -SWA_DOTS <- dots[SWA,] # select 341 which 'in SWA' -# = like the 'select by location' step in QGIS - -# jetzt die selektieren, die "CTS-SWA..." sind - -# Vorher: Replacing NA's with "None" for all SWA_DOTS -# muss aber leider vorher als level hinzugefuegt werden -# get levels and add "None" -levels <- levels(SWA_DOTS@data$datasetNam) -levels[length(levels) + 1] <- "None" -# refactor datasetNam to include "None" as level -# and replace NA with "None" -SWA_DOTS@data$datasetNam <- factor(SWA_DOTS@data$datasetNam, - levels = levels) -SWA_DOTS@data$datasetNam[is.na(SWA_DOTS@data$datasetNam)] <- "None" - -# jetzt endlich unsere 91 CTS auswaehlen -CTS_DOTS<-SWA_DOTS[(SWA_DOTS@data$datasetNam == "CTS-SWA-9-1"),] - -# 4. EXTRACT DEM #------------------------------------------ - -# define circles with a radius of 15 km around SWA_DOTS -# with "dismo::" I make sure the "circles" function works -# by calling it directly from dismo-package -BUFF <- dismo::circles(CTS_DOTS,d=15000,lonlat=T) - -# everybody: draw 10000 random points within BUFF circles -SAMPLE <- spsample(BUFF@polygons, 10000, type='random') -# das dauert am laengsten/takes the longest time - -# extract DEM values at point positions -ElevDots<-extract(elev,SAMPLE,df=T) # "point sampling tool" -hist(ElevDots$CHELSA.DEM) # to see the data - -# 5. SUBSET BY DEM #-------------------------- - -# In CTS-SWA: "Pergularia tomentosa is a light requiring species -# occurring in open places, insolated, and on dry sandy or -# gravelly-sandy-clayey soils, frequently together with Artemisia -# and various semishrublets (suffrutex). In vertical distribution -# it has been found in depressions around the Dead Sea, more or -# less at -250 m and usually it does not grow beyond 900-1200 m -# elevation. The most elevated localities have been found in -# Pakistan at 1330 m and in Afghanistan at 1 500-1600 m." -# -> In this case I choose -250 to 1250m and neglect the outliers. - -# Join the (spatial)SAMPLE and ElevDots data (both n = 10000) -SAMPLEDOTS <- SpatialPointsDataFrame(SAMPLE, ElevDots) -summary(SAMPLEDOTS) -# Subsetting an SPDF in R is easy. -# You can do it just like a normal dataframe. -# remember "CHELSA.DEM" is the elevation -SMPLDTS<-SAMPLEDOTS[(SAMPLEDOTS@data$CHELSA.DEM > -250),] -SMPLDT <-SMPLDTS[(SMPLDTS@data$CHELSA.DEM < 1250),] -hist(SMPLDT@data$CHELSA.DEM, - col="red", - add=T) # to control - -# 6. PLOT THE DATA #----------------------------------------- -plot(SWA) -plot(elev, add=T) -plot(SWA, add=T) -plot(dots,col="black",pch=20,cex=0.5,add=T) -plot(CTS_DOTS, col="red",pch=16,cex=0.5,add=T) -plot(BUFF, cex=5, add=T) #enlarged by 5 -plot(SMPLDT, col="blue", pch=16,cex=0.1,add=T) - -# 7. PLOT 3D FOR FUN #-------------------------------------- -require(rasterVis) -plot3D(elev, - maxpixels=1e6, - zfac=0.1, - drape=NULL, - col=terrain.colors) - - +## Install of additional packages +## pkgs are will only be downloaded, if they aren't installed yet + +# pacman is currently only available for R 3.5 and higher, which isn't available for Linux LTS distros +#if(!require(pacman)) { +# install.packages("pacman") +# library(pacman) +#} + +pkgNeeded <- c("base", "raster", "maptools", "sp", "dismo", "rasterVis", "rgdal") +isInstalled <- sapply(pkgNeeded, require, character.only = TRUE) +print(isInstalled) + + +pacman::p_load("raster", "maptools", "sp", "dismo", "rasterVis", "rgdal") +# notice: when '#' the line is "neutral = Not run" + +# set path to DEM-raster +rasterpath<-"C:/Users/CHORO_04/Onedrive/PRACTICALS/AREALKUNDE/SDM" +# set path to species dots shapefile + +# 2. GET THE DATA #----------------------------------------- +# DEM raster +setwd(rasterpath) +elev <- raster("CHELSA-DEM.tif") #from the defined path + +# Create Study area "SWA" +ext <- extent(elev) # get extent of DEM raster +SWA <- as(ext, 'SpatialPolygons')#convert ext to a SpatialPolygon +proj4string(SWA) <- CRS("+init=epsg:4326") #define KBS WGS84=4326 + +# Species occurrences +setwd("C:/Users/CHORO_04/OneDrive/AREALKUNDE/PRACTICALS/SPEC01") +# the resp. path to species dots +# you only need to change *THIS PATH* and *ELEVATION* range for +# all next species + +dots <- maptools::readShapePoints("DOTS.shp") +proj4string(dots) <- CRS("+init=epsg:4326") #define KBS WGS84=4326 + +# 3. OUT GBIF DOTS #---------------------------------- +SWA_DOTS <- dots[SWA,] # select 341 which 'in SWA' +# = like the 'select by location' step in QGIS + +# jetzt die selektieren, die "CTS-SWA..." sind + +# Vorher: Replacing NA's with "None" for all SWA_DOTS +# muss aber leider vorher als level hinzugefuegt werden +# get levels and add "None" +levels <- levels(SWA_DOTS@data$datasetNam) +levels[length(levels) + 1] <- "None" +# refactor datasetNam to include "None" as level +# and replace NA with "None" +SWA_DOTS@data$datasetNam <- factor(SWA_DOTS@data$datasetNam, + levels = levels) +SWA_DOTS@data$datasetNam[is.na(SWA_DOTS@data$datasetNam)] <- "None" + +# jetzt endlich unsere 91 CTS auswaehlen +CTS_DOTS<-SWA_DOTS[(SWA_DOTS@data$datasetNam == "CTS-SWA-9-1"),] + +# 4. EXTRACT DEM #------------------------------------------ + +# define circles with a radius of 15 km around SWA_DOTS +# with "dismo::" I make sure the "circles" function works +# by calling it directly from dismo-package +BUFF <- dismo::circles(CTS_DOTS,d=15000,lonlat=T) + +# everybody: draw 10000 random points within BUFF circles +SAMPLE <- spsample(BUFF@polygons, 10000, type='random') +# das dauert am laengsten/takes the longest time + +# extract DEM values at point positions +ElevDots<-extract(elev,SAMPLE,df=T) # "point sampling tool" +hist(ElevDots$CHELSA.DEM) # to see the data + +# 5. SUBSET BY DEM #-------------------------- + +# In CTS-SWA: "Pergularia tomentosa is a light requiring species +# occurring in open places, insolated, and on dry sandy or +# gravelly-sandy-clayey soils, frequently together with Artemisia +# and various semishrublets (suffrutex). In vertical distribution +# it has been found in depressions around the Dead Sea, more or +# less at -250 m and usually it does not grow beyond 900-1200 m +# elevation. The most elevated localities have been found in +# Pakistan at 1330 m and in Afghanistan at 1 500-1600 m." +# -> In this case I choose -250 to 1250m and neglect the outliers. + +# Join the (spatial)SAMPLE and ElevDots data (both n = 10000) +SAMPLEDOTS <- SpatialPointsDataFrame(SAMPLE, ElevDots) +summary(SAMPLEDOTS) +# Subsetting an SPDF in R is easy. +# You can do it just like a normal dataframe. +# remember "CHELSA.DEM" is the elevation +SMPLDTS<-SAMPLEDOTS[(SAMPLEDOTS@data$CHELSA.DEM > -250),] +SMPLDT <-SMPLDTS[(SMPLDTS@data$CHELSA.DEM < 1250),] +hist(SMPLDT@data$CHELSA.DEM, + col="red", + add=T) # to control + +# 6. PLOT THE DATA #----------------------------------------- +plot(SWA) +plot(elev, add=T) +plot(SWA, add=T) +plot(dots,col="black",pch=20,cex=0.5,add=T) +plot(CTS_DOTS, col="red",pch=16,cex=0.5,add=T) +plot(BUFF, cex=5, add=T) #enlarged by 5 +plot(SMPLDT, col="blue", pch=16,cex=0.1,add=T) + +# 7. PLOT 3D FOR FUN #-------------------------------------- +require(rasterVis) +plot3D(elev, + maxpixels=1e6, + zfac=0.1, + drape=NULL, + col=terrain.colors) + +