download.file("mml:n latauspalvelun antama osoite",
destfile = "./local_file.zip")
unzip(zipfile = "./local_file.zip",
exdir = "./")
rm("local_file.zip")
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
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
<- "24.425354,60.275131,24.675293,60.360523"
bbox_string # Pilkotaan kukin luku omaksi konponentiksi numeeriseen vektoriin
<- strsplit(bbox_string, split = ",") %>% unlist() %>% as.numeric()
bbox_vector # Annetaan vektorille nimet
names(bbox_vector) <- c("xmin", "ymin", "xmax", "ymax")
# tehdään sf-objekti
<- st_as_sfc(st_bbox(bbox_vector, crs = st_crs(4326))) %>%
bbbox # Käännetään oikeaa CRS:ään
st_transform(3067)
# tehdään lopuksi WKT formaattiin
<- st_as_text(st_geometry(bbbox))
bboxWKT bboxWKT
[1] "POLYGON ((357610.6 6684831, 371429.3 6684318, 371764.6 6693824, 357982 6694337, 357610.6 6684831))"
# Käytetään filterointiin
<- sf::st_read("/mnt/laatokka/local_data/maastotietokanta/mtkmaasto.gpkg",
nuuksio 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
<- st_zm(nuuksio) 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
<- st_transform(nuuksio, 4326)
nuuksio_wgs84 leaflet(nuuksio_wgs84) %>%
addProviderTiles(provider = providers$OpenTopoMap) %>%
addPolygons(color = "orange")