6. Standard tools for spatial data analysis

Foreword on working directory, data and packages

The Data for this tutorial are provided via Github. In order to reproudce the code with your own data, replace the url with your local filepath

If you work locally, you may define the working directory

setwd("YOUR/FILEPATH") 
# this will not work. if you want to use it, define a path on your computer

This should be the place where your data and other folders are stored. more Information can be found here

Code chunks for saving data and results are disabeled in this tutorial (using a #). Please replace (“YOUR/FILEPATH”) with your own path.

The code is written in such a way, that packages will be automatically installed (if not already installed) and loaded. Packages will be stored in the default location of your computer.

#Intro Simple features, i.e. records that come with a geometry, can be manipulated or changed. for the case where these manipulations involve geometries, different functions than with data frames are necessary. This is the content of this tutorial.

We will perform a number of preparational steps in order to create a map of tropical nights in the near future. The area of interest is the Trinational Metropolitan Area Upper Rhine Valley (TMO). The data on tropical nights were modeled by the German Weather Service (DWD) for the Clim’Ability Project.

We will work with the sf package.

if (!require(sf)){install.packages("sf"); library(sf)}
## Lade nötiges Paket: sf
## Warning: Paket 'sf' wurde unter R Version 3.6.2 erstellt
## Linking to GEOS 3.6.1, GDAL 2.2.3, PROJ 4.9.3

Methods for simple feature objects

A complete list of methods for a particular class can be obtained with methods(class = "sf"). To see al objects in the package use ls("package:sf") and for all functions us lsf("package:sf"). Additional functions can be useful from the sp package

methods(class = "sf")
##  [1] $<-                   [                     [[<-                 
##  [4] aggregate             as.data.frame         cbind                
##  [7] coerce                dbDataType            dbWriteTable         
## [10] identify              initialize            merge                
## [13] plot                  print                 rbind                
## [16] show                  slotsFromS3           st_agr               
## [19] st_agr<-              st_area               st_as_sf             
## [22] st_bbox               st_boundary           st_buffer            
## [25] st_cast               st_centroid           st_collection_extract
## [28] st_convex_hull        st_coordinates        st_crop              
## [31] st_crs                st_crs<-              st_difference        
## [34] st_geometry           st_geometry<-         st_interpolate_aw    
## [37] st_intersection       st_intersects         st_is                
## [40] st_join               st_line_merge         st_nearest_points    
## [43] st_node               st_normalize          st_point_on_surface  
## [46] st_polygonize         st_precision          st_segmentize        
## [49] st_set_precision      st_simplify           st_snap              
## [52] st_sym_difference     st_transform          st_triangulate       
## [55] st_union              st_voronoi            st_wrap_dateline     
## [58] st_write              st_zm                
## see '?methods' for accessing help and source code

Data set

We will use data from the Clim’Ability project. More specifically, modeled Climate dat for the near future in a bounding box over the TMO.

# load  shp fromfrom Github
# load from Github
download.file("https://raw.githubusercontent.com/RafHo/teaching/master/angewandte_geodatenverarbeitung/datasource/eca_tr_tasmin_yearDiff_ensmean_EUR-11_rcp45_day_2021-2050minus1971-2000_bas_rhin_neu.zip", destfile = "eca_tr_tasmin_yearDiff_ensmean_EUR-11_rcp45_day_2021-2050minus1971-2000_bas_rhin_neu.zip")

# unzip shapefile
unzip("eca_tr_tasmin_yearDiff_ensmean_EUR-11_rcp45_day_2021-2050minus1971-2000_bas_rhin_neu.zip")

tropical_nights_point_rcp45_2021_2050 <- st_read("eca_tr_tasmin_yearDiff_ensmean_EUR-11_rcp45_day_2021-2050minus1971-2000_bas_rhin_neu.shp")
## Reading layer `eca_tr_tasmin_yearDiff_ensmean_EUR-11_rcp45_day_2021-2050minus1971-2000_bas_rhin_neu' from data source `C:\Users\Nils Riach\Github\spatial_thinking_website\blogdown\content\courses\angewandte_geodatenverarbeitung\eca_tr_tasmin_yearDiff_ensmean_EUR-11_rcp45_day_2021-2050minus1971-2000_bas_rhin_neu.shp' using driver `ESRI Shapefile'
## Simple feature collection with 450 features and 1 field
## geometry type:  POINT
## dimension:      XY
## bbox:           xmin: 279173.6 ymin: 5184926 xmax: 521812.3 ymax: 5501534
## epsg (SRID):    NA
## proj4string:    +proj=utm +zone=32 +ellps=GRS80 +units=m +no_defs

The name of the layer gives us a couple of info on what the data set contains:

  • eca_tr_tasmin_yearDiff
    • tropical nights per year
  • ensmean
    • mean of the model ensemble
  • EUR-11
    • resolution of the data (EUR-11, ~12.5km).
  • rcp45
  • day
    • unit of change signal
  • 2021-2050minus1971-2000
    • time frame of change signal

Explore the data set

plot(tropical_nights_point_rcp45_2021_2050)

summary(tropical_nights_point_rcp45_2021_2050)
##     tr_____                 geometry  
##  Min.   :0.03056   POINT        :450  
##  1st Qu.:1.17569   epsg:NA      :  0  
##  Median :1.71528   +proj=utm ...:  0  
##  Mean   :1.81673                      
##  3rd Qu.:2.30208                      
##  Max.   :4.26389

Prepare the data set

To keep control of the changes we will perform on the data set, we will copy the feature every time we will perform an operation on it. be careful, though, since it is easy to mix up data sets. Usage of understandable name conventions is highly recommended. Especially if you want to reuse the script in the future!

names(tropical_nights_point_rcp45_2021_2050) <- c("trop_nights", "geometry")
tr1 <- tropical_nights_point_rcp45_2021_2050
summary(tropical_nights_point_rcp45_2021_2050)
##   trop_nights               geometry  
##  Min.   :0.03056   POINT        :450  
##  1st Qu.:1.17569   epsg:NA      :  0  
##  Median :1.71528   +proj=utm ...:  0  
##  Mean   :1.81673                      
##  3rd Qu.:2.30208                      
##  Max.   :4.26389

Recap: using the right crs

EPSG code is unknown but proj4string is known.However, the proj4string turns out to be identical with that of EPSG 25832 - ETRS89 / UTM zone 32N. We will set it to that.

st_crs(tr1)
## Coordinate Reference System:
##   No EPSG code
##   proj4string: "+proj=utm +zone=32 +ellps=GRS80 +units=m +no_defs"
tr1_25832 <- st_set_crs(tr1, 25832)

st_crs(tr1_25832)
## Coordinate Reference System:
##   EPSG: 25832 
##   proj4string: "+proj=utm +zone=32 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs"

Points to polygons

Temperture data is normally represented as continous area. In our case, the data is only available as point features. We will create so called voronoi polygons from the points. Since the st_voronoi function is not performing as expected, we will use the voronoi function from the disom package. For that, we must convert our soimple feature to class sp.

# convert to class spatial
tr1_25832_spatial <- as(tr1_25832, 'Spatial')


#load package dismo
if (!require(dismo)){install.packages("dismo"); library(dismo)}
## Lade nötiges Paket: dismo
## Warning: Paket 'dismo' wurde unter R Version 3.6.3 erstellt
## Lade nötiges Paket: raster
## Warning: Paket 'raster' wurde unter R Version 3.6.2 erstellt
## Lade nötiges Paket: sp
## Warning: Paket 'sp' wurde unter R Version 3.6.2 erstellt
if (!require(deldir)){install.packages("deldir"); library(deldir)}
## Lade nötiges Paket: deldir
## Warning: Paket 'deldir' wurde unter R Version 3.6.2 erstellt
## deldir 0.1-25
tr1_25832_voronoi <- dismo::voronoi(tr1_25832_spatial)

plot(tr1_25832_voronoi)
plot(tr1_25832_spatial, add=TRUE)

Clip to area of interest

We will load the area of interest, which is a dissolved shapefile of all the communities within the TMO. We now want to clip our voronoiy polygons the our area of interest.

# add the administrative boundaries
download.file("https://raw.githubusercontent.com/RafHo/teaching/master/angewandte_geodatenverarbeitung/datasource/Administrative_Boundaries_dissolved.zip", destfile = "Administrative_Boundaries_dissolved.zip")

# unzip shapefile
unzip("Administrative_Boundaries_dissolved.zip")

ab <- st_read("Administrative_Boundaries_dissolved.shp")
## Reading layer `Administrative_Boundaries_dissolved' from data source `C:\Users\Nils Riach\Github\spatial_thinking_website\blogdown\content\courses\angewandte_geodatenverarbeitung\Administrative_Boundaries_dissolved.shp' using driver `ESRI Shapefile'
## Simple feature collection with 1 feature and 20 fields
## geometry type:  POLYGON
## dimension:      XY
## bbox:           xmin: 336371.5 ymin: 5214597 xmax: 492299.8 ymax: 5466220
## epsg (SRID):    25832
## proj4string:    +proj=utm +zone=32 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs
str(ab)
## Classes 'sf' and 'data.frame':   1 obs. of  21 variables:
##  $ id        : Factor w/ 1 level "1": 1
##  $ name      : Factor w/ 1 level "Landau in der Pfalz": 1
##  $ type      : Factor w/ 1 level "village": 1
##  $ nuts0     : Factor w/ 1 level "DEU": 1
##  $ n0_name   : Factor w/ 1 level "Deutschland": 1
##  $ nuts1     : Factor w/ 1 level "DEB": 1
##  $ n1_name   : Factor w/ 1 level "Rheinland-Pfalz": 1
##  $ nuts2     : Factor w/ 1 level "DEB3": 1
##  $ n2_name   : Factor w/ 1 level "Rheinhessen-Pfalz": 1
##  $ nuts3     : Factor w/ 1 level "DEB33": 1
##  $ n3_name   : Factor w/ 1 level "Landau in der Pfalz": 1
##  $ id_osm    : Factor w/ 1 level "62391": 1
##  $ id_ags    : Factor w/ 1 level "07313000": 1
##  $ id_insee  : Factor w/ 0 levels: NA
##  $ id_ofs    : Factor w/ 0 levels: NA
##  $ id_order08: Factor w/ 1 level "07313000": 1
##  $ area_km2  : num 77.3
##  $ perim_km  : num 85.9
##  $ pop06     : num 42028
##  $ pop13     : num 43641
##  $ geometry  :sfc_POLYGON of length 1; first list element: List of 2
##   ..$ : num [1:25207, 1:2] 364707 364707 364767 364816 364872 ...
##   ..$ : num [1:92, 1:2] 391059 391261 391311 391375 391464 ...
##   ..- attr(*, "class")= chr  "XY" "POLYGON" "sfg"
##  - attr(*, "sf_column")= chr "geometry"
##  - attr(*, "agr")= Factor w/ 3 levels "constant","aggregate",..: NA NA NA NA NA NA NA NA NA NA ...
##   ..- attr(*, "names")= chr  "id" "name" "type" "nuts0" ...
# plot area of interest
plot(tr1_25832_voronoi)
plot(ab, add=TRUE)
## Warning in plot.sf(ab, add = TRUE): ignoring all but the first attribute

# remove unnecessary columns. Only keep geometry column
ab <- ab[ ,"geometry"]
str(ab)
## Classes 'sf' and 'data.frame':   1 obs. of  1 variable:
##  $ geometry:sfc_POLYGON of length 1; first list element: List of 2
##   ..$ : num [1:25207, 1:2] 364707 364707 364767 364816 364872 ...
##   ..$ : num [1:92, 1:2] 391059 391261 391311 391375 391464 ...
##   ..- attr(*, "class")= chr  "XY" "POLYGON" "sfg"
##  - attr(*, "sf_column")= chr "geometry"
##  - attr(*, "agr")= Factor w/ 3 levels "constant","aggregate",..:
#ab$id <- NULL
str(ab)
## Classes 'sf' and 'data.frame':   1 obs. of  1 variable:
##  $ geometry:sfc_POLYGON of length 1; first list element: List of 2
##   ..$ : num [1:25207, 1:2] 364707 364707 364767 364816 364872 ...
##   ..$ : num [1:92, 1:2] 391059 391261 391311 391375 391464 ...
##   ..- attr(*, "class")= chr  "XY" "POLYGON" "sfg"
##  - attr(*, "sf_column")= chr "geometry"
##  - attr(*, "agr")= Factor w/ 3 levels "constant","aggregate",..:
# now back to an sf object
voronoi_sf <- st_as_sf(tr1_25832_voronoi)

Spatial clipping is a form of spatial subsetting that involves changes to the geometry columns of at least some of the affected features. Clipping is only possible on mor complex features than points: lines polygons and their respective multi-equivalents.

Clipping can be performed by the st_intersection function:

voronoi_clip <- st_intersection( voronoi_sf, ab)
## Warning: attribute variables are assumed to be spatially constant throughout all
## geometries
plot(voronoi_clip)

Comparing changes

Now we will look at the values in the object:

plot(voronoi_clip$trop_nights)

hist(voronoi_clip$trop_nights)

boxplot(voronoi_clip$trop_nights)

# compare the cliped and thje original values
par(mfrow=c(1,2))
boxplot(voronoi_clip$trop_nights, main= "cliped area")
boxplot(tr1$trop_nights, main="original files")

Create map(s)

Lastly we prepare a map with the tmap package

if (!require(tmap)){install.packages("tmap"); library(tmap)}
## Lade nötiges Paket: tmap
## Warning: Paket 'tmap' wurde unter R Version 3.6.3 erstellt
#quick and dirty
qtm(voronoi_clip)

qtm(voronoi_clip, fill = "trop_nights", fill.pallete = "RdYlGn")

# more sophisticated
tm_shape(voronoi_clip) + tm_fill(col = "trop_nights", palette = "Reds", style = "jenks", 
                             title = "Tropical Nights", 
                             legend.hist = TRUE) + 
  #tm_borders f?gt die Grenzen hier hinzu. alpha kl?rt Transparenz 
  tm_borders(alpha= .2, col = "black") + 
  tm_compass() +
  tm_scale_bar() +
  tm_credits("Data: DWD [2019]; 2021-2050 RCP4.5", size = 0.7, position=c("left","bottom"))+
  tm_layout(title = NA,  legend.outside= FALSE, 
            legend.position = c("right", "center"),
            attr.outside= FALSE,  attr.position= c("RIGHT", "BOTTOM"), 
            legend.frame = FALSE, legend.width = -0.35, legend.hist.width = 0.3, legend.hist.height = 0.24,
            inner.margins = c(0.07, 0.03, 0.01, 0.27),
            frame= TRUE,
            design.mode= FALSE)

Now with different classes. This is done by changeign the style argument. Let’s see what changes.

#1st
equal <-tm_shape(voronoi_clip) + tm_fill(col = "trop_nights", palette = "Reds", style = "equal", 
                             title = "Tropical Nights (equal)", 
                             legend.hist = TRUE) + 
  #tm_borders f?gt die Grenzen hier hinzu. alpha kl?rt Transparenz 
  tm_borders(alpha= .2, col = "black") + 
  tm_compass() +
  tm_scale_bar() +
  tm_credits("Data: DWD [2019]", size = 0.7, position=c("left","bottom"))+
  tm_layout(title = NA,  legend.outside= FALSE, 
            legend.position = c("right", "center"),
            attr.outside= FALSE,  attr.position= c("RIGHT", "BOTTOM"), 
            legend.frame = FALSE, legend.width = -0.35, legend.hist.width = 0.3, legend.hist.height = 0.24,
            inner.margins = c(0.07, 0.03, 0.01, 0.27),
            frame= TRUE,
            design.mode= FALSE)
#2nd
pretty <- tm_shape(voronoi_clip) + tm_fill(col = "trop_nights", palette = "Reds", style = "pretty", 
                             title = "Tropical Nights (pretty)", 
                             legend.hist = TRUE) + 
  #tm_borders f?gt die Grenzen hier hinzu. alpha kl?rt Transparenz 
  tm_borders(alpha= .2, col = "black") + 
  tm_compass() +
  tm_scale_bar() +
  tm_credits("Data: DWD [2019]", size = 0.7, position=c("left","bottom"))+
  tm_layout(title = NA,  legend.outside= FALSE, 
            legend.position = c("right", "center"),
            attr.outside= FALSE,  attr.position= c("RIGHT", "BOTTOM"), 
            legend.frame = FALSE, legend.width = -0.35, legend.hist.width = 0.3, legend.hist.height = 0.24,
            inner.margins = c(0.07, 0.03, 0.01, 0.27),
            frame= TRUE,
            design.mode= FALSE)
#3rd
kmeans <- tm_shape(voronoi_clip) + tm_fill(col = "trop_nights", palette = "Reds", style = "kmeans", 
                             title = "Tropical Nights (kmeans)", 
                             legend.hist = TRUE) + 
  #tm_borders f?gt die Grenzen hier hinzu. alpha kl?rt Transparenz 
  tm_borders(alpha= .2, col = "black") + 
  tm_compass() +
  tm_scale_bar() +
  tm_credits("Data: DWD [2019]", size = 0.7, position=c("left","bottom"))+
  tm_layout(title = NA,  legend.outside= FALSE, 
            legend.position = c("right", "center"),
            attr.outside= FALSE,  attr.position= c("RIGHT", "BOTTOM"), 
            legend.frame = FALSE, legend.width = -0.35, legend.hist.width = 0.3, legend.hist.height = 0.24,
            inner.margins = c(0.07, 0.03, 0.01, 0.27),
            frame= TRUE,
            design.mode= FALSE)
#4th
 sd <- tm_shape(voronoi_clip) + tm_fill(col = "trop_nights", palette = "Reds", style = "sd", 
                             title = "Tropical Nights (sd)", 
                             legend.hist = TRUE) + 
  #tm_borders f?gt die Grenzen hier hinzu. alpha kl?rt Transparenz 
  tm_borders(alpha= .2, col = "black") + 
  tm_compass() +
  tm_scale_bar() +
  tm_credits("Data: DWD [2019]", size = 0.7, position=c("left","bottom"))+
  tm_layout(title = NA,  legend.outside= FALSE, 
            legend.position = c("right", "center"),
            attr.outside= FALSE,  attr.position= c("RIGHT", "BOTTOM"), 
            legend.frame = FALSE, legend.width = -0.35, legend.hist.width = 0.3, legend.hist.height = 0.24,
            inner.margins = c(0.07, 0.03, 0.01, 0.27),
            frame= TRUE,
            design.mode= FALSE)

current.mode <- tmap_mode("plot") 
## tmap mode set to plotting
tmap_arrange(equal, pretty, kmeans, sd)

tmap_mode(current.mode)
## tmap mode set to plotting

——-End——-

Previous
Next