Exercise 2_1: Import Vector Data

Introduction to Geospatial Techniques for Social Scientists in R

Author

Stefan Jünger, Anne-Kathrin Stroppe, Dennis Abel

Exercises

Let’s get started with some hands-on exercises and exploring some data! This exercise focuses on importing different vector data geometries from various sources and getting used to simple feature data tables.

In the folder data, you can find the data files prepped for all the exercises which rely on preexisting data sets. However, if you like to play around with your data, feel free to do so!

🏋 Exercise 1

Load the administrative borders of Germany. For this purpose, import the shapefile VG250_STA.shp from the respective folder and assign it to an object named germany.

Plot the newly created object.

💡 Tip

Don’t forget to load the packages sf and set your working directory before starting this exercise.

The R console output already gives you some information on the data you just loaded. You also see that the object germany appeared in your environment. Let’s explore the data set a little bit closely.

🏋 Exercise 2

Explore the data and answer the following questions:

What is the geometry type?
What is the ID of the coordinate reference system?
What is the name of the column containing the information on geometries?
Which attributes are assigned to the geometries?

💡 Tip

You can either extract the information directly from the output or problematically with functions from base R and the sf package (you may want to have a look, e.g., at the help files ?sf::st_geometry or sf::st_crs()).

🏋 Exercise 3

Do you have an idea why there are 11 observations instead of only 1, even though we expected a shapefile with the polygon of Germany? Create a new object that contains only one observation filtered by the variable GF == 4. The data provider, the Federal Agency for Cartography and Geodesy (BKG), has named this variable, which stands for “Geofactor” (German: “Geofaktor”). Maybe you can figure out what this category means after plotting the data.

💡 Tip

You can apply standard R procedures to filter data, such as dplry’s filter() command.

Solutions

✅ Solution 1

library(sf)

germany <- sf::read_sf("./data/VG250_STA.shp") 

plot(germany)

Warning: Paket 'sf' wurde unter R Version 4.4.3 erstellt

Linking to GEOS 3.13.0, GDAL 3.10.1, PROJ 9.5.1; sf_use_s2() is TRUE

Warning: plotting the first 10 out of 23 attributes; use max.plot = 23 to plot
all

✅ Solution 2

sf::st_geometry(germany) # Multipolygon: Set of several Polygons. Each row is a polygon.

Geometry set for 11 features 
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 280371.1 ymin: 5235856 xmax: 921292.4 ymax: 6106244
Projected CRS: ETRS89 / UTM zone 32N
First 5 geometries:

MULTIPOLYGON (((609387.6 5267931, 609423.3 5267...

MULTIPOLYGON (((492810 6029412, 492821.2 602927...

MULTIPOLYGON (((816809.9 5978628, 816854.5 5978...

MULTIPOLYGON (((502125 5972301, 502668.3 597211...

MULTIPOLYGON (((500425.1 5285873, 500485.1 5285...

attr(germany, "sf_column") # geometry

[1] "geometry"

sf::st_crs(germany) # ETRS89, ID = EPSG:25832

Coordinate Reference System:
  User input: ETRS89 / UTM zone 32N 
  wkt:
PROJCRS["ETRS89 / UTM zone 32N",
    BASEGEOGCRS["ETRS89",
        ENSEMBLE["European Terrestrial Reference System 1989 ensemble",
            MEMBER["European Terrestrial Reference Frame 1989"],
            MEMBER["European Terrestrial Reference Frame 1990"],
            MEMBER["European Terrestrial Reference Frame 1991"],
            MEMBER["European Terrestrial Reference Frame 1992"],
            MEMBER["European Terrestrial Reference Frame 1993"],
            MEMBER["European Terrestrial Reference Frame 1994"],
            MEMBER["European Terrestrial Reference Frame 1996"],
            MEMBER["European Terrestrial Reference Frame 1997"],
            MEMBER["European Terrestrial Reference Frame 2000"],
            MEMBER["European Terrestrial Reference Frame 2005"],
            MEMBER["European Terrestrial Reference Frame 2014"],
            MEMBER["European Terrestrial Reference Frame 2020"],
            ELLIPSOID["GRS 1980",6378137,298.257222101,
                LENGTHUNIT["metre",1]],
            ENSEMBLEACCURACY[0.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; Denmark - onshore and offshore; Germany - onshore and offshore; Italy - onshore and offshore; Norway including Svalbard - onshore and offshore; Spain - offshore."],
        BBOX[36.53,6,84.01,12.01]],
    USAGE[
        SCOPE["Pan-European conformal mapping at scales larger than 1:500,000."],
        AREA["Europe between 6°E and 12°E and approximately 36°30'N to 84°N."],
        BBOX[36.53,6,84.01,12.01]],
    ID["EPSG",25832]]

names(germany) # Name, Area Type, Area Size, Population Size ,...

 [1] "ADE"      "GF"       "BSG"      "ARS"      "AGS"      "SDV_ARS" 
 [7] "GEN"      "BEZ"      "IBZ"      "BEM"      "NBD"      "SN_L"    
[13] "SN_R"     "SN_K"     "SN_V1"    "SN_V2"    "SN_G"     "FK_S3"   
[19] "NUTS"     "ARS_0"    "AGS_0"    "WSK"      "DEBKG_ID" "geometry"

✅ Solution 3

# It seems like the shapefile contains not only Germany's land area but also coastal areas and lakes.
# filter 
germany_new <-
  germany |> 
  dplyr::filter(GF == 4)

plot(germany_new)

Warning: plotting the first 10 out of 23 attributes; use max.plot = 23 to plot
all