R script by Welk

DEM2POINT.R

+#  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)
+
+