User Tools

Site Tools


wikistudholland:utwente_helen

Species Distribution Modeling with Swedish Wetland Vegetation

Helen Moor, Stockholm Resilience Centre, Stockholm University contact

INTRODUCTION

General framework of the analysis

Using species distribution modeling (MaxEnt, as appropriate for occurrence only data, Elith et al. 2010), we aim to predict the future distribution of the 100 functionally most important wetland species across Sweden. This projection will afterwards be used to assess potential future changes in regional species pool composition (focus on Norrstörm Drainage Basin, central Sweden), and concurrent effects on trait distributions for proposed effect traits. Ultimately, we will qualitatively discuss climate change impact on regional ecosystem service potential via species composition and trait distribution change (see e.g. Lavorel et al. 2011).

keywords: SDM, MaxEnt, functional traits, wetlands, vegetation, ecosystem services

Project objectives

Here the key objective is to process spatial data (vector and raster) from several sources to create a GRASS database. From that, ultimately create an input raster stack (with common projection, extent, resolution) to be used for species distribution modeling (using MaxEnt). Also, to practice the use and integration of bash, gdal/ogr, grass, and R.

If time, also 1) create pH raster from point data, predicting from soil/bedrock data (using GWR or Random Forests); 2) run MaxEnt on occurrence point data.

METADATA

Vector data

  • Sweden outline (e.g. available from GADM)
  • Norrström Drainage Basin (NDB) outline
  • Soil and Bedrock; CRS: Sweref99 TM (EPSG:3006); format: shp; available from Geological Survey of Sweden
  • pH point data; CRS: RT90 and Sweref99 TM; format: .csv; inbuilt coordinate error 1000m; potentially available upon request; original source: Swedish National Forest Inventory

Raster data

  • Worldclim current: Tmin, Tmax, Prec, Bioclimatic Variables; Resolution: 30 arc sec; CRS: WGS84 (EPSG:4326); format .bil; available: Climate current
  • Worldclim future (HadFGEM2-AO, 2070, rcp60): Tmin, Tmax, Prec, Bioclimatic Variables; Resolution: 30 arc sec; CRS: WGS84 (EPSG:4326); format .bil; available: Climate 2070
  • Heightmap (DEM) Sweden; Resolution: 50m, CRS: Sweref99 TM (EPSG:3006); format: .tif; downloaded via SU licence; original source: Lantmäteriet Heightmap
  • CORINE Landcover; Resolution: 100 m; CRS: EPSG:3035; format: .tif; available from EEA

Other data sources

  • Species point occurrence data, from Swedish National Wetland inventory (Våtmarksinventeringen); from local sqlite-database, but data now available at GBIF

METHOD

Bash script including GRASS commands. The script:

  • downloads data from publicly available sources (where applicable)
  • edits the downloaded data
  • imports the data into GRASS
  • further processes data in GRASS (e.g. calculates slope, aspect, twi from DEM)

Computational implementation

A Grass GISdatabase is first created manually in QGis. Once the database is established, the script downloads and edits data to import into the location europeWGS84, then creates a new location swedenSWEREF and proceeds to import data into this location.
Note: wget is not standard for MacOSX; see here for how to install it. Else: use curl.

Worldclim data - Grass database - Location Europe

- buildgrassdb.sh
#!/usr/bin/bash
######################################################################################### 
# Author Helen Moor
# Purpose: download and reformat worldclim climate data and other spatial data, 
# to compile a GRASS GIS database for Europe (Worldclim) and Sweden (all data)
# adapted to MacOSX 9, MacBook
#########################################################################################
 
 
#### 1. Create GRASS database
## prepare GRASS:
# MacOS 9: to be able to open grass in commandline, first:
# - add to .bash-profile: 
#	export PATH="/usr/X11R6/bin:$PATH"
#	alias grass="open -a GRASS-6.4.app"
# - open GRASS in existing GISDATABASE, and type the following to set -text to default: 
#	g.gui -un gui=text
 
 
##################################################################################
 
#### 2. Get data and import data into Grassdb:
 
