# install.packages("sf")
library(sf)
library(tidyverse)R Tutorial: sf
Dieses Tutorial basiert auf Arbeiten von BMARSYS-Studierenden im begleitenden Seminar zum DS2-Modul und wurde redaktionell überarbeitet.
Das R-Paket sf (Simple Features) ist das moderne Standardpaket für räumliche Vektordaten in R. Es hat die älteren Pakete sp, rgdal und rgeos abgelöst und integriert sich nahtlos in das tidyverse — insbesondere in dplyr für Datenmanipulation und ggplot2 für Visualisierung. „Simple Features” ist dabei ein internationaler Standard (ISO 19125) für die Repräsentation von Geometrien wie Punkten, Linien und Polygonen.
Installation
sf-Objekte verstehen
sf verwendet drei hierarchisch aufgebaute Klassen, die aufeinander aufbauen:
- sfg (Simple Feature Geometry) — Die Geometrie eines einzelnen Features, z.B. ein Punkt (
POINT), eine Linie (LINESTRING) oder ein Polygon (POLYGON) - sfc (Simple Feature Column) — Eine Listenspalte, die mehrere sfg-Objekte zusammenfasst, also die Geometrien aller Features in einem Datensatz. Hier wird auch das Koordinatenreferenzsystem (CRS) gespeichert.
- sf — Ein normaler Data Frame mit einer zusätzlichen sfc-Geometriespalte (standardmäßig
geometrygenannt). Das ist die Klasse, mit der man im Alltag arbeitet.
Alle sf-Funktionen beginnen mit dem Präfix st_ (für „spatial type”). Eine Übersicht aller verfügbaren Methoden für sf-Objekte:
methods(class = "sf") [1] [ [[<-
[3] [<- $<-
[5] aggregate anti_join
[7] arrange as.data.frame
[9] cbind coerce
[11] dbDataType dbWriteTable
[13] distinct dplyr_reconstruct
[15] drop_na duplicated
[17] filter full_join
[19] gather group_by
[21] group_split identify
[23] initialize inner_join
[25] left_join merge
[27] mutate nest
[29] pivot_longer pivot_wider
[31] plot points
[33] print rbind
[35] rename_with rename
[37] right_join rowwise
[39] sample_frac sample_n
[41] select semi_join
[43] separate_rows separate
[45] show slice
[47] slotsFromS3 spread
[49] st_agr st_agr<-
[51] st_area st_as_s2
[53] st_as_sf st_as_sfc
[55] st_bbox st_boundary
[57] st_break_antimeridian st_buffer
[59] st_cast st_centroid
[61] st_collection_extract st_concave_hull
[63] st_convex_hull st_coordinates
[65] st_crop st_crs
[67] st_crs<- st_difference
[69] st_drop_geometry st_exterior_ring
[71] st_filter st_geometry
[73] st_geometry<- st_inscribed_circle
[75] st_interpolate_aw st_intersection
[77] st_intersects st_is_full
[79] st_is_valid st_is
[81] st_join st_line_merge
[83] st_m_range st_make_valid
[85] st_minimum_bounding_circle st_minimum_rotated_rectangle
[87] st_nearest_points st_node
[89] st_normalize st_point_on_surface
[91] st_polygonize st_precision
[93] st_reverse st_sample
[95] st_segmentize st_set_precision
[97] st_shift_longitude st_simplify
[99] st_snap st_sym_difference
[101] st_transform st_triangulate_constrained
[103] st_triangulate st_union
[105] st_voronoi st_wrap_dateline
[107] st_write st_z_range
[109] st_zm summarise
[111] text transform
[113] transmute ungroup
[115] unite unnest
see '?methods' for accessing help and source codeDaten einlesen und erkunden
sf kann eine Vielzahl räumlicher Formate lesen — darunter Shapefiles (.shp), GeoJSON (.geojson), GeoPackage (.gpkg) und viele weitere über die GDAL-Bibliothek. Die zentrale Funktion ist st_read(). Für dieses Tutorial verwenden wir den im sf-Paket enthaltenen Beispieldatensatz nc (Bezirke in North Carolina mit demographischen Daten):
# Beispieldatensatz: Bezirke in North Carolina (im sf-Paket enthalten)
nc <- st_read(system.file("shape/nc.shp", package = "sf"))Reading layer `nc' from data source
`/Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/sf/shape/nc.shp'
using driver `ESRI Shapefile'
Simple feature collection with 100 features and 14 fields
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: -84.32385 ymin: 33.88199 xmax: -75.45698 ymax: 36.58965
Geodetic CRS: NAD27# Klasse prüfen — sf erweitert die data.frame-Klasse
class(nc)[1] "sf" "data.frame"# Spalten anzeigen
colnames(nc) [1] "AREA" "PERIMETER" "CNTY_" "CNTY_ID" "NAME" "FIPS"
[7] "FIPSNO" "CRESS_ID" "BIR74" "SID74" "NWBIR74" "BIR79"
[13] "SID79" "NWBIR79" "geometry" # Erste Zeilen
head(nc, 5)Simple feature collection with 5 features and 14 fields
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: -81.74107 ymin: 36.07282 xmax: -75.77316 ymax: 36.58965
Geodetic CRS: NAD27
AREA PERIMETER CNTY_ CNTY_ID NAME FIPS FIPSNO CRESS_ID BIR74 SID74
1 0.114 1.442 1825 1825 Ashe 37009 37009 5 1091 1
2 0.061 1.231 1827 1827 Alleghany 37005 37005 3 487 0
3 0.143 1.630 1828 1828 Surry 37171 37171 86 3188 5
4 0.070 2.968 1831 1831 Currituck 37053 37053 27 508 1
5 0.153 2.206 1832 1832 Northampton 37131 37131 66 1421 9
NWBIR74 BIR79 SID79 NWBIR79 geometry
1 10 1364 0 19 MULTIPOLYGON (((-81.47276 3...
2 10 542 3 12 MULTIPOLYGON (((-81.23989 3...
3 208 3616 6 260 MULTIPOLYGON (((-80.45634 3...
4 123 830 2 145 MULTIPOLYGON (((-76.00897 3...
5 1066 1606 3 1197 MULTIPOLYGON (((-77.21767 3...Wie man sieht, ist ein sf-Objekt im Wesentlichen ein Data Frame mit einer zusätzlichen Spalte geometry, die die räumliche Information enthält. Alle normalen Operationen auf Data Frames funktionieren weiterhin.
Schnellvisualisierung
Die base-R-Funktion plot() erstellt für sf-Objekte automatisch eine Kartenansicht. Standardmäßig werden bis zu 6 Attribute als separate Karten dargestellt. Durch Angabe eines einzelnen Attributs in eckigen Klammern wird nur eine Karte erzeugt:
# Standardmäßig werden bis zu 6 Attribute geplottet
plot(nc)
# Einzelnes Attribut
plot(nc["AREA"])
Koordinatenreferenzsysteme
Ein Koordinatenreferenzsystem (CRS) definiert, wie Koordinaten auf der Erdoberfläche zu interpretieren sind. Es gibt geographische CRS (Länge/Breite in Grad) und projizierte CRS (x/y in Metern). Die Wahl des CRS beeinflusst, wie Entfernungen und Flächen berechnet werden.
CRS abfragen
st_crs(nc)Coordinate Reference System:
User input: NAD27
wkt:
GEOGCRS["NAD27",
DATUM["North American Datum 1927",
ELLIPSOID["Clarke 1866",6378206.4,294.978698213898,
LENGTHUNIT["metre",1]]],
PRIMEM["Greenwich",0,
ANGLEUNIT["degree",0.0174532925199433]],
CS[ellipsoidal,2],
AXIS["latitude",north,
ORDER[1],
ANGLEUNIT["degree",0.0174532925199433]],
AXIS["longitude",east,
ORDER[2],
ANGLEUNIT["degree",0.0174532925199433]],
ID["EPSG",4267]]CRS transformieren
Mit st_transform() werden Daten in ein anderes Koordinatensystem umprojiziert. Dies ist z.B. nötig, wenn man verschiedene Datensätze mit unterschiedlichen CRS kombinieren möchte oder wenn man Flächen/Entfernungen in Metern berechnen will:
# Transformation in ETRS89/UTM Zone 32N (EPSG: 25832)
# nc_utm <- st_transform(nc, 25832)
# Häufige EPSG-Codes:
# 4326 — WGS84 (Länge/Breite) — der GPS-Standard
# 25832 — ETRS89/UTM Zone 32N (Europa)
# 3857 — WGS84 Pseudo-Mercator (Google Maps, OpenStreetMap)CRS in ggplot2
In ggplot2 kann das CRS direkt in coord_sf() gesetzt werden, um die Karte in einer bestimmten Projektion darzustellen — die zugrunde liegenden Daten werden dabei nicht verändert:
ggplot(nc) +
geom_sf(aes(fill = AREA)) +
scale_fill_viridis_c() +
coord_sf(crs = st_crs(3857)) +
ggtitle("WGS 84/Pseudo-Mercator") +
theme_bw()
Räumliche Operationen
sf bietet zahlreiche räumliche Operationen, die direkt auf sf-Objekte angewendet werden können. Die wichtigsten sind hier als Übersicht zusammengestellt. Die auskommentierten Befehle sind Beispiele, die zwei Layer (a, b) benötigen:
# Verschneidung (Intersection) — gemeinsamer Bereich zweier Layer
# result <- st_intersection(a, b)
# Vereinigung (Union) — Zusammenführen zu einem Objekt
# result <- st_union(a, b)
# Pufferung (Buffer) — z.B. 1000 Meter um jedes Feature
# result <- st_buffer(data, dist = 1000)
# Differenz — was in a, aber nicht in b liegt
# result <- st_difference(a, b)
# Schwerpunkt (Centroid) — geometrischer Mittelpunkt jedes Polygons
# centroids <- st_centroid(nc)
# Fläche berechnen (Ergebnis in m², abhängig vom CRS)
st_area(nc)Units: [m^2]
[1] 1137107793 610916077 1423145355 694378925 1520366979 967504822
[7] 615794941 903423919 1179065710 1232475139 1136017416 1524295167
[13] 1426763054 1085709751 718024778 1893655681 524303669 1986581059
[19] 812132036 626215554 439637846 640597398 863142124 1276325061
[25] 1083947009 1697657775 1109799786 1800353048 1036247721 770426970
[31] 1422972995 585145178 1311460371 1224942117 800163805 1186288078
[37] 2194374294 1179004039 1550151186 690514844 665066784 1457728244
[43] 1340416729 1005633561 988219530 1163804357 2019609428 1810365923
[49] 944348527 1350014516 1685059736 1068120639 1691385005 2082034143
[55] 1447025244 943796075 2045470574 1420873777 707648814 653349704
[61] 1471057561 1436128964 1550970115 1186032312 788508058 1265459073
[67] 1829451696 1447903974 918204712 1312725482 1043733633 961860879
[73] 781909574 1046090580 986760532 917758923 601335294 1321974824
[79] 2438120829 833576485 1210382282 1738664778 1228776807 1648541762
[85] 1400697543 995179656 1678005426 2072031752 1228366621 519232890
[91] 1785013769 808690576 1978885855 2439935278 1264198838 2289052992
[97] 2181566551 2450830549 430798470 2166454052# Begrenzungsrahmen (Bounding Box)
st_bbox(nc) xmin ymin xmax ymax
-84.32385 33.88199 -75.45698 36.58965
Tipp
Für Flächen- und Entfernungsberechnungen sollte ein projiziertes CRS (z.B. UTM) verwendet werden, da geographische Koordinaten (Grad) zu ungenauen Ergebnissen führen.
Visualisierung mit ggplot2
Die tidyverse-Integration von sf macht die Visualisierung mit ggplot2 besonders einfach. Die Funktion geom_sf() erkennt automatisch die Geometriespalte — aes(x = ..., y = ...) muss nicht angegeben werden. Farben, Rahmen und andere Aesthetics können wie gewohnt über aes() gesteuert werden:
ggplot(nc) +
geom_sf(aes(fill = AREA)) +
scale_fill_viridis_c("Fläche") +
theme_bw()
Mehrere Layer kombinieren
Wie in ggplot2 üblich, können mehrere geom_sf()-Layer übereinandergelegt werden. Jeder Layer kann aus einem anderen sf-Objekt stammen. Im folgenden Beispiel werden die Bezirksgrenzen mit ihren Schwerpunkten (Centroids) kombiniert:
# Schwerpunkte berechnen
nc_centroids <- st_centroid(nc)
ggplot() +
geom_sf(data = nc, aes(fill = AREA)) +
geom_sf(data = nc_centroids, color = "red", size = 0.5) +
scale_fill_viridis_c("Fläche") +
theme_bw()
Daten manipulieren mit dplyr
Ein großer Vorteil von sf gegenüber den Vorgängerpaketen: Da sf-Objekte Data Frames sind, funktionieren alle dplyr-Verben (filter(), mutate(), select(), group_by(), summarise() etc.) direkt. Die Geometriespalte wird dabei automatisch beibehalten — auch wenn sie nicht explizit ausgewählt wird. Hier ein Beispiel, das neue Variablen berechnet, filtert und das Ergebnis direkt an ggplot2 weitergibt:
nc %>%
mutate(density = BIR74 / AREA) %>%
select(NAME, density) %>%
filter(density > 1000) %>%
ggplot() +
geom_sf(aes(fill = density)) +
scale_fill_viridis_c() +
theme_bw()
Tipp
summarise() auf gruppierte sf-Objekte führt automatisch ein st_union() der Geometrien durch — so lassen sich z.B. Bezirke zu Regionen zusammenfassen.