MML:n maastotietokannan geopackage-tiedoston lukeminen R:ään suodattaen SQL:n ja spatiaalisen operaatioiden avulla

geopackage
SQL
spatial
data
R
maastotietokanta
mml
Tekijä

Markus Kainu

Julkaistu

26. tammikuuta 2023

Maanmittauslaitoksen maastotietokantaa voi käyttää eri tavoin. Tuoreessa blogissa Tilastokeskuksen ja Maanmittauslaitoksen OCP API Features rajapintojen hyödyntäminen R-kielellä esittelin maastotietokannan käyttöä rajapinnan kautta.

geopackage-tietokantatiedosto on yksi vaihtoehto, joka esitellään Maanmittauslaitokset verkkosivuilla: Maastotietokannan GeoPackage-versio.

Maastotietokannan GeoPackage-tietokantatiedostot ovat ladattavissa Avoimien aineistojen tiedostopalvelusta (Valitse tuote: Maastotietokanta > GeoPackage korkeus tai GeoPackage maasto). Asioinnin päätteeksi saat sähköpostiin latauslinkin, josta saat zip-pakatun tiedoston.

Tässä esimerkissä käytetään maasto versiota, joka lataamisen ja purkamisen päätteksi on levyllä nimellä mtkmaasto.gpkg. Tiedoston koko on tosiaan 71 gigatavua ja siksi se on liian iso ladattavaksi kokonaan tietokoneen muistiin. Tässä postauksessa esitellään miten pilkkoa tiedosto pienemmäksi SQL-hakulauseiden tai spatiaalisen rajausten avulla.

R:llä tiedosto ladataan ja puretaan seuraavasti

download.file("mml:n latauspalvelun antama osoite", 
              destfile = "./local_file.zip")
unzip(zipfile = "./local_file.zip", 
      exdir = "./")
rm("local_file.zip")

Sitten listataan geopackage-tiedoston layerit, joita tässä on 122.

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


st_layers("./mtkmaasto.gpkg")

Tässä olen erityisesti kiinnostunut järvistä, mutta kuten tiedetään niitä Suomessa riittää. Siksi haluan aluksi rajata haun pelkästään Nuuksion alueen järviin. Tehdään siis http://bboxfinder.com-palvelussa neliskanttinen laatikko Nuuksion päälle (ks. kuva) ja valitaan siitä laatikon ääripisteiden koordinaatit: 24.425354,60.275131,24.675293,60.360523.

Tätä kahden koordinaattipisteen tieto tulee muokata ETRS89 / TM35FIN(E,N) koordinaattijärjestelään sekä käsitellä Well-known text representation of geometry -muotoon eli WKT:ksi.

# Otetaan ensin string
bbox_string <- "24.425354,60.275131,24.675293,60.360523"
# Pilkotaan kukin luku omaksi konponentiksi numeeriseen vektoriin
bbox_vector <- strsplit(bbox_string, split = ",") %>% unlist() %>% as.numeric()
# Annetaan vektorille nimet
names(bbox_vector) <- c("xmin", "ymin", "xmax", "ymax")  
# tehdään sf-objekti
bbbox <- st_as_sfc(st_bbox(bbox_vector,  crs = st_crs(4326))) %>%
  # Käännetään oikeaa CRS:ään 
  st_transform(3067)
# tehdään lopuksi WKT formaattiin
bboxWKT <- st_as_text(st_geometry(bbbox))
bboxWKT
[1] "POLYGON ((357610.6 6684831, 371429.3 6684318, 371764.6 6693824, 357982 6694337, 357610.6 6684831))"
# Käytetään filterointiin
nuuksio <- sf::st_read("/mnt/laatokka/local_data/maastotietokanta/mtkmaasto.gpkg", 
                      layer = "jarvi", 
                      wkt_filter = bboxWKT
                      )
Reading layer `jarvi' from data source 
  `/mnt/laatokka/local_data/maastotietokanta/mtkmaasto.gpkg' 
  using driver `GPKG'
Simple feature collection with 164 features and 10 fields
Geometry type: POLYGON
Dimension:     XYZ
Bounding box:  xmin: 357636.3 ymin: 6681899 xmax: 371491.1 ymax: 6695818
z_range:       zmin: -999.999 zmax: 97.3
Projected CRS: ETRS89 / TM35FIN(E,N)
# poistetaan Z-ulottuvuus
nuuksio <- st_zm(nuuksio)

Tästä nähdään että bboxiin osui yhteensä 164 järveä. Samalla voidaan katsoa mitä muuttujia jarvi-kerros pitää sisällään. Erityisesti kiinnostaa muuttuja keskikorkeus SQL-filtterointiä varten.

Mutta piirretään ensin vielä vuorovaikutteilen leaflet-kartta järvistä niin että järven väri on oranssi.

library(leaflet)
# Muokataan takaisin WGS84 koska Leaflet
nuuksio_wgs84 <- st_transform(nuuksio, 4326)
leaflet(nuuksio_wgs84) %>% 
  addProviderTiles(provider = providers$OpenTopoMap) %>% 
  addPolygons(color = "orange")