Intro

Lentolaserkeilauksessa tarkasteltavan kohteen kolmiulotteista rakennetta kuvaava pistepilvi tuotetaan ilma-aluksesta kohteen yläpuolelta. Laserkeilaus perustuu suunnattuun laseretäisyysmittaukseen: Laserkeilaimelta lähetetty laserpulssi osuu kohteeseen, ja osa siitä heijastuu paluukaikuna takaisin laserkeilaimelle. Lähetetyn pulssin ja rekisteröidyn paluukaiun välisestä aikaerosta voidaan laskea pulssin kulkema matka. Kun keilaimen tarkka sijainti sekä lähetetyn laserpulssin suunta tiedetään, voidaan kohteen kolmiulotteinen sijainti määrittää. Metsien kaukokartoituksessa hyödynnettävä lentolaserkeilausaineisto kerätään lentolinja kerrallaan. Lentolinjan leveyteen vaikuttavat keilauksen avauskulma sekä lentokorkeus. Mitä suurempi on keilauksen avauskulma, sitä suurempi alue saadaan keilattua samalta lentokorkeudelta. Lentokorkeutta nostamalla keilauslinjan leveys maastossa kasvaa, mutta samalla tuotettavan pistepilven tiheys laskee. Lisätietoa tässä projektissa käytetyistä ALS-aineistoista löytyy aineistokuvauksissa (http://www.helsinki.fi/~korpela/FOR254/html/aineisto.html#kp).

Lentolaserkeilauksessa jokaista lähetettyä laserpulssia vastaa yksi tai useampi paluukaiku, ja paluukaikujen lukumäärä riippuu tarkasteltavan kohteen ominaisuuksista. Jokainen paluukaiku puolestaan rekisteröidään pistepilveen erillisinä pisteinä. Kun laserpulssi osuu latvustoon, useimmiten osa pulssista heijastuu latvan yläosista ja osa latvuston sisältä ja alemmista latvuskerroksista. Latvusaukoissa laserpulssi heijastuu takaisin keilaimelle viimeistään maanpinnasta. ALS-aineiston laserpisteet voidaankin luokitella paluukaiun järjestysnumeron perusteella ensimmäisiin (first), keskimmäisiin (intermediate) ja viimeisiin (last echoes) kaikuihin. Ensimmäisiä kaikuja kuvaavat pisteet kuvaavat parhaiten latvuston pituutta ja tiheyttä, kun taas viimeiset kaiut kuvaavat useimmiten maanpinnan korkeutta. ALS-aineistossa pisteiden korkeudet on useimmiten määritelty suhteessa sovittuun nollatasoon, yleensä merenpinnantasoon. Viimeisten kaikujen perusteella voidaan muodostaa maanpinnan korkeutta merenpinnasta kuvaava maastomalli, jonka avulla pistepilven pisteiden z-koordinaatit voidaan normalisoida korkeuksiksi maanpinnasta. Korkeusnormalisoidusta ALS-pistepilven avulla saadaan tietoa yksittäisten puiden, puuryhmien ja metsiköiden rakenteesta. Pistepilvestä voidaan automaattisesti tunnistaa yksittäisiä puita, mitata puuston pituutta ja tiheyttä sekä ennustaa puuston rakennetta kuvaavia puustotunnuksia (esim. White ym. 2013, Kotivuori ym. 2016). Näitä varten ALS-aineisto jaetaan xy-tasossa hilaruutuihin, ja ruutuihin osuneista pisteistä lasketaan pisteiden tiheys- ja korkeusjakaumaa kuvaavia laserpiirteitä, joista muodostetaan edelleen rastereita.

ALS-aineistojen käsittelyssä ja analysoinnissa käytettiin LAStools-ohjelmistoa (http://rapidlasso.com/LAStools) sekä R-ohjelmistoa (https://www.r-project.org/) R-studio-ympäristössä (https://www.rstudio.com/). LAStools on laserkeilausaineistojen tehokkaaseen prosessointiin kehitetty ohjelmistopaketti, joka koostuu sekä vapaasti hyödynnettävistä että lisensoiduista työkaluista. R on puolestaan avoimeen lähdekoodiin perustuva ohjelmisto, joka monipuolisten työkalupakettien myötä soveltuu erinomaisesti opetus- ja tutkimuskäyttöön sekä data-analyyseihin. Erityisesti R-pakettien sp (https://cran.r-project.org/web/packages/sp/index.html)ja raster (https://cran.r-project.org/web/packages/raster/) ansiosta R:llä voidaan käsitellä tehokkaasti myös paikkatietoaineistoja, kuten ALS-piirrerastereita ja shapefile-tiedostoja.

Piirteiden laskenta

Tässä dokumentissa käydään läpi ALS-piirteiden laskennan vaiheet LAStools-työkaluja ja R:ää hyödyntämällä. R-skripti yhden ALS-vuosikerran piirteiden laskemiseen löytyy tämän projektin OneDrive-kansiosta (https://helsinkifi-my.sharepoint.com/:f:/g/personal/korpela_ad_helsinki_fi/Ei0CRRUz5NRGqoqkUnf3F7sB_XPF8eviGxE4OXGtCgtgsw)

Ensin otetaan käyttöön tarvittavat R-paketit, jotka tulee ensimmäistä käyttökertaa varten asentaa install.packages()-komennolla.

library(sp); library(raster); library(rgdal); library(rgl);library(maptools);library(gridExtra);library(rasterVis)

LAStools-työkaluja voidaan ajaa R:llä komentorivin kautta. Tätä varten määritellään tiedostopolku LAStools-työkaluihin ja muodostetaan funktio, jolla LAStools-työkalukomentoja ajetaan f(“LAStools-komento”):

ld <- "C:/LAStools/bin/" 

f <- function(las_string){
  command <- paste0(ld, las_string)
  print(command)
  system(command,wait=T)
}

Eri ajankohtina kerätty ALS-aineisto on tallennettuna 500 m x 500 m tiileihin, joten n. 2 km x 4 km kokoiselle alueelle ulottuvaa n. 300 ha mallimetsälöä varten tarvittiin yhteensä 32 erillistä las-tiedostoa jokaista ALS-vuosikertaa kohden. ALS-piirteiden laskentaa varten nämä yhden vuosikerran 500 m x 500 m las-palat ensin yhdistetään käyttämällä lasmerge-työkalua (http://www.cs.unc.edu/~isenburg/lastools/download/lasmerge_README.txt).

# Tässä oletuksena, että työskentelykansiossa on vain yhden ajankohdan .las-tiedostot, jotka yhdistetään
f(paste0("lasmerge -i *.las -o merge.las"))

Seuraavaksi yhdistetty pistepilvi rajataan tarkastelualueen rajojen sisään lasclip-työkalulla (http://www.cs.unc.edu/~isenburg/lastools/download/lasclip_README.txt) rajoja kuvaavan shapefilen avulla. Pisteet, jotka ovat shapefilen ulkopuolella, poistetaan. Näin saadaan kohdennettua analyysit tarkastelualueelle ja samalla hieman pienennettyä käsiteltävän datan määrää.

f(paste0("lasclip -i merge.las -poly vjsuo_rajat.shp -o clip.las"))

Kuva 1. Las-tiedostojen yhdistäminen ja rajaus mallimetsälön ALS-paikkatietoanalyysejä varten

Yhdistetystä ja rajatusta pistepilvestä lasketaan ALS-piirteet lascanopy-työkalulla (http://www.cs.unc.edu/~isenburg/lastools/download/lascanopy_README.txt). Piirteet lasketaan pulssien ensimmäisistä/ainoista kaiuista, ja keilauksen avauskulma rajoitetaan 15 asteeseen. Ensimmäiset kaiut ovat muita kaikuja stabiilimpia kuvaamaan latvuston pituutta, tiheyttä ja niiden vaihtelua, sillä ne heijastuvat takaisin keilaimelle heti pulssin osuessa latvuston pintaan. Latvuston tiheyttä kuvataan tässä projektissa latvuston läpäisevien (maahan osuvien, h < 2 m) kaikujen osuutena kaikista yli 2 m korkeudelta tulleista kaiuista. Suurilla keilauskulmilla pulssien läpäisemättömyys ei johdu ainoastaan latvuston tiheydestä vaan myös mittausgeometriasta: kun pulssi osuu latvustoon riittävän suuressa kulmassa, se ei läpäise latvustoa yhtä hyvin kuin suoraan ylhäältä latvustoon osuva pulssi, vaan heijastuu takaisin keilaimelle latvuston alemmista osista.

Oheisessa esimerkissä ALS-piirteet lasketaan 16 m hilaruutuihin, ja tiheys- ja korkeuspiirteiden laskennassa käytettävä korkeusraja on 2 m. Tyypillisiä laserpiirteitä ovat latvuston peittävyys (cov), pisteiden korkeusjakauman 95 %:n prosenttipistet sekä keskiarvo ja maksimiarvo. Piirrerasterit tallennetaan .tif-tiedostoina ETRS-TM35FIN-koordinaatistoon.

f(paste0("lascanopy -i clip.las -step 16 -first_only -keep_scan_angle -15 15 -cov -cover_cutoff 2.00005 -fractions -height_cutoff 2.00001 -max -avg -p 95 -etrs89 -otif -utm 35north"))

Muodostetaan lopuksi piirrekartat, jotka kuvaavat latvuston peittävyyttä (cov) ja pituutta (p95).

# Väripaletti kartan piirtoon
mycols <- rasterTheme(region=colorRampPalette(brewer.pal(9,'Greens'))(100))
# Latvuspeitto
p1<-levelplot(raster("clip_cov.tif")*100, margin = FALSE, main = "cov (%) 2016, 16 m x 16 m", par.settings = mycols, zlim = c(0,100)) + layer(sp.lines(vjsuo_rajat, lwd=0.8, col='black'))
# Pituus
p2<-levelplot(raster("clip_p95.tif"), margin = FALSE, main = "p95 (m) 2016, 16 m x 16 m", par.settings = mycols, zlim = c(0,35)) + layer(sp.lines(vjsuo_rajat, lwd=0.8, col='black'))
# Tulostetaan kartat vierekkäin
grid.arrange(p1, p2, ncol=2)

Kuva 2. Kartat latvuston tiheyttä ja pituutta kuvaavista laserpiirteistä mallimetsälön alueella.

Viitteet

Kotivuori, E., Korhonen, L. and Packalen, P., 2016. Nationwide airborne laser scanning based models for volume, biomass and dominant height in Finland. Silva Fennica, 50(4), p.28.

White, J.C., Wulder, M.A., Varhola, A., Vastaranta, M., Coops, N.C., Cook, B.D., Pitt, D. and Woods, M., 2013. A best practices guide for generating forest inventory attributes from airborne laser scanning data using an area-based approach. The Forestry Chronicle, 89(6), pp.722-723.