[8]{.chapter-number}  [Centroids of transport zones]{.chapter-title}

8 Centroids of transport zones

In this section we will calculate the geometric and the weighted centroids of transport zones.

8.1 Geometric centroids

Taking the Municipalities_geo data from the previous section, we will calculate the geometric centroids, using the st_centroid() function.

Code

library(dplyr)
library(sf)
library(mapview)

Municipalities_geo = st_read("data/Municipalities_geo.gpkg", quiet = TRUE)

Centroids_geo = st_centroid(Municipalities_geo)

This creates points at the geometric center of each polygon.

Code

mapview(Centroids_geo)

Code

mapview(Centroids_geo) + mapview(Municipalities_geo, alpha.regions = 0) # both maps, with full transparency in polygons

But… is this the best way to represent the center of a transport zone?

These results may be biased by the shape of the polygons, and not represent where activities are. Example: lakes, forests, etc.

To overcome this, we can use weighted centroids.

8.2 Weighted centroids

We will weight centroids of transport zones by population and by number of buildings.

For this, we will need the census data (INE 2022).

Code

Census = st_read("data/census.gpkg", quiet = TRUE)

mapview(Census |> filter(Municipality == "Lisboa"), zcol = "Population")

It was not that easy to estimate weighted centroids with R, as it is with GIS software. But there is this new package centr that can help us (Zomorrodi 2024).

We need to specify the group we want to calculate the mean centroids, and the weight variable we want to use.

Code

# test
library(centr)
Centroid_pop = Census |> 
  mean_center(group = "Municipality", weight = "Population")

We can do the same for buildings.

Code

Centroid_build = Census |> 
  mean_center(group = "Municipality", weight = "Buildings")

8.3 Compare centroids in a map

8.3.1 Interactive map

Code

mapview(Centroids_geo, col.region = "blue") +
  mapview(Centroid_pop, col.region = "red") +
  mapview(Centroid_build, col.region = "black") +
  mapview(Municipalities_geo, alpha.regions = 0) # polygon limits

See how the building, population and geometric centroids of Sintra are apart, from closer to Tagus, to the rural area.

8.3.2 Static map

To produce the same map, using only plot() and st_geometry(), we need to make sure that the geometries have the same crs.

Code

st_crs(Centroids_geo) # 4326 WGS84
st_crs(Centroid_pop) # 3763 Portugal TM06

So, we need to transform the geometries to the same crs.

Code

Centroid_pop = st_transform(Centroid_pop, crs = 4326)
Centroid_build = st_transform(Centroid_build, crs = 4326)

Now, to use plot() we incrementally add layers to the plot.

Code

# Plot the Municipalities_geo polygons first (with no fill)
plot(st_geometry(Municipalities_geo), col = NA, border = "black")

# Add the Centroids_geo points in blue
plot(st_geometry(Centroids_geo), col = "blue", pch = 16, add = TRUE) # add!

# Add the Centroid_pop points in red
plot(st_geometry(Centroid_pop), col = "red", pch = 16, add = TRUE)

# Add the Centroid_build points in black
plot(st_geometry(Centroid_build), col = "black", pch = 16, add = TRUE)

Static map of different centroids of Municipalities

In the next section we will use these centroids to calculate the desire lines between them.

--- format: pdf: prefer-html: true --- # Centroids of transport zones In this section we will calculate the geometric and the weighted centroids of transport zones. ## Geometric centroids Taking the `Municipalities_geo` data from the previous section, we will calculate the geometric centroids, using the `st_centroid()` function. ```{r geocentroid} #| message: false #| warning: false library(dplyr) library(sf) library(mapview) Municipalities_geo = st_read("data/Municipalities_geo.gpkg", quiet = TRUE) Centroids_geo = st_centroid(Municipalities_geo) ``` This creates points at the geometric center of each polygon. ```{r} #| message: false #| fig-format: png mapview(Centroids_geo) mapview(Centroids_geo) + mapview(Municipalities_geo, alpha.regions = 0) # both maps, with full transparency in polygons ``` But... is this the best way to represent the center of a transport zone? These results may be biased by the shape of the polygons, and not represent where activities are. Example: lakes, forests, etc. To overcome this, we can use **weighted centroids**. ## Weighted centroids We will weight centroids of transport zones by population and by number of buildings. For this, we will need the census data [@INEcensus]. ```{r} #| eval: false #| include: false # data prep CENSUS = st_read("D:/GIS/MQAT/original/BGRI21_LISBOA.gpkg") CENSUS = CENSUS |> select(DTMN21, SUBSECCAO, N_INDIVIDUOS, N_EDIFICIOS_CLASSICOS) |> rename(Mun_code = "DTMN21", UID = "SUBSECCAO", Population = "N_INDIVIDUOS", Buildings = "N_EDIFICIOS_CLASSICOS") CENSUS = CENSUS |> left_join(Municipalities |> select(Mun_code, Municipality) |> distinct()) |> select(Municipality, UID, Population, Buildings) st_write(CENSUS, "data/census.gpkg") ``` ```{r getcensus} #| message: false #| fig-format: png Census = st_read("data/census.gpkg", quiet = TRUE) mapview(Census |> filter(Municipality == "Lisboa"), zcol = "Population") ``` It was not that easy to estimate weighted centroids with R, as it is with GIS software. But there is this new package [`centr`](https://ryanzomorrodi.github.io/centr) that can help us [@centr]. We need to specify the **group** we want to calculate the **mean centroids**, and the **weight variable** we want to use. ```{r} # test library(centr) Centroid_pop = Census |> mean_center(group = "Municipality", weight = "Population") ``` We can do the same for buildings. ```{r} Centroid_build = Census |> mean_center(group = "Municipality", weight = "Buildings") ``` ## Compare centroids in a map ### Interactive map ```{r} #| fig-format: png mapview(Centroids_geo, col.region = "blue") + mapview(Centroid_pop, col.region = "red") + mapview(Centroid_build, col.region = "black") + mapview(Municipalities_geo, alpha.regions = 0) # polygon limits ``` See how the building, population and geometric centroids of Sintra are apart, from closer to Tagus, to the rural area. ### Static map To produce the same map, using only `plot()` and `st_geometry()`, we need to make sure that the geometries have the same crs. ```{r} #| eval: false st_crs(Centroids_geo) # 4326 WGS84 st_crs(Centroid_pop) # 3763 Portugal TM06 ``` So, we need to transform the geometries to the same crs. ```{r} Centroid_pop = st_transform(Centroid_pop, crs = 4326) Centroid_build = st_transform(Centroid_build, crs = 4326) ``` Now, to use `plot()` we incrementally add layers to the plot. ```{r} #| fig.dpi: 300 #| fig-cap: "Static map of different centroids of Municipalities" #| fig-format: png # Plot the Municipalities_geo polygons first (with no fill) plot(st_geometry(Municipalities_geo), col = NA, border = "black") # Add the Centroids_geo points in blue plot(st_geometry(Centroids_geo), col = "blue", pch = 16, add = TRUE) # add! # Add the Centroid_pop points in red plot(st_geometry(Centroid_pop), col = "red", pch = 16, add = TRUE) # Add the Centroid_build points in black plot(st_geometry(Centroid_build), col = "black", pch = 16, add = TRUE) ``` In the [next section](desire-lines.qmd) we will use these centroids to calculate the desire lines between them.