######### Location EUROPE (for Worldclim data only)
# enter grass in GISDATABASE: GRASSDB, Location: europeWGS84, MAPSET:PERMANENT
grass ~/Documents/GISdata/GRASSDB/europeWGS84/PERMANENT
# check:
g.gisenv
g.region -p # extent is europe (set up in qgis)
# change resolution to 30 arcsecs:
g.region nsres=0:00:30 ewres=0:00:30  
 
#############################################
######### WORLDCLIM DATA 
# download, unzip and save raster data from worldclim 
# (.bil, 30 arc sec resolution, epsg 4326):
 
# cd to input working directory:
cd ~/Documents/GISdata/Worldclim
 
# create folders to hold downloaded data:
mkdir worldclim_current_tmp  worldclim_future_tmp
 
######### CURRENT CONDITIONS:
#####  1. download 5 zip files for current climate: tmin, tmax, prec, bio1-9, bio10-19
cd ~/Documents/GISdata/Worldclim/worldclim_current_tmp
 
## download, unzip, fix header, remove .zip:
 
for layer in tmin tmax prec bio1-9 bio10-19; do	
	curl -L -O "http://biogeo.ucdavis.edu/data/climate/worldclim/1_4/grid/cur/"$layer"_30s_bil.zip"  
	unzip $layer"_30s_bil" 
 
	# if else statement for different filenaming with respect to zipfilename
	if [ "$layer" == tmin ] || [ "$layer" == tmax ] || [ "$layer" == prec ] ; then 
 
	# fix the problem with the missing header: 
        # (http://duncanjg.wordpress.com/2010/08/09/import-worldclim-data-into-grass/):
		for i in $(seq 1 12); do 
			echo -e "PIXELTYPE     SIGNEDINT\r" >>$layer"_"$i.hdr; # append "PIXELTYPE SIGNEDINT" to each hdr file
		done 
 
	elif [ "$layer" == bio1-9 ]; then
		for i in $(seq 1 9); do 
			echo -e "PIXELTYPE     SIGNEDINT\r" >>"bio_"$i.hdr;
		done 
	else
		for i in $(seq 10 19); do 
			echo -e "PIXELTYPE     SIGNEDINT\r" >>"bio_"$i.hdr;
		done 
 
	fi	
	rm $layer"_30s_bil.zip"
done
 
# check what is there:
gdalinfo $layer"_1.bil"
 
##### 2. cut all .bil files to europe extent (same as in location europeWGS84); output as Gtiff 
 
for layer in tmin tmax prec bio; do	
	if [ "$layer" == bio ]; then
		for i in $(seq 1 19); do
			gdalwarp -te -25.1 34.9 35 71.3 $layer"_"$i.bil eu_$layer"_"$i.tif
		done
	else 
		for i in $(seq 1 12); do
			gdalwarp -te -25.1 34.9 35 71.3 $layer"_"$i.bil eu_$layer"_"$i.tif
		done
	fi
done
 
# check:
gdalinfo "eu_"$layer"_1.tif"  
 
## remove original bil/hdr files:
rm *.hdr; rm *.bil
 
##### 3. import the tifs into GRASS
 
# tmin, tmax, prec: 12 files for 12 months
for layer in tmin tmax prec bio; do	
	for i in $(seq 1 12); do r.in.gdal -o input=eu_$layer"_"$i.tif output=$layer$i; done
done
# biovars: 19 files
layer=bio
for i in $(seq 1 19); do r.in.gdal -o input=eu_$layer"_"$i.tif output=$layer$i; done
 
# check:
g.list type=rast
r.info map=tmin1
 
 
######### FUTURE CONDITIONS:
# Climate scenario: 2070 (70), GCM: HadGEM2-AO (HD), RCP: rcp60 (60)
# cd to folder for future:
cd ~/Documents/GISdata/Worldclim/worldclim_future_tmp/
mkdir -p europetmp 
 
##### 1. download and unzip (these files are in .tif already; type=float32)
 
