First steps with SF
Jesper
20/6/2021
Learning the SF package
SF Vignette 1
In this note I take some first steps in learning to work with spatial objects in R using the sf-package. There are off course different spatial objects: Points, Lines and Polygons, to mention just the basics. In this note I use polygons as example. They are often important for map-making exercises which is also going to be one of the examples given here:
- Load a shape file with multipolygons
- Merge data to the sf object
- Plot it.
Secondly, I want to take a look at some of the basic methods available for manipulation of sf-objects in particular the polygons (so points and lines must wait for another note).
SF objects are data.frames with list column
The objects are data.frames with at least one column being a list column. Compared to a matrix, the benefit of a data.frame, is that it can store different types of data in each column. However, compared to lists, data.frames are still limited because they have a single number of rows, hence each column is equally long.
So how can you initialize a data.frame where one column is a list column? First how not to do it
a <- 1:3
b <- letters[1:3]
c <- list(a=c(1,2,3),b=c(2,3,4),c=c(3,4,5))
data.frame(a,b,c)## a b a.1 b.1 c
## 1 1 a 1 2 3
## 2 2 b 2 3 4
## 3 3 c 3 4 5secondly how to it using the operator
data.frame(a,b,I(c))## a b c
## a 1 a 1, 2, 3
## b 2 b 2, 3, 4
## c 3 c 3, 4, 5as explained here https://stackoverflow.com/questions/9547518/create-a-data-frame-where-a-column-is-a-list/9547594.
Reading a shapefile
path_shape <- "C://Users//np83zg//OneDrive - Aalborg Universitet//Skrivebord//cars-master//rspace_osm//Learning SF//shapes"
setwd(path_shape)
library(sf)## Warning: package 'sf' was built under R version 4.0.5## Linking to GEOS 3.9.0, GDAL 3.2.1, PROJ 7.2.1 mun_shp <- st_read("municipalities_DK.shp",options = "ENCODING=UTF-8")## options: ENCODING=UTF-8
## Reading layer `municipalities_DK' from data source `C:\Users\np83zg\OneDrive - Aalborg Universitet\Skrivebord\cars-master\rspace_osm\Learning SF\shapes\municipalities_DK.shp' using driver `ESRI Shapefile'
## Simple feature collection with 98 features and 2 fields
## Geometry type: MULTIPOLYGON
## Dimension: XYZ
## Bounding box: xmin: 441676 ymin: 6049782 xmax: 892787 ymax: 6402217
## z_range: zmin: -999 zmax: 0
## Projected CRS: ETRS89 / UTM zone 32N # There are problems with danish names ... don't know how to solve
# Using rgdal instead
# If reading a shapefile, the data source name (dsn= argument) is the
# folder (directory) where the shapefile is, and the layer is the name
# of the shapefile (without the .shp extension).
library(rgdal)## Warning: package 'rgdal' was built under R version 4.0.5## Loading required package: sp## Warning: package 'sp' was built under R version 4.0.5## rgdal: version: 1.5-23, (SVN revision 1121)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 3.2.1, released 2020/12/29
## Path to GDAL shared files: C:/Users/np83zg/OneDrive - Aalborg Universitet/Dokumenter/R/win-library/4.0/rgdal/gdal
## GDAL binary built with GEOS: TRUE
## Loaded PROJ runtime: Rel. 7.2.1, January 1st, 2021, [PJ_VERSION: 721]
## Path to PROJ shared files: C:/Users/np83zg/OneDrive - Aalborg Universitet/Dokumenter/R/win-library/4.0/rgdal/proj
## PROJ CDN enabled: FALSE
## Linking to sp version:1.4-5
## To mute warnings of possible GDAL/OSR exportToProj4() degradation,
## use options("rgdal_show_exportToProj4_warnings"="none") before loading rgdal.
## Overwritten PROJ_LIB was C:/Users/np83zg/OneDrive - Aalborg Universitet/Dokumenter/R/win-library/4.0/rgdal/proj foo <- readOGR(dsn=".", layer="municipalities_DK")## Warning in OGRSpatialRef(dsn, layer, morphFromESRI = morphFromESRI, dumpSRS =
## dumpSRS, : Discarded datum European_Terrestrial_Reference_System_1989 in Proj4
## definition: +proj=utm +zone=32 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m
## +no_defs## OGR data source with driver: ESRI Shapefile
## Source: "C:\Users\np83zg\OneDrive - Aalborg Universitet\Skrivebord\cars-master\rspace_osm\Learning SF\shapes", layer: "municipalities_DK"
## with 98 features
## It has 2 fields## Warning in readOGR(dsn = ".", layer = "municipalities_DK"): Z-dimension
## discarded # This returns S4 onject. You can get slotnames using
slotNames(foo)## [1] "data" "polygons" "plotOrder" "bbox" "proj4string" # The access slots using
foo@data## muni_name ID_muni
## 0 København 101
## 1 Frederiksberg 147
## 2 Ballerup 151
## 3 Brøndby 153
## 4 Dragør 155
## 5 Gentofte 157
## 6 Gladsaxe 159
## 7 Glostrup 161
## 8 Herlev 163
## 9 Albertslund 165
## 10 Hvidovre 167
## 11 Høje Taastrup 169
## 12 Lyngby-Taarbæk 173
## 13 Rødovre 175
## 14 Ishøj 183
## 15 TÃ¥rnby 185
## 16 Vallensbæk 187
## 17 Furesø 190
## 18 Allerød 201
## 19 Fredensborg 210
## 20 Helsingør 217
## 21 Hillerød 219
## 22 Hørsholm 223
## 23 Rudersdal 230
## 24 Egedal 240
## 25 Frederikssund 250
## 26 Greve 253
## 27 Køge 259
## 28 Halsnæs 260
## 29 Roskilde 265
## 30 Solrød 269
## 31 Gribskov 270
## 32 Odsherred 306
## 33 Holbæk 316
## 34 Faxe 320
## 35 Kalundborg 326
## 36 Ringsted 329
## 37 Slagelse 330
## 38 Stevns 336
## 39 Sorø 340
## 40 Lejre 350
## 41 Lolland 360
## 42 Næstved 370
## 43 Guldborgsund 376
## 44 Vordingborg 390
## 45 Bornholm 400
## 46 Middelfart 410
## 47 Assens 420
## 48 Faaborg-Midtfyn 430
## 49 Kerteminde 440
## 50 Nyborg 450
## 51 Odense 461
## 52 Svendborg 479
## 53 Nordfyns 480
## 54 Langeland 482
## 55 Ærø 492
## 56 Haderslev 510
## 57 Billund 530
## 58 Sønderborg 540
## 59 Tønder 550
## 60 Esbjerg 561
## 61 Fanø 563
## 62 Varde 573
## 63 Vejen 575
## 64 Aabenraa 580
## 65 Fredericia 607
## 66 Horsens 615
## 67 Kolding 621
## 68 Vejle 630
## 69 Herning 657
## 70 Holstebro 661
## 71 Lemvig 665
## 72 Struer 671
## 73 Syddjurs 706
## 74 Norddjurs 707
## 75 Favrskov 710
## 76 Odder 727
## 77 Randers 730
## 78 Silkeborg 740
## 79 Samsø 741
## 80 Skanderborg 746
## 81 Aarhus 751
## 82 Ikast-Brande 756
## 83 Ringkøbing-Skjern 760
## 84 Hedensted 766
## 85 Morsø 773
## 86 Skive 779
## 87 Thisted 787
## 88 Viborg 791
## 89 Brønderslev 810
## 90 Frederikshavn 813
## 91 Vesthimmerlands 820
## 92 Læsø 825
## 93 Rebild 840
## 94 Mariagerfjord 846
## 95 Jammerbugt 849
## 96 Aalborg 851
## 97 Hjørring 860 slot(foo,"data")## muni_name ID_muni
## 0 København 101
## 1 Frederiksberg 147
## 2 Ballerup 151
## 3 Brøndby 153
## 4 Dragør 155
## 5 Gentofte 157
## 6 Gladsaxe 159
## 7 Glostrup 161
## 8 Herlev 163
## 9 Albertslund 165
## 10 Hvidovre 167
## 11 Høje Taastrup 169
## 12 Lyngby-Taarbæk 173
## 13 Rødovre 175
## 14 Ishøj 183
## 15 TÃ¥rnby 185
## 16 Vallensbæk 187
## 17 Furesø 190
## 18 Allerød 201
## 19 Fredensborg 210
## 20 Helsingør 217
## 21 Hillerød 219
## 22 Hørsholm 223
## 23 Rudersdal 230
## 24 Egedal 240
## 25 Frederikssund 250
## 26 Greve 253
## 27 Køge 259
## 28 Halsnæs 260
## 29 Roskilde 265
## 30 Solrød 269
## 31 Gribskov 270
## 32 Odsherred 306
## 33 Holbæk 316
## 34 Faxe 320
## 35 Kalundborg 326
## 36 Ringsted 329
## 37 Slagelse 330
## 38 Stevns 336
## 39 Sorø 340
## 40 Lejre 350
## 41 Lolland 360
## 42 Næstved 370
## 43 Guldborgsund 376
## 44 Vordingborg 390
## 45 Bornholm 400
## 46 Middelfart 410
## 47 Assens 420
## 48 Faaborg-Midtfyn 430
## 49 Kerteminde 440
## 50 Nyborg 450
## 51 Odense 461
## 52 Svendborg 479
## 53 Nordfyns 480
## 54 Langeland 482
## 55 Ærø 492
## 56 Haderslev 510
## 57 Billund 530
## 58 Sønderborg 540
## 59 Tønder 550
## 60 Esbjerg 561
## 61 Fanø 563
## 62 Varde 573
## 63 Vejen 575
## 64 Aabenraa 580
## 65 Fredericia 607
## 66 Horsens 615
## 67 Kolding 621
## 68 Vejle 630
## 69 Herning 657
## 70 Holstebro 661
## 71 Lemvig 665
## 72 Struer 671
## 73 Syddjurs 706
## 74 Norddjurs 707
## 75 Favrskov 710
## 76 Odder 727
## 77 Randers 730
## 78 Silkeborg 740
## 79 Samsø 741
## 80 Skanderborg 746
## 81 Aarhus 751
## 82 Ikast-Brande 756
## 83 Ringkøbing-Skjern 760
## 84 Hedensted 766
## 85 Morsø 773
## 86 Skive 779
## 87 Thisted 787
## 88 Viborg 791
## 89 Brønderslev 810
## 90 Frederikshavn 813
## 91 Vesthimmerlands 820
## 92 Læsø 825
## 93 Rebild 840
## 94 Mariagerfjord 846
## 95 Jammerbugt 849
## 96 Aalborg 851
## 97 Hjørring 860Merging data to the sf-object
The sf-object is a data.frame. Each row is an sfg. This is what would normally be thought of as an observation. So one row one observation. In this case each observation is associated with an area in space defined by the polygon. Each polygon have an identifying number.
Data is merged using the merge command merge(x,y), where x and y are data.frames. Personally I prefer to work with data.tables but the sf-object cannot be transformed to a data.table without loosing its sf-class. This means that it cannot be plotted anymore using the sf-plot method.
# Load data
path_shape <- "C://Users//np83zg//OneDrive - Aalborg Universitet//Skrivebord//cars-master//rspace_osm//Learning SF//shapes"
setwd(path_shape)
mun_shp <- st_read("municipalities_DK.shp",options = "ENCODING=UTF-8")## options: ENCODING=UTF-8
## Reading layer `municipalities_DK' from data source `C:\Users\np83zg\OneDrive - Aalborg Universitet\Skrivebord\cars-master\rspace_osm\Learning SF\shapes\municipalities_DK.shp' using driver `ESRI Shapefile'
## Simple feature collection with 98 features and 2 fields
## Geometry type: MULTIPOLYGON
## Dimension: XYZ
## Bounding box: xmin: 441676 ymin: 6049782 xmax: 892787 ymax: 6402217
## z_range: zmin: -999 zmax: 0
## Projected CRS: ETRS89 / UTM zone 32N library(data.table)## Warning: package 'data.table' was built under R version 4.0.5 kom_names <- read.table("kom_names.txt",header=TRUE)
kom_names <- as.data.table(kom_names)
kom_names[,kom_ny:=NULL]
setkey(kom_names,"kom")
plot(mun_shp)
mun_shp <- as.data.table(mun_shp)
class(mun_shp)## [1] "data.table" "data.frame" mun_shp <- merge(mun_shp,kom_names,by.x="ID_muni",by.y="kom")
mun_shp[,muni_name:=NULL]You can plot before transforming to data.table but not after. For merging you have to transform if you wanna depend on data.table’s data-management capabilities. However, then the object is no longer sf and cannot be plotted unless transformed back to sf. But it not necessary to transform to merge because you simply merge data.frames.
Load shape, merge and plot
path_main <- "C://Users//np83zg//OneDrive - Aalborg Universitet//Skrivebord//cars-master//rspace_osm//Learning SF"
path_shape <- "C://Users//np83zg//OneDrive - Aalborg Universitet//Skrivebord//cars-master//rspace_osm//Learning SF//shapes"
setwd(path_shape)
mun_shp <- st_read("municipalities_DK.shp",options = "ENCODING=UTF-8")## options: ENCODING=UTF-8
## Reading layer `municipalities_DK' from data source `C:\Users\np83zg\OneDrive - Aalborg Universitet\Skrivebord\cars-master\rspace_osm\Learning SF\shapes\municipalities_DK.shp' using driver `ESRI Shapefile'
## Simple feature collection with 98 features and 2 fields
## Geometry type: MULTIPOLYGON
## Dimension: XYZ
## Bounding box: xmin: 441676 ymin: 6049782 xmax: 892787 ymax: 6402217
## z_range: zmin: -999 zmax: 0
## Projected CRS: ETRS89 / UTM zone 32N library(data.table)
kom_names <- read.table("kom_names.txt",header=TRUE)
kom_names <- as.data.table(kom_names)
kom_names[,kom_ny:=NULL]
setkey(kom_names,"kom")
mun_shp <- merge(mun_shp,kom_names,by.x="ID_muni",by.y="kom")
library(tmap)## Warning: package 'tmap' was built under R version 4.0.5 mun_shape <- mun_shp
mun_shape$dn <- abs(rnorm(98))
tm_shape(mun_shape) + tm_polygons("dn",palette="RdBu",title="Population growth rate") +
tm_shape(mun_shape) + tm_borders("white", lwd = 1.2) 
setwd(path_main)SF-methods
There are several methods for working with sf objects
methods(class = "sf")## [1] $<- [ [[<-
## [4] aggregate as.data.frame cbind
## [7] coerce dbDataType dbWriteTable
## [10] filter identify initialize
## [13] merge plot print
## [16] rbind show slotsFromS3
## [19] st_agr st_agr<- st_area
## [22] st_as_s2 st_as_sf st_bbox
## [25] st_boundary st_buffer st_cast
## [28] st_centroid st_collection_extract st_convex_hull
## [31] st_coordinates st_crop st_crs
## [34] st_crs<- st_difference st_filter
## [37] st_geometry st_geometry<- st_inscribed_circle
## [40] st_interpolate_aw st_intersection st_intersects
## [43] st_is st_is_valid st_join
## [46] st_line_merge st_m_range st_make_valid
## [49] st_nearest_points st_node st_normalize
## [52] st_point_on_surface st_polygonize st_precision
## [55] st_reverse st_sample st_segmentize
## [58] st_set_precision st_shift_longitude st_simplify
## [61] st_snap st_sym_difference st_transform
## [64] st_triangulate st_union st_voronoi
## [67] st_wrap_dateline st_write st_z_range
## [70] st_zm transform
## see '?methods' for accessing help and source code