diff --git a/DEM2POINT.R b/DEM2POINT.R new file mode 100644 index 0000000000000000000000000000000000000000..6aaca5320f826f5dfd8544918f3b7bcc61881220 --- /dev/null +++ b/DEM2POINT.R @@ -0,0 +1,110 @@ +# 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) + +