for layer in tn tx pr bi; do	
	curl -L -O "http://biogeo.ucdavis.edu/data/climate/cmip5/30s/hd60"$layer"70.zip"
	unzip "hd60"$layer"70.zip"
	rm "hd60"$layer"70.zip"
	# clip to Europe extent
	if [ "$layer" == bi ] ; then  
		for i in $(seq 1 19); do
			gdalwarp -te -25.1 34.9 35 71.3 "hd60"$layer"70"$i".tif" europetmp/"eu_hd60"$layer"70"$i".tif"
		done
	else
		for i in $(seq 1 12); do
			gdalwarp -te -25.1 34.9 35 71.3 "hd60"$layer"70"$i".tif" europetmp/"eu_hd60"$layer"70"$i".tif"
		done
	fi
	rm "hd60"$layer"70"*.tif  # clean up on the go
done
 
## remove original bil/hdr files:
rm *.hdr; rm *.bil
 
##### 2. import the tifs into GRASS
 
# tmin, tmax, prec: 12 files for 12 months 
for layer in tn tx pr; do	
	for i in $(seq 1 12); do r.in.gdal -o input=europetmp/eu_hd60$layer"70"$i.tif output=70hd60$layer$i; done
done
# biovars: 19 files
layer=bi
for i in $(seq 1 19); do r.in.gdal -o input=europetmp/eu_hd60$layer"70"$i.tif output=70hd60$layer$i; done
 
# check:
g.list type=rast
r.info map=tmin1

Create new Location for Sweden

The script continues to create a new location for Sweden (extent follows a polygon outline of Sweden). Resolution is set to 50m, corresponding to highest resolution of all input rasters. Projection Sweref99TM (EPSG:3006). Note: grass does not know EPSG:3006, but recognises equivalence to UTM33N, EPSG:32633, and stores data using this CRS.

- buildgrassdb_continued.sh
 
#############################################
######### create new Location: SWEDEN
 
 
### 1. create new location:
# IMPORTANT: set resolution for region to 50 m from start 
# (else: on the fly nearest neighbour resampling)
# easiest with an input raster: 
# use gdal to clip and reproject one of my inputs to Sweden:
 
# get poly outline of Sweden:
#curl -L -O "http://biogeo.ucdavis.edu/data/gadm2/shp/SWE_adm.zip"
# or use own:
ogrinfo -so -al ~/Documents/GISdata/Sweden_outline/Sweden_WGS84.shp
 
## clip one of the Europe rasters to this outline:
cd ~/Documents/GISdata/Worldclim/worldclim_current_tmp
mkdir swedentmp
# clip European extent to outline of Sweden only (using polygon outline) 
# AND reproject to Sweref99TM (epsg3006) AND resample to 50m resolution:
gdalwarp -dstnodata -9999 -t_srs EPSG:3006 -tr 50 50 -cutline \
~/Documents/GISdata/Sweden_outline/Sweden_WGS84.shp -crop_to_cutline \
eu_tmin_1.tif swedentmp/swe_templ.tif
 
gdalinfo swedentmp/swe_templ.tif 
 
# create new location from this raster in existing GRASSDB: 
r.in.gdal input=~/Documents/GISdata/Worldclim/worldclim_current_tmp/swedentmp/swe_templ.tif output=swe_templ location=swedenSWEREF --overwrite 
 
# switch to new location:
g.mapset mapset=PERMANENT location=swedenSWEREF gisdbase=~/Documents/GISdata/GRASSDB
 
# remove the raster used for creating location (in grass and on HD):
g.remove rast=swe_templ 
rm swedentmp/swe_templ.tif
 
# switch to new location: not sure this works
# g.gisenv set="LOCATION_NAME=swedenSWEREF" ### seems to work, but prompt still says euopreWGS84
# exit grass: (safer)
exit

Worldclim Sweden

Continue to populate the new location with data for Sweden. First: Worldclim data:

- buildgrassdb_continued.sh
 
#############################################
### BUILD UP DATABASE FOR SWEDEN:
 
#  reopen GRASS in swedenSWEREF location, mapset PERMANENT:
grass ~/Documents/GISdata/GRASSDB/swedenSWEREF/PERMANENT
 
########## WORLDCLIM - import to grass
## CURRENT:
cd ~/Documents/GISdata/Worldclim/worldclim_current_tmp/
mkdir -p swedentmp  
 
