GIS och Kartor i R
Jun 10, 2020
Filip Wästberg
5 minute read

Geografiska informationssystem eller GIS är ett stort område inom dataanalys som traditionellt endast varit möjligt i mjukvaruprogram med höga licenskostnader. Men med den senaste tidens explosion av mjukvara som är Open Source är det idag möjligt att göra minst lika avancerade analyser helt kostnadsfritt. Här tänkte jag visa hur vi kan använda paket i R för att göra GIS.

Vi utgår från kartdata från SCB, som ni hittar här. Ladda ner Shape-filen som är en zippad fil, unzippa den och välj sedan geografiska indelning du är ute efter.

För att kunna läsa in kartor i R använder vi oss av RStudio, skapar ett nytt projekt (File -> New Project…) och lägger in hela mappen med filer för den geografiska indelning du valt.

Du kan så klart göra det här programmatiskt också genom:

download.file("https://www.scb.se/contentassets/3443fea3fa6640f7a57ea15d9a372d33/shape_svenska.zip", destfile = "shape_svenska.zip")

unzip("shape_svenska.zip", exdir = "shape_svenska")

unzip("shape_svenska/LanRT90.zip", exdir = "lan_rt90")

Då får vi dessa filer i mappen lan_rt90.

list.files("lan_rt90")
## [1] "Lan_RT90_region.dbf" "Lan_RT90_region.prj" "Lan_RT90_region.shp"
## [4] "Lan_RT90_region.shx"

Viktigt att notera är att du behöver alla filer i mappen för att ladda in geografisk data till R. Men du behöver bara ladda in en av filerna, nämligen den som är Shape (.shp).

För att läsa in och bearbeta kartfiler använder vi paktet sf(Simple Features) som gör det enkelt att läsa in och manipulera kartfiler.

Vi laddar in shape-filen med read_sf() och ser då att varje län har en geometry, en multipolygon, som är gränserna för länet.

library(sf)
## Linking to GEOS 3.8.1, GDAL 3.1.1, PROJ 6.3.1
lan <- read_sf("lan_rt90/Lan_RT90_region.shp")

lan
## Simple feature collection with 21 features and 2 fields
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: 1230810 ymin: 6137234 xmax: 1880916 ymax: 7669583
## projected CRS:  RT90 2.5 gon V
## # A tibble: 21 x 3
##    LnKod LnNamn                                                         geometry
##    <chr> <chr>                                                <MULTIPOLYGON [m]>
##  1 01    Stockholms  (((1581896 6569894, 1584367 6572819, 1585410 6571048, 1587…
##  2 03    Uppsala     (((1588105 6611272, 1587491 6610219, 1591541 6601881, 1591…
##  3 04    Södermanla… (((1508793 6538144, 1503218 6541877, 1502168 6544820, 1499…
##  4 05    Östergötla… (((1448643 6445587, 1446285 6441023, 1438051 6439351, 1437…
##  5 06    Jönköpings  (((1451287 6423652, 1452508 6419347, 1456711 6414227, 1455…
##  6 07    Kronobergs  (((1455396 6340554, 1456478 6342310, 1460661 6342030, 1458…
##  7 08    Kalmar      (((1495818 6327165, 1491423 6329015, 1489201 6333906, 1483…
##  8 09    Gotlands    (((1658877 6415526, 1662237 6416751, 1664956 6417647, 1666…
##  9 10    Blekinge    (((1423434 6230844, 1422094 6233619, 1418608 6233013, 1415…
## 10 12    Skåne       (((1370410 6257557, 1375642 6259726, 1380641 6259554, 1384…
## # … with 11 more rows

För att visualisera det här krävs inget mer än att använda ggplot2 tillsammans med geom_sf() på data.

library(tidyverse)
ggplot(lan) +
  geom_sf()

Med geom_sf() får vi ett koordinatsystem med longitude och latitude, men vill vi bara ha en “clean” karta kan vi exempelvis använda theme_map() från paketet ggthemes.

library(ggthemes)
ggplot(lan) +
  geom_sf() +
  theme_map()

