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 4.1.3 erstellt
## Linking to GEOS 3.10.2, GDAL 3.4.1, PROJ 7.2.1; sf_use_s2() is TRUE

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] $<-                          [                           
##  [3] [[<-                         aggregate                   
##  [5] as.data.frame                cbind                       
##  [7] coerce                       dbDataType                  
##  [9] dbWriteTable                 filter                      
## [11] identify                     initialize                  
## [13] merge                        plot                        
## [15] print                        rbind                       
## [17] show                         slotsFromS3                 
## [19] st_agr                       st_agr<-                    
## [21] st_area                      st_as_s2                    
## [23] st_as_sf                     st_as_sfc                   
## [25] st_bbox                      st_boundary                 
## [27] st_buffer                    st_cast                     
## [29] st_centroid                  st_collection_extract       
## [31] st_convex_hull               st_coordinates              
## [33] st_crop                      st_crs                      
## [35] st_crs<-                     st_difference               
## [37] st_drop_geometry             st_filter                   
## [39] st_geometry                  st_geometry<-               
## [41] st_inscribed_circle          st_interpolate_aw           
## [43] st_intersection              st_intersects               
## [45] st_is                        st_is_valid                 
## [47] st_join                      st_line_merge               
## [49] st_m_range                   st_make_valid               
## [51] st_minimum_rotated_rectangle st_nearest_points           
## [53] st_node                      st_normalize                
## [55] st_point_on_surface          st_polygonize               
## [57] st_precision                 st_reverse                  
## [59] st_sample                    st_segmentize               
## [61] st_set_precision             st_shift_longitude          
## [63] st_simplify                  st_snap                     
## [65] st_sym_difference            st_transform                
## [67] st_triangulate               st_union                    
## [69] st_voronoi                   st_wrap_dateline            
## [71] st_write                     st_z_range                  
## [73] st_zm                        transform                   
## 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 `V:\github\nils_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
## Bounding box:  xmin: 279173.6 ymin: 5184926 xmax: 521812.3 ymax: 5501534
## Projected CRS: UTM_Zone_32_Northern_Hemisphere

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:
##   User input: UTM_Zone_32_Northern_Hemisphere 
##   wkt:
## PROJCRS["UTM_Zone_32_Northern_Hemisphere",
##     BASEGEOGCRS["GCS_GRS 1980(IUGG, 1980)",
##         DATUM["D_unknown",
##             ELLIPSOID["GRS80",6378137,298.257222101,
##                 LENGTHUNIT["metre",1,
##                     ID["EPSG",9001]]]],
##         PRIMEM["Greenwich",0,
##             ANGLEUNIT["Degree",0.0174532925199433]]],
##     CONVERSION["UTM zone 32N",
##         METHOD["Transverse Mercator",
##             ID["EPSG",9807]],
##         PARAMETER["Latitude of natural origin",0,
##             ANGLEUNIT["Degree",0.0174532925199433],
##             ID["EPSG",8801]],
##         PARAMETER["Longitude of natural origin",9,
##             ANGLEUNIT["Degree",0.0174532925199433],
##             ID["EPSG",8802]],
##         PARAMETER["Scale factor at natural origin",0.9996,
##             SCALEUNIT["unity",1],
##             ID["EPSG",8805]],
##         PARAMETER["False easting",500000,
##             LENGTHUNIT["metre",1],
##             ID["EPSG",8806]],
##         PARAMETER["False northing",0,
##             LENGTHUNIT["metre",1],
##             ID["EPSG",8807]],
##         ID["EPSG",16032]],
##     CS[Cartesian,2],
##         AXIS["(E)",east,
##             ORDER[1],
##             LENGTHUNIT["metre",1,
##                 ID["EPSG",9001]]],
##         AXIS["(N)",north,
##             ORDER[2],
##             LENGTHUNIT["metre",1,
##                 ID["EPSG",9001]]]]
tr1_25832 <- st_set_crs(tr1, 25832)
## Warning: st_crs<- : replacing crs does not reproject data; use st_transform for
## that
st_crs(tr1_25832)
## Coordinate Reference System:
##   User input: EPSG:25832 
##   wkt:
## PROJCRS["ETRS89 / UTM zone 32N",
##     BASEGEOGCRS["ETRS89",
##         DATUM["European Terrestrial Reference System 1989",
##             ELLIPSOID["GRS 1980",6378137,298.257222101,
##                 LENGTHUNIT["metre",1]]],
##         PRIMEM["Greenwich",0,
##             ANGLEUNIT["degree",0.0174532925199433]],
##         ID["EPSG",4258]],
##     CONVERSION["UTM zone 32N",
##         METHOD["Transverse Mercator",
##             ID["EPSG",9807]],
##         PARAMETER["Latitude of natural origin",0,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8801]],
##         PARAMETER["Longitude of natural origin",9,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8802]],
##         PARAMETER["Scale factor at natural origin",0.9996,
##             SCALEUNIT["unity",1],
##             ID["EPSG",8805]],
##         PARAMETER["False easting",500000,
##             LENGTHUNIT["metre",1],
##             ID["EPSG",8806]],
##         PARAMETER["False northing",0,
##             LENGTHUNIT["metre",1],
##             ID["EPSG",8807]]],
##     CS[Cartesian,2],
##         AXIS["(E)",east,
##             ORDER[1],
##             LENGTHUNIT["metre",1]],
##         AXIS["(N)",north,
##             ORDER[2],
##             LENGTHUNIT["metre",1]],
##     USAGE[
##         SCOPE["Engineering survey, topographic mapping."],
##         AREA["Europe between 6°E and 12°E: Austria; Belgium; Denmark - onshore and offshore; Germany - onshore and offshore; Norway including - onshore and offshore; Spain - offshore."],
##         BBOX[38.76,6,83.92,12]],
##     ID["EPSG",25832]]

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 4.1.3 erstellt
## Lade nötiges Paket: raster
## Warning: Paket 'raster' wurde unter R Version 4.1.3 erstellt
## Lade nötiges Paket: sp
## Warning: Paket 'sp' wurde unter R Version 4.1.3 erstellt
if (!require(deldir)){install.packages("deldir"); library(deldir)}
## Lade nötiges Paket: deldir
## deldir 1.0-6      Nickname: "Mendacious Cosmonaut"
## 
##      The syntax of deldir() has had an important change. 
##      The arguments have been re-ordered (the first three 
##      are now "x, y, z") and some arguments have been 
##      eliminated.  The handling of the z ("tags") 
##      argument has been improved.
##  
##      The "dummy points" facility has been removed. 
##      This facility was a historical artefact, was really 
##      of no use to anyone, and had hung around much too 
##      long.  Since there are no longer any "dummy points", 
##      the structure of the value returned by deldir() has 
##      changed slightly.  The arguments of plot.deldir() 
##      have been adjusted accordingly; e.g. the character 
##      string "wpoints" ("which points") has been 
##      replaced by the logical scalar "showpoints". 
##      The user should consult the help files.
tr1_25832_voronoi <- dismo::voronoi(tr1_25832_spatial)
## Warning in proj4string(x): CRS object has comment, which is lost in output; in tests, see
## https://cran.r-project.org/web/packages/sp/vignettes/CRS_warnings.html
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj
## = prefer_proj): Discarded datum Unknown based on GRS80 ellipsoid in Proj4
## definition
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 
##   `V:\github\nils_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
## Bounding box:  xmin: 336371.5 ymin: 5214597 xmax: 492299.8 ymax: 5466220
## Projected CRS: ETRS89 / UTM zone 32N
str(ab)
## Classes 'sf' and 'data.frame':   1 obs. of  21 variables:
##  $ id        : chr "1"
##  $ name      : chr "Landau in der Pfalz"
##  $ type      : chr "village"
##  $ nuts0     : chr "DEU"
##  $ n0_name   : chr "Deutschland"
##  $ nuts1     : chr "DEB"
##  $ n1_name   : chr "Rheinland-Pfalz"
##  $ nuts2     : chr "DEB3"
##  $ n2_name   : chr "Rheinhessen-Pfalz"
##  $ nuts3     : chr "DEB33"
##  $ n3_name   : chr "Landau in der Pfalz"
##  $ id_osm    : chr "62391"
##  $ id_ags    : chr "07313000"
##  $ id_insee  : chr NA
##  $ id_ofs    : chr NA
##  $ id_order08: chr "07313000"
##  $ 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 [1:3] "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 [1:20] "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 [1:3] "XY" "POLYGON" "sfg"
##  - attr(*, "sf_column")= chr "geometry"
##  - attr(*, "agr")= Factor w/ 3 levels "constant","aggregate",..: 
##   ..- attr(*, "names")= chr(0)
#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 [1:3] "XY" "POLYGON" "sfg"
##  - attr(*, "sf_column")= chr "geometry"
##  - attr(*, "agr")= Factor w/ 3 levels "constant","aggregate",..: 
##   ..- attr(*, "names")= chr(0)
# 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. First of all, the sf_objects make sure that the epsg code is the same or defined:

st_crs(voronoi_sf) <- 25832 
## Warning: st_crs<- : replacing crs does not reproject data; use st_transform for
## that
st_crs(ab) <- 25832

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
#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)
## Warning: argument design.mode not used anymore. Please use the stand-alone
## function tmap_design_mode() to set the design mode globally.

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)
## Warning: argument design.mode not used anymore. Please use the stand-alone
## function tmap_design_mode() to set the design mode globally.
#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)
## Warning: argument design.mode not used anymore. Please use the stand-alone
## function tmap_design_mode() to set the design mode globally.
#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)
## Warning: argument design.mode not used anymore. Please use the stand-alone
## function tmap_design_mode() to set the design mode globally.
#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)
## Warning: argument design.mode not used anymore. Please use the stand-alone
## function tmap_design_mode() to set the design mode globally.
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