# import all clipped and reprojected rasters into grass:
 
for layer in tmin tmax prec; do
	for i in $(seq 1 12); do
	# clip European extent to outline of Sweden only (using polygon outline) 
        # AND reproject to Sweref99TM (epsg3006):
	gdalwarp -dstnodata -9999 -t_srs EPSG:3006 \
          -cutline ~/Documents/GISdata/Sweden_outline/Sweden_WGS84.shp -crop_to_cutline \
          "eu_"$layer"_"$i".tif" swedentmp/swe_$layer$i.tif
	## import to GRASS
	r.in.gdal input=~/Documents/GISdata/Worldclim/worldclim_current_tmp/swedentmp/"swe_"$layer$i".tif" output="swe_"$layer$i --overwrite
	done
done
 
layer=bio
for i in $(seq 1 19); do
	# clip European extent to outline of Sweden only (using polygon outline) 
        # AND reproject to Sweref99TM (epsg3006):
	gdalwarp -dstnodata -9999 -t_srs EPSG:3006 -cutline ~/Documents/GISdata/Sweden_outline/Sweden_WGS84.shp \
          -crop_to_cutline "eu_"$layer"_"$i".tif" swedentmp/swe_$layer$i.tif
	## import to GRASS:
	r.in.gdal input=~/Documents/GISdata/Worldclim/worldclim_current_tmp/swedentmp/"swe_"$layer$i".tif" output="swe_"$layer$i --overwrite
done
 
# remove the tif files:
#rm -r swedentmp
#rm -r europetmp
 
## FUTURE:
cd ~/Documents/GISdata/Worldclim/worldclim_future_tmp/
mkdir -p swedentmp
 
for layer in tn tx pr; do
	for i in $(seq 1 12); do
	# clip European extent to outline of Sweden only (using polygon outline) 
        # AND reproject to Sweref99TM (epsg3006)
	gdalwarp -dstnodata -9999 -t_srs EPSG:3006 -cutline ~/Documents/GISdata/Sweden_outline/Sweden_WGS84.shp \
          -crop_to_cutline europetmp/"eu_hd60"$layer"70"$i".tif" swedentmp/"swe_70hd60"$layer$i.tif
	## import to GRASS
	r.in.gdal input=~/Documents/GISdata/Worldclim/worldclim_future_tmp/swedentmp/"swe_70hd60"$layer$i".tif" \                
          output="swe_70hd60"$layer$i --overwrite
	done
done
 
layer=bi
for i in $(seq 1 19); do
	# clip European extent to outline of Sweden only (using polygon outline) 
        # AND reproject to Sweref99TM (epsg3006):
	gdalwarp -dstnodata -9999 -t_srs EPSG:3006 -cutline ~/Documents/GISdata/Sweden_outline/Sweden_WGS84.shp \
          -crop_to_cutline europetmp/"eu_hd60"$layer"70"$i".tif" swedentmp/"swe_70hd60"$layer$i.tif
	## import to GRASS:
	r.in.gdal input=~/Documents/GISdata/Worldclim/worldclim_future_tmp/swedentmp/"swe_70hd60"$layer$i".tif" \ 
           output="swe_70hd60"$layer$i --overwrite
done
 
# remove the tif files:
#rm -r swedentmp
#rm -r europetmp

Land cover data

Continue to populate the new location with additional data for Sweden (continue working in GRASS). Second: CORINE land cover data:

- buildgrassdb_continued.sh
 
 
################################################
# OTHER DATASETS to add to Sweden:
################################################
######### CORINE landcover 100x100m from EEA, 
# epsg3035 Type=Byte (44 categories, with color table)
 
cd ~/Documents/GISdata/EEA_CORINE/
curl -L -O "http://www.eea.europa.eu/data-and-maps/data/corine-land-cover-2006-raster-2/clc-2006-100m/g100_06.zip" # needs -L to follow links in case file has been moved
unzip g100_06.zip
gdalinfo g100_06_version2011/g100_06.tif # Europe, proj ETRS_1989_LAEA (epsg3035), 100m, colour table 
 
