Topographic Database

Krzysztof Dyba, Jakub Nowosad

Definition

Topographic Database (pl. Baza Danych Obiektów Topograficznych) is a vector (object) database containing the spatial location of topographic objects with their characteristics for Poland. The content and detail of the database correspond to the topographic map in the scale 1:10000. The thematic scope includes information on water network, communication network, land cover, buildings and technical structures, utility infrastructure, land use, protected areas, territorial division units, and other objects. The database is available in the Geography Markup Language (GML) format. The source of its data comes from:

Purpose

The purpose of this vignette is to perform spatial operations on vector data from the Topographic Database. We focus on four cases, taking into account different types of geometry, i.e. point, line, and polygon, and their attributes. Also, we show how they can be visualized.

Analysis

# attach packages
library(sf)
library(rgugik)

Our analysis area is the bieszczadzki county located in the Subcarpathian (podkarpackie) voivodeship. It is the farthest south area in Poland that also has the lowest population density (19 people on km2).

Database

We start by downloading the topographic database for our county using the topodb_download() function.

# 22.4 MB
topodb_download("bieszczadzki", outdir = "./data")

If you run into any problem with the download, remember that you can pass another download method from download.file() as a function argument.

topodb_download(req_df, outdir = "./data", method = "wget")

The downloaded database consists of many files in the GML format. All the data necessary for analyzes can be found in the data/PL.PZGiK.332.1801/BDOT10k/ location. A brief description of the structure of this database can be found here.

First, let’s load the file with administrative units (“PL.PZGiK.332.1801__OT_ADJA_A.xml”) using the sf package and its read_sf() function.

territory = read_sf("data/PL.PZGiK.332.1801/BDOT10k/PL.PZGiK.332.1801__OT_ADJA_A.xml")

The file contains administrative units at various levels, let’s choose the lowest level, i.e. communes. There are three types of communes in this dataset, specified in the rodzaj column: urban (GM), rural (GW), and urban-rural (Gmw). Let’s select them.

communes = territory[territory$rodzaj %in% c("GM", "GW", "Gmw"), "rodzaj"]
table(communes$rodzaj)
## 
## Gmw  GW 
##   1   2

We can see that the bieszczadzki county consists of two rural communes and one urban-rural commune. Let’s visualize it.

plot(communes, axes = TRUE, main = "Bieszczadzki county")

Lengths and categories of roads

In the first task, we calculate the lengths of roads, taking into account their categories. Road data can be found in the “PL.PZGiK.332.1801__OT_SKDR_L.xml” file.

roads = read_sf("data/PL.PZGiK.332.1801/BDOT10k/PL.PZGiK.332.1801__OT_SKDR_L.xml")

Let’s plot them. We use the plot() function again, but this time we combine the two layers into one image. The first layer (in the background) contains roads and to add another layer, we have to set the argument reset to "FALSE". Then we can add a second layer with the territory borders by setting the add argument to "TRUE".

plot(roads["katZarzadzania"], main = "Road category", reset = FALSE)
plot(st_geometry(territory), add = TRUE)

We have six road categories related to the managing entity. Those are: national (K), voivodeship (W), county (P), communal (G), institutional (Z), and other (I) roads.

We use the st_length() function to find the length of each object in the table. Next, we create a data frame consisting of the road category and its length. Then we aggregate this data frame and calculate the sum of the lengths for each category.

length_roads = st_length(roads)
length_roads = data.frame(length = length_roads,
                          class = as.factor(roads$katZarzadzania))
length_roads = aggregate(length ~ class, data = length_roads, FUN = sum)

The results are given in meters - let’s convert them into kilometers.

# convert to [km]
length_roads$length = units::set_units(length_roads$length, "km")

Let’s also change the names of the categories.

road_class = c("communal", "other", "national", "county", "voivodeship",
               "institutional")
levels(length_roads$class) = road_class

Now we can see the results. The other type of roads dominates that consists mainly of non-public roads.

length_roads
##           class          length
## 1      communal  186.17329 [km]
## 2         other 2819.55806 [km]
## 3      national   19.24996 [km]
## 4        county  187.17870 [km]
## 5   voivodeship  105.67567 [km]
## 6 institutional   10.44893 [km]

We can also calculate the total length of the roads in this area.

sum(length_roads$length)
## 3328.285 [km]

The result is about 3328.285 km. Another aspect of the data that we can investigate is the density of the road network. We need to calculate the total area first, and then divide the total road length by the total area.

communes_area = sum(st_area(communes))
communes_area = units::set_units(communes_area, "km2")
density = sum(length_roads$length)/communes_area
density = units::set_units(density, "km/km2")
density
## 2.921031 [km/km2]

The road density is about 2.92 km/km2.

Roads through the rivers

Another dataset included in the Topographic Database contains rivers for this area (“PL.PZGiK.332.1801__OT_SWRS_L.xml”).

rivers = read_sf("data/PL.PZGiK.332.1801/BDOT10k/PL.PZGiK.332.1801__OT_SWRS_L.xml")
rivers = rivers[rivers$rodzaj == "Rz", ] # select only rivers

Rivers are divided into smaller sections with different parameters, such as river width or data source. Let’s merge sections from the same rivers into a single feature (geometry) - we can use the attribute with an ID (idMPHP) for that purpose. We can also give each river a category number by creating a sequence from 1 to n using the seq_len function.

rivers = aggregate(rivers[, c("geometry", "idMPHP")],
                   list(rivers$idMPHP),
                   sum)
rivers$idMPHP = seq_len(length(unique(rivers$idMPHP)))
rivers$idMPHP = as.factor(rivers$idMPHP)

Let’s visualize the rivers’ courses.

plot(rivers["idMPHP"], main = "Rivers", reset = FALSE)
plot(st_geometry(territory), add = TRUE)

With rivers and roads, we can designate points of intersection that symbolize bridges and crossings. We can use the st_intersection() function for this.

bridges = st_geometry(st_intersection(rivers, roads))
length(bridges)
## [1] 81

We get 81 such points. Let’s plot them.

# use 'dev.off()' to reset previous plot
plot(st_geometry(rivers), main = "Bridges and crossings", col = "blue")
plot(st_geometry(territory), add = TRUE)
plot(bridges, add = TRUE, pch = 20)