(Acest articol a fost publicat pentru prima dată pe r.iresmi.netși cu amabilitate a contribuit la R-bloggeri). (Puteți raporta problema legată de conținutul acestei pagini aici)
Doriți să vă distribuiți conținutul pe R-bloggeri? dați clic aici dacă aveți un blog, sau aici dacă nu aveți.
Ziua 1 din 30DayMapChallenge: « Puncte » (anterior).
IGN scanează Franța folosind LIDAR. Mai mult de jumătate din țară este disponibilă pentru descărcare în prezent. Să ne jucăm cu punctele din jurul Aiguille du Midi…
library(readr)
library(dplyr)
library(purrr)
library(fs)
library(lidR)
dir_create("tiles")
Selectarea a 4 pătrate pe harta de selecție permite generarea unui catalog de piese de 1✕1 km de descărcat.
read_lines("liste_dalle.txt") |>
tibble(url = _) |>
mutate(file = path("tiles", path_file(url))) |>
pwalk((url, file) download.file(url, file))
Citiți și manipulați fișierele LAZ cu ajutorul pachetului {lidR}. Păstrăm doar o regiune de interes la 500 m în jurul vârfului.
ctg <- readLAScatalog("tiles")
roi <- clip_circle(ctg, x = 1001400, y = 6538400, radius = 500)
Norul de puncte a fost deja clasificat așa că putem extrage doar pământul (deși la această altitudine nu există vegetație).
ground <- filter_ground(roi) plot(ground)
Afișăm 38 de milioane de puncte (50 de puncte/m²)!
Și apoi putem crea un model digital de teren și o hartă de deal…
dtm <- rasterize_terrain(roi, 2, tin(), pkg = "terra")
dtm_prod <- terra::terrain(dtm, v = c("slope", "aspect"), unit = "radians")
dtm_hillshade <- terra::shade(slope = dtm_prod$slope, aspect = dtm_prod$aspect)
plot(dtm_hillshade, col = gray(0:50/50), legend = FALSE)
