Author

Dominic Royé

Published

November 1, 2019

Modified

December 30, 2024

The General Directorate for the Cadastre of Spain has spatial information of the all buildings except for the Basque Country and Navarra. This data set is part of the implementation of INSPIRE, the Space Information Infrastructure in Europe. More information can be found here. We will use the links (urls) in ATOM format, which is an RSS type for web feeds, allowing us to obtain the download link for each municipality. Since 2022, a package to access the API is directly available CatastRo.

Packages

Package Description
tidyverse Collection of packages (visualization, manipulation): ggplot2, dplyr, purrr, etc.
sf Simple Feature: import, export and manipulate vector data
CatastRo Catastro Spain API
tmap Easy creation of thematic maps
classInt Create univariate class intervals
# install the packages if necessary

if (!require("tidyverse")) install.packages("tidyverse")
if (!require("CatastRo")) install.packages("CatastRo")
if (!require("sf")) install.packages("sf")
if (!require("tmap")) install.packages("tmap")
if (!require("classInt")) install.packages("classInt")

# load packages

library(sf)
library(CatastRo)
library(tidyverse)
library(classInt)
library(tmap)

Access the data

To import the building data we use the catr_atom_get_buildings() function.

buildings_val <- catr_atom_get_buildings("Valencia", to = "Valencia")

buildings_val[,1:5]
Simple feature collection with 36346 features and 5 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 720570.9 ymin: 4351286 xmax: 734981.9 ymax: 4382906
Projected CRS: ETRS89 / UTM zone 30N
First 10 features:
                      gml_id               lowerCorner              upperCorner
1  ES.SDGC.BU.000100100YJ27F  725589.5805 4377132.9425 725605.5005 4377146.6725
2  ES.SDGC.BU.000100100YJ28A  720570.9393 4382157.3445  720646.359 4382259.3929
3  ES.SDGC.BU.000100100YJ36A  730582.5135 4362033.4165  730589.793 4362043.9665
4  ES.SDGC.BU.000100200YJ27F    725711.36 4377093.8715     725745.63 4377124.46
5  ES.SDGC.BU.000100200YJ28A   720851.8475 4381937.212  720866.847 4381955.6425
6  ES.SDGC.BU.000100300YJ27F      725336.89 4376989.78     725386.83 4377084.18
7  ES.SDGC.BU.000100300YJ36C  729967.6535 4364320.2035 729979.9435 4364324.2535
8  ES.SDGC.BU.000100400YJ27F  724956.7325 4376895.1345  724971.625 4376909.0335
9  ES.SDGC.BU.000100400YJ36C   729953.5525 4364259.855  730039.851 4364319.6035
10 ES.SDGC.BU.000200100YJ26F 726653.43835 4367240.6706  726673.662 4367258.1575
   beginLifespanVersion conditionOfConstruction                       geometry