I längden vill vi dock fylla den här kartan med någon slags information. Vi vänder oss då till SCB med paketet pxweb, som är ett paket för att hämta data från SCB direkt till R. Koden nedan är i huvudsak genererad från funktionen pxweb_interactive(), som vi använt för att interaktivt prata med SCB:s API. Nedan har jag tagit ner data för antalet personer eftergymnasial utbildning per län samt befolkning, för att kunna visualisera andelen i regionen med eftergymnasial utbildning.

Är du intresserad av hur man hämtar data från SCB till R har jag skrivit ett inlägg om det här. Är du inte intresserad av hur vi hämtar data kan du bara hoppa över det här kodblocket.

library(pxweb)
library(janitor)

pxweb_query_list <- list(
  "Region"= lan$LnKod, ## Här använder jag länskoder från shape-filen för att få ner alla län
  "UtbildningsNiva"=c("6"),
  "Alder"=c("16","17","18","19","20","21","22","23","24","25","26","27","28","29","30","31","32","33","34","35","36","37","38","39","40","41","42","43","44","45","46","47","48","49","50","51","52","53","54","55","56","57","58","59","60","61","62","63","64","65","66","67","68","69","70","71","72","73","74","75","76","77","78","79","80","81","82","83","84","85","86","87","88","89","90","91","92","93","94","95+"),
  "Kon"=c("1","2"),
  "ContentsCode"=c("000000I2"),
  "Tid"=c("2019"))

# Download data 
px_utb <- pxweb_get(url = "http://api.scb.se/OV0104/v1/doris/sv/ssd/UF/UF0506/UtbBefRegionR",
            query = pxweb_query_list)

# Convert to data.frame 
utb <- as.data.frame(px_utb, column.name.type = "text", variable.value.type = "text")

## Vi är bara ute efter att visualisera totalen per region, inte på åldersnivå, därför summerar jag per region
tbl_utb <- utb %>% 
  janitor::clean_names() %>% 
  group_by(region) %>% 
  summarise(eftergymn = sum(befolkning_16_95_ar))

pxweb_query_list <-  list(
  "Region"= lan$LnKod,
  "Alder"=c("tot"),
  "Kon"=c("1","2"),
  "ContentsCode"=c("BE0101A9"),
  "Tid"=c("2019"))

# Download data 
px_bef <- pxweb_get(url = "http://api.scb.se/OV0104/v1/doris/sv/ssd/BE/BE0101/BE0101A/FolkmangdNov",
            query = pxweb_query_list)

# Convert to data.frame 
bef <- as.data.frame(px_bef, column.name.type = "text", variable.value.type = "text")

tbl_bef <- bef %>% 
  clean_names() %>% 
  group_by(region) %>% 
  summarise(befolkning = sum(folkmangden_den_1_november))

tbl_eftergymn <- tbl_bef %>% 
  left_join(tbl_utb) %>% 
  mutate(procent_eftergymn = round(eftergymn /befolkning, 3))

Nu har vi data som vi kan joina på vår geografiska data. Ett problem är dock att vi i data från SCB:s API inte fått med länskoder. Det löser jag genom att rensa bort strängen " län" med str_replace(). Formatetet på länsnamn är då detsamma och vi kan joina på den geografiska information och visualisera andel personer med eftergymnasial utbildning per län.

tbl_eftergymn <- tbl_eftergymn %>% 
  mutate(region = str_replace(region, " län", "")) %>% 
  left_join(lan, by = c("region" = "LnNamn"))

ggplot(tbl_eftergymn, aes(fill = procent_eftergymn)) +
  geom_sf(aes(geometry = geometry)) +
  theme_map() +
  theme(legend.position = c(1, 0.15)) +
  scale_fill_viridis_c(option = "magma",
                       labels = scales::percent_format(accuracy = 1),
                       direction = -1) +
  labs(
    title = "Andel eftergymnasial utbildning per län",
    caption = "Källa: SCB",
    fill = "Eftergymnasial\nutbildning"
  )

Visualiseringen kan vi sedan spara som en png-fil.

ggsave("scb-karta-eftergymn.png")
## Saving 7 x 5 in image