R-ohjelmointi.org

Tilastotieteellistä ohjelmointia R-kielellä

Laserkeilausaineistot ja korkeusmallit R:ssä rayshader-paketilla höystettynä

Laserkeilausaineisto kuvaa maanpintaa ja sen muotoja. Laserkeilausaineistoa voidaan käyttää esimerkiksi korkeusmallien muodostamiseen. Laserkeilausaineistoja on saatavilla toistaiseksi vain osasta maata avoimena datana Maanmittauslaitoksen karttapalvelusta. Helsingin alueen laserkeilausaineistoa voi ladata myös Helsingin karttapalvelusta.

Laserkeilausaineiston käyttö varjostusanalyysiin

Laserkeilausaineistoja voi ladata R:ään paketilla lidR. Aineiston perusteella voidaan sitten helposti tutkia esimerkiksi kaupungin varjoisuutta rayshader-paketin tarjoamilla työkaluilla. Katsotaanpa esimerkkiä Helsingin Kampista. Yllä olevasta Helsingin karttapalvelusta ladattiin tiedosto rgb_672496b.laz. Kyseinen karttalehti kattaa mm. Kampin metroaseman seudun.

Aineisto saadaan ladattua R:n objektiin las vaikkapa seuraavasti:

library(lidR)
library(reshape2)
library(rayshader)
 
setwd("C:/users/lenovo/desktop/lidar")
 
las <- readLAS("rgb_672496b_kamppi.laz")

Laserkeilausaineiston perusteella voidaan luoda esimerkiksi interaktiivinen kuva rgl-kirjastoa käyttäen:

plot(las)

Kuvassa ylhäällä ja keskellä on Kampin kauppakeskus ja Tennispalatsin aukio.

Laserkeilausaineiston perusteella on mahdollista tutkia esimerkiksi varjostusta. Seuraava koodi laskee rakennusten aiheuttaman varjostuksen päivän mittaan säteenseurantaa:

hmean = grid_metrics(las, mean(Z), res=1.1)
hmean2 <-  as.matrix(dcast(hmean, X~Y))
amb = ambient_shade(heightmap = as.matrix(hmean2), sunbreaks = 60, maxsearch = 100)
plot_map(amb)

Muodostuvassa kuvassa varjon määrää kuvataan valkoinen-musta -asteikolla, jossa varjoisimmat alueet ovat kaikkein tummimpia:

Kuva on nyt eri asennossa kuin edellinen, ja Kampin kauppakeskus on kuvassa alhaalla ja keskellä.

Korkeusmalli ja meren pinta

Korkeusmalli on helpoin tuoda R:ään elevatr-pakettiä käyttäen. Se ei käytä edellä luettua laserkeilausaineistoa, vaan hakee malli julkisesta rajapinnasta. Korkeusmallin avulla voidaan sitten tutkia vaikkapa, miten merenpinnan nousu vaikuttaisi Helsingin keskustaan. Aineiston hakemiseen tarvitaan koordinaatit ja niitä vastaava projektiotieto:

library(elevatr)
library(rayshader)
library(sp)
 
loc_df <- data.frame(x=24.95, y=60.18)
 
prj_dd <- "+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs"
 
x <- get_elev_raster(locations = loc_df, prj = prj_dd, z=12, src="aws", api_key=NULL)

Aineisto voidaan plotata seuraavasti:

plot(x)

Kuvassa oikealla alhaalla on Helsingin keskusta ja ylhäällä oikealle muun muassa Kallio. Google Mapsissä vastaava kohta on suunnilleen https://www.google.com/maps/@60.1712453,24.907984,14z.

Kyseiseen korkeusmalliin voidaan lisätä vesi esimerkiksi seuraavasti:

elevation <- x
elmat2 = matrix(raster::extract(elevation,raster::extent(elevation),buffer=10000),nrow=ncol(elevation),ncol=nrow(elevation))
 
s1 <- sphere_shade(elmat2)
s2 <- add_water(s1, detect_water(elmat2, cutoff=0.90, min_area=700))
 
plot_map(s2)

Tuloksena on seuraavanlainen kuva. Ei täydellinen, mutta vesialueet ovat karkealla tasolla sinne päin:

Käyttämällä 3D-plottia (rgl-pakettia käyttäen), saadaan kuvaan veden pinta määritellyllä tasolla (tässä mustalla värillä). Seuraavassa koodissa luodaan ensin kuva nykyisellä merenpinnan tasolla, ja sitten metriä korkeammalla tasolla:

plot_3d(s2, elmat2, water=TRUE, waterdepth=18, shadow=T, watercolor="black", zoom=0.75, wateralpha=1)
snapshot3d("C:/users/lenovo/desktop/rgl_18m.png")
 
plot_3d(s2, elmat2, water=TRUE, waterdepth=19, shadow=T, watercolor="black", zoom=0.75, wateralpha=1)
snapshot3d("C:/users/lenovo/desktop/rgl_19m.png")

Kuvat voidaan yhdistää vierekkäin magick-paketilla vertailun helpottamiseksi:

library(magick)
img1 <- image_read("C:/users/lenovo/desktop/rgl_18m.png")
img2 <- image_read("C:/users/lenovo/desktop/rgl_19m.png")
image_write(image_append(c(img1, img2)), "C:/users/lenovo/desktop/cmbn.png")

Tuloksena muodostuvasta kuvasta on nähtävissä, ettei merenpinnan nopeakaan kohoaminen metrillä näyttäisi muuttavat Helsingin rantaviivaa erityisen radikaalisti, vaikka sillä toki olisi vaikutuksia mm. Kauppatoriin (kuvassa oikealla alhaalla, noin 1/3 alareunasta ylöspäin):

Miksikö juuri metri? Se on IPCC:n arvio merenpinnan noususta vuoteen 2100 mennessä. Tosin, samassa ajassa maanpinta kohoaa noin 30 cm, joten vaikutus ei ehkä Suomessa ole niin suuri. Wikipedian artikkelia kokoaa yhteen arvioiden tiedot.

Edit 2018-10-15: Lisätietoja rayshader-paketistä löytyy muun muassa seuraavista Tyler Morgan-Wallin blogiartikkeleista:
http://www.tylermw.com/3d-maps-with-rayshader/
http://www.tylermw.com/making-beautiful-maps/
http://www.tylermw.com/throwing-shade/
http://www.tylermw.com/3d-printing-rayshader/


Category