## reproject, clip, and import for Sweden:
cd g100_06_version2011/
# clip European extent to outline of Sweden only (using polygon outline) 
# AND reproject to Sweref99TM (epsg3006):
gdalwarp -dstnodata -9999 -t_srs EPSG:3006 -cutline ~/Documents/GISdata/Sweden_outline/Sweden_WGS84.shp \
  -crop_to_cutline g100_06.tif swe_CORINE_100m.tif
## import to GRASS
r.in.gdal input=~/Documents/GISdata/EEA_CORINE/g100_06_version2011/swe_CORINE_100m.tif \
  output=CORINE_EEA --overwrite

Soil and bedrock types

Continue to populate the new location with additional data for Sweden (continue working in GRASS). Third: categorical soil type data (SGU):

- buildgrassdb_continued.sh
 
#################################################
######### SOIL DATA SGU  
# vector data, Sweref99TM
cd ~/Documents/GISdata/SGU2013/SGU2014/
 
# download and unzip:
for name in jord berggrund; do
curl -L -O http://www.sgu.se/dokument/service_sgu_publ/$name"1M.zip"
unzip $name"1M.zip" # unzips into folder
done
 
# rasterize : (are already in Sweref, 1:1M)
gdal_rasterize -a JG2 -a_nodata -9999 -tr 1000.0 1000.0 -l jona_jg2 \
  -ot Int32 jord/data/jona_jg2.shp jord/data/soil1km.tif
gdal_rasterize -a BRG -a_nodata -9999 -tr 1000.0 1000.0 -l berggrundsytor \
  -ot Int32 berggrund/berggrundsytor.shp berggrund/bedrock1km.tif
 
# clip to Sweden outline:
gdalwarp -dstnodata -9999 -t_srs EPSG:3006 -cutline ~/Documents/GISdata/Sweden_outline/Sweden_WGS84.shp \
  -crop_to_cutline jord/data/soil1km.tif jord/data/soil1km_cropped.tif
gdalwarp -dstnodata -9999 -t_srs EPSG:3006 -cutline ~/Documents/GISdata/Sweden_outline/Sweden_WGS84.shp \
  -crop_to_cutline berggrund/bedrock1km.tif berggrund/bedrock1km_cropped.tif
 
# clean up:
rm jord/data/soil1km.tif
rm berggrund/bedrock1km.tif
 
## create new mapset within swedenSWEREF: SGUsoil 
# (switches there automatically accord. to g.gisenv but not the prompt)
 
g.mapset -c mapset=SGUsoil location=swedenSWEREF gisdbase=~/Documents/GISdata/GRASSDB
 
# switch to new mapset: ("a fairly radical maneuvre")
#g.mapset mapset=SGUsoil
# or exit and re-enter, to be on safe side:
exit
grass ~/Documents/GISdata/GRASSDB/swedenSWEREF/SGUsoil
 
cd ~/Documents/GISdata/SGU2013/SGU2014/
 
# import rasterized and clipped tifs to grass:
r.in.gdal input=jord/data/soil1km_cropped.tif output=SGU_soil_1km
r.in.gdal input=berggrund/bedrock1km_cropped.tif output=SGU_bedrock_1km
 
exit

Vector data

Continue to populate the new location with additional data for Sweden (continue working in GRASS). Fourth: vector data (outlines of Sweden and the Norrström Drainage Basin):

- buildgrassdb_continued.sh
 
#########################################################################################
######### VECTOR DATA Sweden and Norrstrom drainage basin:
grass ~/Documents/GISdata/GRASSDB/swedenSWEREF/PERMANENT
 
##import vector sweden outline into grass
v.in.ogr -oe dsn=~/Documents/GISdata/Sweden_outline/Sweden_SWEREF99.shp output=vect_swe_boundary
g.list type=vect  # worked
v.info map=vect_swe_boundary
 
## and Norrstrom outline:
v.in.ogr -oe dsn=~/Documents/GISdata/NDB/NDB_border.shp output=vect_NDB_boundary

Digital Elevation Model

Continue to populate the new location with additional data for Sweden (continue working in GRASS). Fifth: Digital Elevation Model (DEM) of Sweden, resolution 50m, Sweref99TM. From the DEM, calculate additional rasters of aspect, slope, and topographic wetness index (twi) in GRASS:

