Skip to content
Snippets Groups Projects
Commit 5988319b authored by Stefan Kranz's avatar Stefan Kranz
Browse files

check if pkgs are already installed

parent 3a4d2238
No related branches found
No related tags found
No related merge requests found
# 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)
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment