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

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")

Kauniaisten järvet saa haettua seuraavasti:

kunta <- geofi::get_municipalities() %>% filter(municipality_name_fi == "Kauniainen")
kunta_wkt <- kunta %>% 
  st_geometry() %>% # convert to sfc
  st_as_text() # convert to well known text

# Käytetään filterointiin
kunta_jarvet <- sf::st_read("/mnt/laatokka/local_data/maastotietokanta/mtkmaasto.gpkg", 
                      layer = "jarvi", 
                      wkt_filter = kunta_wkt
                      )

Reading layer `jarvi' from data source 
  `/mnt/laatokka/local_data/maastotietokanta/mtkmaasto.gpkg' 
  using driver `GPKG'
Simple feature collection with 7 features and 10 fields
Geometry type: POLYGON
Dimension:     XYZ
Bounding box:  xmin: 373202.3 ymin: 6676719 xmax: 376942.9 ymax: 6683543
z_range:       zmin: -999.999 zmax: 43.9
Projected CRS: ETRS89 / TM35FIN(E,N)

# poistetaan Z-ulottuvuus
kunta_jarvet <- st_zm(kunta_jarvet)
mapview(kunta_jarvet)

Nyt kun saimme Nuuksion järvet kivasti kartalle, heräsi mielenkiinto piirtää kartalle Suomen korkeimmalla olevat järvet. jarvi-layerillä on järven keskikorkeus ja jos haluamme vain yli 700m ja alle 1300m korkeudessa olevat järvet niin voimme rajata ne SQL-lausekkeella seuraavasti.

korkeat <- sf::st_read("/mnt/laatokka/local_data/maastotietokanta/mtkmaasto.gpkg", 
                      layer = "jarvi", 
                      query= 'SELECT * FROM jarvi WHERE keskikorkeus > 700000 AND keskikorkeus < 1300000'
                      )

Reading query `SELECT * FROM jarvi WHERE keskikorkeus > 700000 AND keskikorkeus < 1300000'
from data source 
  `/mnt/laatokka/local_data/maastotietokanta/mtkmaasto.gpkg' 
  using driver `GPKG'
Simple feature collection with 1020 features and 10 fields
Geometry type: POLYGON
Dimension:     XYZ
Bounding box:  xmin: 245499.5 ymin: 7630551 xmax: 296196 ymax: 7699607
z_range:       zmin: -999.999 zmax: 1191
Projected CRS: ETRS89 / TM35FIN(E,N)

# poistetaan Z-ulottuvuus
korkeat_wgs84 <- st_zm(korkeat) %>%
  st_transform(4326)
leaflet(korkeat_wgs84) %>% 
  addProviderTiles(provider = providers$OpenTopoMap) %>% 
  addPolygons(color = "orange")