- buildgrassdb_continued.sh
 
#########################################################################################
#########  DEM 50m res, from Lantmateriet:
#grass ~/Documents/GISdata/GRASSDB/swedenSWEREF/PERMANENT
 
cd ~/Documents/GISdata/SwedenDEM_Lantmateriet_Andrew/
gdalinfo SwedenDEM50mSWEREF.tif  # epsg 3006, 50mx50m, from tiles, 
# slightly larger than the vector outline (should be masked by region)
 
# import:
r.in.gdal input=~/Documents/GISdata/SwedenDEM_Lantmateriet_Andrew/SwedenDEM50mSWEREF.tif output=DEM_50M
 
## create new mapset for work with DEM: (switches automatically to new)
g.mapset -c mapset=DEMwork location=swedenSWEREF gisdbase=~/Documents/GISdata/GRASSDB
 
# copy DEM into working mapset DEMwork:  #not required, have access to PERMANENT
# g.copy rast='DEM_50M@PERMANENT',DEM_50M # NOT NECESSARY, can use any raster in location
 
## from DEM, calculate slope, aspect and topographic wetness index: 
 
# slope (degrees) and aspect (E N W S): 
r.slope.aspect -a elevation=DEM_50M slope=DEM50M_slope aspect=DEM50M_aspect
 
# topographic wetness index (ln(a/tan(beta)):
r.topidx input=DEM_50M output=DEM50M_twi
 
## :this took a long time, fills up memory and slows down exponentially;
## tried using GNU parallel (installed via homebrew), but aborted after warning
#parallel r.topidx input=DEM_50M output=DEM50M_twi
 
 
# export all 3 rasters as tif:
#cd ~/Documents/GISdata/SwedenDEM_Lantmateriet_Andrew/
mkdir -p derived/
 
# note: had to add createopt="TFW=YES" to export readable tifs 
 
for layer in aspect slope twi; do
r.out.gdal input=DEM50M_$layer format=GTiff type=Float32 output=DEM50m_$layer.tif nodata=-9999 createopt="TFW=YES" 
 
# the GTiff output from GRASS has dimensions 0x0, only readable because of TFW file; 
# therefore warp into proper tif, cut to polygon and give it some metatdata:
gdalwarp -dstnodata -9999 -s_srs EPSG:32633 -t_srs EPSG:3006 -ot Float32 -dstnodata -9999 -tr 50 50 \
  -cutline ~/Documents/GISdata/Sweden_outline/Sweden_WGS84.shp -crop_to_cutline DEM50m_$layer.tif \
   derived/DEM50m_$layer.tif
rm DEM50m_$layer.tif; rm DEM50m_$layer.tfw
done
 
#########################################################################################

To be able to use bedrock types as input to a random forest model, reclassify the raster. Originally contains 146 classes, but the randomForest function (in R package randomForest) cannot handle >32 levels of categorical data.

Reclassification was done manually via the following steps:

  • load original shape file into QGis
  • export attribute table as csv

Using the following R script:

  • import the csv
  • subset to unique categories, and export the csv
  • manually create aggregate classes and add to csv (OpenOffice)
  • re-import edited csv and create reclassification matrix from it
  • reclassify the raster
  • export as tif

Finally add to the GRASS database.

- reclassifyBedrock.R
# in QGis, exported attribute table of shape file to csv file.
 
# upload this csv :
rocks <- read.csv(file="~/Documents/GISdata/SGU2013/SGU2014/berggrund/berggrundsytor.csv")
 
# what's in it?:
names(rocks) # BRG is code, Lithology is label
 
# subset to code and label:
rocks <- rocks[,c(1,3)]
 
# subset to unique values:
rocks <- unique(rocks)
 
# sort alphabetically:
rocks <- with(rocks,  rocks[order(LITHOLOGY) , ])
 
# export table to csv::
write.csv(rocks, file="~/Documents/GISdata/SGU2013/SGU2014/berggrund/categories.csv", row.names=F)
 
