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:
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.
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).
We start by downloading the topographic database for our county using
the topodb_download()
function.
If you run into any problem with the download, remember that you can
pass another download method from download.file()
as a
function argument.
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.
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.
##
## 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.
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.
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.
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.
## 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.
## 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.
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.
With rivers and roads, we can designate points of intersection that
symbolize bridges and crossings. We can use the
st_intersection()
function for this.
## [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)