1   2008-10-20T00:00:00              functional MULTIPOLYGON (((725594.9 43...
2   2022-03-10T00:00:00              functional MULTIPOLYGON (((720582.3 43...
3   2006-01-18T00:00:00              functional MULTIPOLYGON (((730586.1 43...
4   2016-10-04T00:00:00              functional MULTIPOLYGON (((725718.3 43...
5   2008-10-20T00:00:00              functional MULTIPOLYGON (((720860.3 43...
6   2016-03-15T00:00:00              functional MULTIPOLYGON (((725376.5 43...
7   2008-10-20T00:00:00                    ruin MULTIPOLYGON (((729967.7 43...
8   2016-03-15T00:00:00              functional MULTIPOLYGON (((724958.6 43...
9   2008-10-20T00:00:00              functional MULTIPOLYGON (((729984.8 43...
10  2017-01-31T00:00:00              functional MULTIPOLYGON (((726661.4 43...

Data preparation

We only have to convert the column of the construction year (beginning) into a Date class. The date column contains some dates in --01-01 format, which does not correspond to any recognizable date. Therefore, we replace the first - with 0000.

buildings_val <- mutate(buildings_val,
  beginning = str_replace(beginning, "^-", "0000") |>
    ymd_hms() |> as_date()
)
Warning: There was 1 warning in `stopifnot()`.
ℹ In argument: `beginning = as_date(ymd_hms(str_replace(beginning, "^-",
  "0000")))`.
Caused by warning:
!  6 failed to parse.

Distribution chart

Before creating the maps of the construction years, which will reflect urban growth, we will make a graph of distribution of the beginning variable. We can clearly identify periods of urban expansion. We will use the ggplot2 package with the geometry of geom_density() for this purpose.

# limit the period after 1750
filter(buildings_val, beginning >= "1750-01-01") |>
  ggplot(aes(beginning)) +
  geom_density(fill = "#2166ac", alpha = 0.7) +
  scale_x_date(
    date_breaks = "20 year",
    date_labels = "%Y"
  ) +
  labs(y = NULL, x = NULL, title = "Evolution of urban development") +
  theme_minimal(base_family = "Montserrat") +
  theme(panel.grid.minor = element_blank())
Warning in grid.Call(C_stringMetric, as.graphicsAnnot(x$label)): font family
not found in Windows font database
Warning in grid.Call(C_stringMetric, as.graphicsAnnot(x$label)): font family
not found in Windows font database
Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
font family not found in Windows font database
Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
font family not found in Windows font database

Buffer of 2,5 km for Valencia

To visualize better the distribution of urban growth, we limit the map to a radius of 2.5 km from the city center. Therefore, we use the geocode_OSM() function of the tmaptools package to obtain the coordinates of Valencia in class sf. Then we project the points to the system we use for the buildings (EPSG: 25830). The st_crs() function returns the coordinate system of a spatial object sf. Finally, we create with the function st_buffer() a buffer with 2500 m and the intersection with our building data. It is also possible to create a buffer in the form of a rectangle indicating the style with the argument endCapStyle =" SQUARE ".

# get the coordinates of Valencia
ciudad_point <- tmaptools::geocode_OSM("Valencia", as.sf = TRUE)

#  project the points
ciudad_point <- st_transform(ciudad_point, st_crs(buildings_val))

# create the buffer
point_bf <- st_buffer(ciudad_point, 2500) # radius of 2500 m


# get the intersection between the buffer and the building
buildings_val25 <- st_intersection(buildings_val, point_bf)
Warning: attribute variables are assumed to be spatially constant throughout
all geometries

Prepare data for mapping

We categorize the year into 15 groups using quartiles. It is also possible to modify the number of classes or the applied method (eg jenks, fisher, etc), you can find more details in the help ?classIntervals.

# find 15 classes
br <- classIntervals(year(buildings_val25$beginning), 15, "quantile")
Warning in classIntervals(year(buildings_val25$beginning), 15, "quantile"): var
has missing values, omitted in finding classes
# create labels
lab <- names(print(br, under = "<", over = ">", cutlabels = FALSE))
style: quantile
     < 1890 1890 - 1913 1913 - 1926 1926 - 1930 1930 - 1940 1940 - 1950 
        929        1368        1145         356        1687        1041 
1950 - 1958 1958 - 1962 1962 - 1966 1966 - 1970 1970 - 1973 1973 - 1978 
       1447        1025        1220        1156        1155        1190 
1978 - 1989 1989 - 2000      > 2000 
       1234        1132        1216 
# categorize the year
buildings_val25 <- mutate(buildings_val25,
  yr_cl = cut(year(beginning),
    br$brks,
    labels = lab,
    include.lowest = TRUE
  )
)

Map of Valencia

For the mapping, we will use the tmap package. It is an interesting alternative to ggplot2. It is a package of functions specialized in creating thematic maps. The philosophy of the package follows the same as in ggplot2, creating multiple layers with different functions, which always start with tm_ *and combine with +. Building a map with tmap always starts with tm_shape(), where the data, we want to draw, is defined. Then we add the corresponding geometry to the data type (tm_polygon(), tm_border(), tm_dots() or even tm_raster()). The tm_layout() function help us to configure the map style.

When we need more colors than the maximum allowed by RColorBrewer, we can pass the colors to the colorRampPalette() function. This function interpolates a set of given colors.

# colours
col_spec <- RColorBrewer::brewer.pal(11, "Spectral")

# colour ramp function
col_spec_fun <- colorRampPalette(col_spec)

# create the final map
tm_shape(buildings_val25) +
  tm_polygons("yr_cl",
    border.col = "transparent",
    palette = col_spec_fun(15), # adapt to the number of classes
    textNA = "Without data",
    title = ""
  ) +
  tm_layout(
    bg.color = "black",
    outer.bg.color = "black",
    legend.outside = TRUE,
    legend.text.color = "white",
    legend.text.fontfamily = "Montserrat",
    panel.label.fontfamily = "Montserrat",
    panel.label.color = "white",
    panel.label.bg.color = "black",
    panel.label.size = 5,
    panel.label.fontface = "bold"
  )
Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
font family not found in Windows font database

We can export our map using the function tmap_save("name.png", dpi = 300). I recommend using the dpi = 300 argument for a good image quality.

An alternative way to the tmap package is ggplot2.

# create the final map
ggplot(buildings_val25) +
  geom_sf(aes(fill = yr_cl), colour = NA) +
  scale_fill_manual(values = col_spec_fun(15), na.translate = FALSE) + # adapt to the number of classes
  labs(title = "VALÈNCIA", fill = NULL) +
  theme_void(base_family = "Montserrat") +
  theme(
    panel.background = element_rect(fill = "black"),
    plot.background = element_rect(fill = "black"),
    legend.justification = .5,
    legend.text = element_text(colour = "white", size = 12),
    legend.key.height = unit(3, "lines"),
    legend.key.width = unit(.5, "lines"),
    plot.title = element_text(
      colour = "white", hjust = .5, size = 60,
      margin = margin(t = 30)
    ),
    plot.caption = element_text(
      colour = "white",
      margin = margin(b = 20), hjust = .5, size = 16
    ),
    plot.margin = margin(r = 40, l = 40)
  )
Warning in grid.Call(C_stringMetric, as.graphicsAnnot(x$label)): font family
not found in Windows font database
Warning in grid.Call(C_stringMetric, as.graphicsAnnot(x$label)): font family
not found in Windows font database
Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
font family not found in Windows font database
Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
font family not found in Windows font database

To export the result of ggplot we can use the function ggsave("name.png").

Dynamic map with leaflet

A very interesting advantage is the tmap_leaflet() function of the tmap package to easily pass a map created in the same frame to leaflet.

# tmap object
m <- tm_shape(buildings_val25) +
  tm_polygons("yr_cl",
    border.col = "transparent",
    palette = col_spec_fun(15), # adapt to the number of classes
    textNA = "Without data",
    title = ""
  )


# dynamic map
tmap_leaflet(m)
Back to top

Reuse

Citation

For attribution, please cite this work as:
Royé, Dominic. 2019. “Visualize Urban Growth.” November 1, 2019. https://dominicroye.github.io/blog/visualize-urban-growth/.
Buy Me A Coffee