close terrain dem (#49)

* raster digest working in principle

* reliéfy se vcelku líbí

* tests complete locally

* vignette updated

* tests rewrite

* typo tuning

.gitignore

@@ -8,4 +8,5 @@ inst/doc
 /data-raw/vodstvo/*
 /data-raw/TRANS/*
 /data-raw/HYDRO/*
+/data-raw/eu_dem*
 *.log

DESCRIPTION

@@ -1,8 +1,8 @@
 Package: RCzechia
 Type: Package
 Title: Spatial Objects of the Czech Republic
-Version: 1.9.5
-Date: 2022-11-08
+Version: 1.10.0
+Date: 2022-11-25
 Authors@R: 
     person("Jindra", "Lacko", , "jindra.lacko@gmail.com", role = c("aut", "cre"),
     comment = c(ORCID = "0000-0002-0375-5156"))
@@ -21,7 +21,7 @@ Imports:
   jsonlite,
   magrittr
 Suggests:
-  raster,
+  terra,
   readxl,
   units,
   ggplot2,
@@ -37,11 +37,12 @@ Suggests:
   dplyr,
   lwgeom,
   s2,
-  covr
+  covr,
+  tidyterra
 SystemRequirements: 
   GDAL (>= 2.2.3), 
   GEOS (>= 3.6.2),
   PROJ (>= 4.9.3)
-RoxygenNote: 7.2.1
+RoxygenNote: 7.2.2
 VignetteBuilder: knitr
 NeedsCompilation: no

NEWS.md

@@ -1,3 +1,9 @@
+## version 1.10.0 (2022-11-25)
+
+ - rasters - vyskopis("actual") and vyskopis("rayshaded") - are based on EU-DEM data, instead of the former ARC ČR 500
+
+ - dependency on {raster} replaced by newer {terra} package
+
 ## version 1.9.4 (2022-10-06)
 
  - fixes incorrect behavior of Vltava in Prague

R/downloader.R

@@ -24,8 +24,13 @@
     curl::curl_download(url = remote_file, destfile = local_file, quiet = T)
    } # /if - local file exists
 
-  local_df <- readRDS(local_file)
+  # everything except rasters
+  if(tools::file_ext(local_file) == "rds") local_df <- readRDS(local_file)
 
+  # rasters, and rasters only
+  if(tools::file_ext(local_file) == "tif") local_df <- terra::rast(local_file)
+
+  # serve the result back
   local_df
 
 } # /function

R/vyskopis.R

@@ -1,42 +1,65 @@
 #' Vyskopis
 #'
-#' Terrain of the Czech Republic as a {raster} package object.
+#' Terrain of the Czech Republic as a {terra} package object.
 #'
-#' The function returns a raster file of either actual relief (values are meters above sea level) or rayshaded relief (created via highly recommended {rayshader} package).
+#' The function returns a raster file of either actual terrain (values are meters above sea level) or rayshaded relief..
 #'
-#' The raster is 5084 by 3403 cells, meaning each pixel is about 90 × 90 meters. It works the best at level of country or regions, at the level of a city or lower it may be somewhat grainy.
+#' The raster is created from EU DEM 1.1 file by Copernicus Land Monitoring service. The original file has pixel resolution 25×25 meters, which is too detailed for purposes of the package and was downsampled by factor of 4.
+#'
+#' The extent of the raster file is bouding box of the Czech Republic; this is a change to prior versions in order to better facilitate use of the raster in natural sciences context. To preserve compatibility optional argument `cropped` has been created, defaulting to `TRUE` (i.e. behaviour before v1.10.0).
 #'
 #' Due to package size constraints both versions are stored externally (and a working internet connection is required to use the package).
 #'
-#' The data is current to September 2016. Downloaded size of the rayshaded raster is 8.4 MB, actual raster is 31.4 MB.
+#' The data is current to year 2011 (but it is not expected to materially change over time). Downloaded size of the rasters is 52 MB.
+#'
+#' @param format Should the function return actual terrain (meters above sea level) or shaded relief (rayshaded). Allowed values are "actual" and "rayshaded".
 #'
-#' @param format Should the function return actual relief (meters above sea level) or shaded relief (rayshaded). Allowed values are "actual" and "rayshaded".
+#' @param cropped Should the raster provide data over Czech Republic's bounding box (cropped = FALSE) or just actual borders (cropped = TRUE). Defaults to TRUE to preserve compatiblity with earlier versions.
 #'
-#' @return \code{raster} package RasterLayer.
+#' @return \code{terra} package SpatRaster
 #'
-#' @source © ArcČR, ARCDATA PRAHA, ZÚ, ČSÚ, 2016 \url{https://www.arcdata.cz/produkty/geograficka-data/arccr-4}
+#' @source Copernicus Land Monitoring service, with funding by the European Union. \url{https://land.copernicus.eu/imagery-in-situ/eu-dem}
 #'
 #' @examples
 #' \donttest{
-#' library(raster)
+#' library(terra)
+#'
+#' # original extent - bounding box over Czech Republic
+#' original_extent <- vyskopis("rayshaded", cropped = FALSE)
+#'
+#' plot(original_extent, col = gray.colors(16))
 #'
-#' relief <- vyskopis("rayshaded")
+#' # add plot of country borders, for context
+#' plot(RCzechia::republika(),
+#'      border = "red",
+#'      col = NA,
+#'      add = TRUE)
+#'
+#' # cropped to size - default behaviour
+#' cropped_extent <- vyskopis("rayshaded")
+#'
+#' plot(cropped_extent, col = gray.colors(16))
 #'
-#' plot(relief, col = gray.colors(16))
 #' }
 #' @export
 
 
-vyskopis <- function(format = "rayshaded") {
+vyskopis <- function(format = "rayshaded", cropped = TRUE) {
   if (!is.element(format, c("actual", "rayshaded"))) {
     stop(paste(format, "is not a valid format; recognized values are \"actual\" or \"rayshaded\"."))
-  } # /if - valid resolution
+  } # /if - valid format
+  if (!is.element(cropped, c(TRUE, FALSE))) {
+    stop(paste(cropped, "is not a valid value for cropped; recognized values are \"TRUE\" or \"FALSE\"."))
+  } # /if - valid valid cropped
 
   if (format == "rayshaded") {
-    result <- .downloader("Vyskopis-stiny.rds")
+    result <- .downloader("vyskopis-shaded-dem.tif")
   } else {
-    result <- .downloader("Vyskopis-vyska.rds")
+    result <- .downloader("vyskopis-dem.tif")
   } # /if - download of result
 
+  # need to handle CRAN policy
+  if(cropped & !is.null(result)) result <- terra::crop(result, terra::vect(RCzechia::republika("low")), mask = TRUE)
+
   result
 } # /function

README.md

@@ -48,12 +48,12 @@ You can also get the latest development version by running `remotes::install_git
 administrative:
 
 * **republika**: borders of the Czech Republic as a polygon
-* **kraje**: 14 regions of the Czech Republic & Prague.  
-* **okresy**: 76 districts (LAU1 areas) of the Czech Republic + Prague (legally not *a district* but *the capital*).  
-* **orp_polygony** 205 municipalities with extended powers (in Czech: obce s rozšířenou působností) + Prague (legally not *a city* but *the capital*).  
-* **obce_polygony**: 6.258 municipalities of the Czech Republic.  
-* **obce_body** the same as obce_polygony, but centroids instead of polygons.  
-* **casti**: primarily 57 city parts of Prague, but also of other cities with defined parts (Brno, Ostrava and other).  
+* **kraje**: 14 regions of the Czech Republic & Prague  
+* **okresy**: 76 districts (LAU1 areas) of the Czech Republic + Prague (legally not *a district* but *the capital*)  
+* **orp_polygony** 205 municipalities with extended powers (in Czech: obce s rozšířenou působností) + Prague (legally not *a city* but *the capital*)  
+* **obce_polygony**: 6.258 municipalities of the Czech Republic  
+* **obce_body** the same as obce_polygony, but centroids instead of polygons  
+* **casti**: primarily 57 city parts of Prague, but also of other cities with defined parts (Brno, Ostrava and other)  
 * **senat_obvody**: 81 senate districts (volební obvody senátu)
 * **volebni_okrsky**: 14.761 general election districts (volební okrsky)
 * **zip_codes**: 2.671 ZIP code areas (poštovní směrovací čísla / PSČ)
@@ -67,7 +67,7 @@ natural:
 * **chr_uzemi**: protected natural areas (chráněná území)
 * **lesy**: woodland areas (more than 30 ha in area)
 * **KFME_grid**: KFME grid cells (faunistické čtverce)
-* **vyskopis**: terrain of the Czech republic as a {raster} package object
+* **vyskopis**: terrain of the Czech republic as a {terra} package raster object
 
 
 All objects are implemented as functions returning data frames, so must be followed by brackets (i.e. `hranice <- republika()`).

data-raw/digest-rasters.R

@@ -0,0 +1,35 @@
+# spojí DEM data z https://land.copernicus.eu/imagery-in-situ/eu-dem/eu-dem-v1.1
+# a ořízne dle potřeby
+
+library(terra)
+library(dplyr)
+library(sf)
+
+dolni <- terra::rast("./data-raw/eu_dem_v11_E40N20.TIF")
+horni <- terra::rast("./data-raw/eu_dem_v11_E40N30.TIF")
+
+celek <- terra::merge(horni, dolni)
+
+# okolí republiky (ne jej)
+maska <- RCzechia::republika() %>%
+  st_bbox() %>%
+  st_as_sfc() %>%
+  st_transform(st_crs(celek))
+
+cast <- crop(celek, maska, mask = T)
+
+cast <- terra::aggregate(cast, fact = 4) # 25 m >> 100 m pixel size
+
+# https://dieghernan.github.io/202210_tidyterra-hillshade/
+slope <- terrain(cast, "slope", unit = "radians")
+aspect <- terrain(cast, "aspect", unit = "radians")
+hill <- shade(slope, aspect, 30, 270)
+
+# normalize names
+names(hill) <- "shades"
+
+cast <- terra::project(cast, "EPSG:4326")
+hill <- terra::project(hill, "EPSG:4326")
+
+writeRaster(cast, "./data-backup/vyskopis-dem.tif", overwrite = T)
+writeRaster(hill, "./data-backup/vyskopis-shaded-dem.tif", overwrite = T)

tests/testthat/test-1-relief.R

@@ -1,5 +1,6 @@
 library(dplyr)
 library(httr)
+library(terra)
 
 test_that("reliéf", {
 
@@ -9,25 +10,33 @@ test_that("reliéf", {
   expect_message(vyskopis(), "internet") # zpráva o chybějícím internetu
   Sys.setenv("NETWORK_UP" = TRUE)
 
-  expect_s4_class(vyskopis(), "RasterLayer")
-  expect_s4_class(vyskopis("actual"), "RasterLayer")
-  expect_s4_class(vyskopis("rayshaded"), "RasterLayer")
+  expect_s4_class(vyskopis(), "SpatRaster")
+  expect_s4_class(vyskopis("actual"), "SpatRaster")
+  expect_s4_class(vyskopis("rayshaded"), "SpatRaster")
+  expect_s4_class(vyskopis("actual", FALSE), "SpatRaster")
+  expect_s4_class(vyskopis("rayshaded", FALSE), "SpatRaster")
 
-  # test velikosti
-  expect_equal(vyskopis()@ncols, 5084) # sloupce jsou všechny
-  expect_equal(vyskopis("actual")@ncols, 5084) # sloupce jsou všechny
-  expect_equal(vyskopis("rayshaded")@ncols, 5084) # sloupce jsou všechny
+  # test rozsahu
+  expect_equal(terra::expanse(vyskopis()), units::drop_units(st_area(republika("low"))), tolerance = 1/100)
+  expect_equal(terra::expanse(vyskopis("actual")), units::drop_units(st_area(republika("low"))), tolerance = 1/100)
+  expect_equal(terra::expanse(vyskopis("rayshaded")), units::drop_units(st_area(republika("low"))), tolerance = 1/100)
 
-  expect_equal(vyskopis()@nrows, 3403) # řádky jsou všechny
-  expect_equal(vyskopis("actual")@nrows, 3403) # řádky jsou všechny
-  expect_equal(vyskopis("rayshaded")@nrows, 3403) # řádky jsou všechny
+  # oříznutý raster je menší než surový
+  expect_gt(terra::expanse(vyskopis(cropped = F)), terra::expanse(vyskopis(cropped = T)))
+  expect_gt(terra::expanse(vyskopis("rayshaded", cropped = F)), terra::expanse(vyskopis("rayshaded", cropped = T)))
+  expect_gt(terra::expanse(vyskopis("actual", cropped = F)), terra::expanse(vyskopis("actual", cropped = T)))
 
   # test projekce - WGS84 pure & unadultered
-  expect_equal(projection(crs(vyskopis())), "+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs +towgs84=0,0,0")
-  expect_equal(projection(crs(vyskopis("actual"))), "+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs +towgs84=0,0,0")
-  expect_equal(projection(crs(vyskopis("rayshaded"))), "+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs +towgs84=0,0,0")
+  expect_equal(st_crs(vyskopis())$input, "WGS 84")
+  expect_equal(st_crs(vyskopis("actual"))$input, "WGS 84")
+  expect_equal(st_crs(vyskopis("rayshaded"))$input, "WGS 84")
+  expect_equal(st_crs(vyskopis("actual", FALSE))$input, "WGS 84")
+  expect_equal(st_crs(vyskopis("rayshaded", FALSE))$input, "WGS 84")
 
   # očekávaná chyba
   expect_error(vyskopis("bflm")) # neznámé rozlišení - očekávám actual / rayshaded
+  expect_error(vyskopis("actual", cropped = "bflm")) # cropped usí být boolean
+  expect_error(vyskopis("actual", cropped = NA))
+
 })
 

vignette.Rmd.orig

@@ -277,17 +277,17 @@ This example covers the second option.
 library(RCzechia)
 library(ggplot2)
 library(dplyr)
-library(raster, exclude = "select") # beware of colision with dplyr::select!!
+library(terra)
 
-# ggplot does not play nice with {raster} package; a data frame is required
-relief <- vyskopis("rayshaded") %>%
-  as("SpatialPixelsDataFrame") %>% # old style format - {sp}
-  as_tibble()
+# terrain cropped to "Czechia proper"
+relief <- vyskopis("rayshaded", cropped = TRUE)
 
 # report results
 ggplot() +
-  geom_raster(data = relief, aes(x = x, y  = y, alpha = -raytraced), # relief
-              fill = "gray30",  show.legend = F) + # no legend is necessary
+  tidyterra::geom_spatraster(data = relief) +
+  scale_fill_gradientn(colors = hcl.colors(50, "Grays"), # 50 shades of Gray...
+                       na.value = NA,
+                       guide = "none") +
   geom_sf(data = subset(RCzechia::reky(), Major == T), # major rivers
           color = "steelblue", alpha = .7) +
   labs(title = "Czech Rivers & Their Basins") +