# manually add aggregated classes - numbered and with label - in OpenOffice,
# then reload:
rocks <- read.csv(file="~/Documents/GISdata/SGU2013/SGU2014/berggrund/categories.csv")
 
# make rcl matrix: col1 "is", col2 "becomes":
reclass <- as.matrix(rocks[,c(1,2)])
 
# cd to dir with raster file, and reclassify original tif based on rcl:
setwd("~/Documents/GISdata/ph_predictors_sweref_50m")
reclassify(rock, reclass, filename='rock_reclassified')
 
system('ls') # written to .gri file - cannot use gdalwarp (unknown format)
 
# load that instead of old rock layer, to export then as tif:
rock <- raster("~/Documents/GISdata/ph_predictors_sweref_50m/rock_reclassified.gri")
 
# check values:
as.matrix(table(values(rock))) # classes 1:15 plus 999 for unknown
 
# make factor to work with it in R:
#rock <- as.factor(rock)
 
# export the reclassified raster to tif:
writeRaster(rock, "rock_reclassified.tif", format="GTiff")
 
# check if it worked:
system('gdalinfo rock_reclassified.tif')
 
# clean up:
rm(rocks, reclass)

Add the reclassified bedrock raster to the GRASS database:

- addBedrock_reclass.sh
 
cd ~/Documents/GISdata/ph_predictors_sweref_50m/
 
# warp to Sweref, data type: Int32
gdalwarp -dstnodata 0 -t_srs EPSG:3006 -ot Int32 -tr 50 50 rock_reclassified.tif rock_reclassified_sweref.tif
 
# import to GRASS, swedenSWEREF, SGUsoil:
grass ~/Documents/GISdata/GRASSDB/swedenSWEREF/SGUsoil
 
r.in.gdal input=~/Documents/GISdata/ph_predictors_sweref_50m/rock_reclassified_sweref.tif \ 
  output=SGU_bedrock_1km_reclassed
 
exit

Note: R for spatial data

4 different ways of loading raster data into R. From my experience, option Nr 4 is by far the best.

- loadRaster.R
##### 1. directly from GRASS:----------------------------------------------------
## in terminal, start grass:
# grass ~/Documents/GISdata/GRASSDB/swedenSWEREF/SGUsoil/
## open R from within a running GRASS session:
# open -a RStudio.app &
 
library(spgrass6)
 
# allows the use of grass commands from within R:
system('g.list type=rast')
 
## read grass raster: readRAST6()
 
soil1 <- readRAST6(vname="SGU_soil_1km")    
 
# takes longer than loading from a file 
# returns large SpatialGridDataFrame (eats up memory)
 
# vname can be vector of multiple files, 
# cat="SGU_soil_1km" should allow factor, but gives error for integer
 
system('exit')
 
 
##### 2. indirectly from GRASS:-------------------------------------------------------
 
library(raster)
 
### raster()- from cellhd in grass db: 
soil2 = raster("~/Documents/GISdata/GRASSDB/swedenSWEREF/SGUsoil/cellhd/SGU_soil_1km")
 
# error - doesnt recognize file format
 
 
##### 3. rgdal via tif:---------------------------------------------------------------
 
library(rgdal)
 
### readGDAL(): 
soil3 <- readGDAL(fname="~/Documents/GISdata/ph_predictors_sweref_50m/SGU_soil_50m.tif")
 
# returns large SpatialGridDataFrame (eats up memory)
# speed similar to readRAST6()
 
 
##### 4. raster() via tif: ---------------------------------------------------------------
 
library(raster)
 
### raster():
soil <- raster("~/Documents/GISdata/ph_predictors_sweref_50m/SGU_soil_50m.tif", values=T) 
 
# uses very little memory - creates link to raster, virtual representation

RESULTS

A neat GRASSDATABASE - stores a huge amount of raster data in just 5 GB, and is easily copy-pasted between machines.

http://spatial-ecology.net/dokuwiki/lib/exe/fetch.php?media=wikistudholland:grassdblayout.png

2 DO: create input raster stack to MaxEnt; run MaxEnt on all predictors

wikistudholland/utwente_helen.txt · Last modified: 2021/01/20 20:36 (external edit)