Created
January 2, 2024 19:06
-
-
Save sumn2u/32dcc73f8dcad5a60322cf51637ffe45 to your computer and use it in GitHub Desktop.
Slovenia Forest Map
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| ############################################# | |
| Slovenia | |
| ############################################# | |
| # 1. PACKAGES | |
| # require(remotes) | |
| # install_version("rayshader", version = "0.35.7", repos = "http://cran.us.r-project.org") | |
| libs <- c( | |
| "giscoR", "terra", "elevatr", | |
| "png", "rayshader", "magick" | |
| ) | |
| installed_libs <- libs %in% rownames( | |
| installed.packages() | |
| ) | |
| if(any(installed_libs == F)){ | |
| install.packages( | |
| libs[!installed_libs] | |
| ) | |
| } | |
| invisible(lapply( | |
| libs, library, | |
| character.only = T | |
| )) | |
| # 1. COUNTRY BOUNDARIES | |
| #---------------------- | |
| country_sf <- giscoR::gisco_get_countries( | |
| country = "SI", | |
| resolution = "1" | |
| ) | |
| # 2. FOREST TYPE RASTER | |
| #---------------------- | |
| url <- "https://s3-eu-west-1.amazonaws.com/vito.landcover.global/v3.0.1/2019/E000N60/E000N60_PROBAV_LC100_global_v3.0.1_2019-nrt_Forest-Type-layer_EPSG-4326.tif" | |
| download.file( | |
| url = url, | |
| destfile = basename(url), | |
| mode = "wb" | |
| ) | |
| forest_type <- terra::rast( | |
| basename(url) | |
| ) | |
| vals <- terra::values( | |
| forest_type, | |
| dataframe = T | |
| ) | |
| names(vals) | |
| names(vals)[1] <- "value" | |
| unique(vals$value) | |
| # 3. CROP FOREST TYPE RASTER | |
| #--------------------------- | |
| crs_lambert <- | |
| "+proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 +datum=WGS84 +units=m +no_frfs" | |
| country_forest_type <- terra::crop( | |
| forest_type, | |
| terra::vect(country_sf), | |
| snap = "in", | |
| mask = T | |
| ) |> | |
| terra::project(crs_lambert) | |
| terra::plot(country_forest_type) | |
| # 4. FOREST TYPE RASTER TO IMAGE | |
| #------------------------------- | |
| cols <- c( | |
| "#073b4c", | |
| "#7cd8ff", | |
| "#ffcf6a", | |
| "#e73f69" | |
| ) | |
| from <- c(0:1, 4:5) | |
| to <- t(col2rgb( | |
| cols | |
| )) | |
| forest_terra <- na.omit( | |
| country_forest_type | |
| ) | |
| forest_type_image <- terra::subst( | |
| forest_terra, | |
| from, | |
| to, | |
| names = cols | |
| ) | |
| terra::plotRGB(forest_type_image) | |
| img_file <- "slovenia-forest-image.png" | |
| terra::writeRaster( | |
| forest_type_image, | |
| img_file, | |
| overwrite = T, | |
| NAflag = 255 | |
| ) | |
| img <- png::readPNG(img_file) | |
| img | |
| # 5. COUNTRY ELEVATION RASTER | |
| #---------------------------- | |
| elev <- elevatr::get_elev_raster( | |
| locations = country_sf, #warn about locations | |
| z = 10, clip = "locations" | |
| ) | |
| elev_lambert <- elev |> | |
| terra::rast() |> | |
| terra::project(crs_lambert) | |
| elmat <- rayshader::raster_to_matrix( | |
| elev_lambert | |
| ) | |
| # 6. RENDER SCENE | |
| #---------------- | |
| h <- nrow(elev_lambert) | |
| w <- ncol(elev_lambert) | |
| elmat |> | |
| rayshader::height_shade( | |
| texture = colorRampPalette( | |
| "white" | |
| )(512) | |
| ) |> | |
| rayshader::add_overlay( | |
| img, | |
| alphalayer = .9, | |
| alphacolor = "white" | |
| ) |> | |
| rayshader::add_shadow( | |
| rayshader::lamb_shade( | |
| elmat, | |
| zscale = 50, | |
| sunaltitude = 90, | |
| sunangle = 315, | |
| ), max_darken = .25 | |
| ) |> | |
| rayshader::add_shadow( | |
| rayshader::texture_shade( | |
| elmat, | |
| detail = .95, | |
| brightness = 90, #warn | |
| contrast = 80, | |
| ), max_darken = .1 | |
| ) |> | |
| rayshader::plot_3d( | |
| elmat, | |
| zscale = 5, | |
| solid = F, | |
| shadow = T, | |
| shadow_darkness = 1, | |
| background = "white", | |
| windowsize = c( | |
| w / 5, h / 5 | |
| ), | |
| zoom = .5, | |
| phi = 85, | |
| theta = 0 | |
| ) | |
| # 7. RENDER OBJECT | |
| #----------------- | |
| tempfilehdr = tempfile(fileext = ".hdr") | |
| tempfilehdr | |
| download.file("https://github.com/milos-agathon/3d-forest-type-map/raw/main/air_museum_playground_4k.hdr",tempfilehdr) | |
| rayshader::render_highquality( | |
| filename = "forest-type-slovenia-3d.png", | |
| preview = T, | |
| light = F, | |
| environment_light = tempfilehdr, | |
| intensity_env = 2, | |
| rotate_env = 90, | |
| interactive = F, | |
| parallel = T, | |
| width = w, height = h | |
| ) | |
| # 8. MAKE LEGEND | |
| #--------------- | |
| png("my_legend.png") | |
| par(family = "mono") | |
| plot( | |
| NULL, xaxt = "n", | |
| yaxt = "n", bty = "n", | |
| ylab = "", xlab = "", | |
| xlim = 0:1, ylim = 0:1, | |
| xaxs = "i", yaxs = "i" | |
| ) | |
| legend( | |
| "center", | |
| legend = c( | |
| "Unknown", | |
| "Evergreen needle leaf", | |
| "Deciduous broad leaf", | |
| "Mixed" | |
| ), | |
| pch = 16, | |
| pt.cex = 3, | |
| cex = 1.5, | |
| bty = "n", | |
| col = cols | |
| ) | |
| dev.off() | |
| # 9. FINAL MAP | |
| #------------- | |
| forest_img <- magick::image_read( | |
| "forest-type-slovenia-3d.png" | |
| ) | |
| my_legend <- magick::image_read( | |
| "my_legend.png" | |
| ) | |
| my_legend_scaled <- magick::image_scale( | |
| magick::image_background( | |
| my_legend, "none"), 2000 | |
| ) |> | |
| magick::image_transparent("white") | |
| p <- magick::image_composite( | |
| magick::image_scale( | |
| forest_img, | |
| "x4000" | |
| ), | |
| my_legend_scaled, | |
| offset = "+100+0" | |
| ) | |
| magick::image_write( | |
| p, | |
| "final-map.png" | |
| ) |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment