Prepare Windows for GDAL

Since GDAL is used for the modifications of the EU-DEM dataset it needs to be defined for windows, where GDAL is located. In this case it is part of the OSGeo4W64 installation.

library(magrittr)
## Define for windows where GDAL is located
gdalpath = "C:/OSGeo4W64/bin"
if(.Platform$OS.type == "windows"){
  gdal.dir <- shortPathName(gdalpath)
  gdalwarp <- paste0(gdal.dir, "/gdalwarp.exe")
}else{
  gdalwarp = "gdalwarp"
}

Define shape of the study area by administrative and natural units

library(sf)
sh <- st_read("../data/raw_data/vector/vg250_0101.utm32s.shape.ebenen/vg250_ebenen/VG250_LAN.shp") %>%
  dplyr::filter(is.element(GEN,c("Hamburg" ,"Schleswig-Holstein")),
                GF == 4)
#> Reading layer `VG250_LAN' from data source `E:\Scales of Transformation\GitLab\b1\data\raw_data\vector\vg250_0101.utm32s.shape.ebenen\vg250_ebenen\VG250_LAN.shp' using driver `ESRI Shapefile'
#> Simple feature collection with 35 features and 23 fields
#> geometry type:  MULTIPOLYGON
#> dimension:      XY
#> bbox:           xmin: 280371.1 ymin: 5235856 xmax: 921292.4 ymax: 6106244
#> epsg (SRID):    NA
#> proj4string:    +proj=utm +zone=32 +ellps=GRS80 +units=m +no_defs

marsh <- st_read("../data/raw_data/vector/Natural_space/nat09_WGS84.shp") %>%
  dplyr::filter(grepl("arsch|Nordfrisische|Geestinsel",NAME)) %>% 
  st_transform(st_crs(sh))
#> Reading layer `nat09_WGS84' from data source `E:\Scales of Transformation\GitLab\b1\data\raw_data\vector\Natural_space\nat09_WGS84.shp' using driver `ESRI Shapefile'
#> Simple feature collection with 516 features and 8 fields
#> geometry type:  POLYGON
#> dimension:      XY
#> bbox:           xmin: 5.865461 ymin: 47.26964 xmax: 15.03949 ymax: 55.06563
#> epsg (SRID):    4326
#> proj4string:    +proj=longlat +datum=WGS84 +no_defs

sh_nomarsh <- st_difference(st_union(sh),st_union(marsh)) %>% 
  st_cast("POLYGON") %>% 
  .[order(st_area(.),decreasing = TRUE)[1:2]]

st_write(obj = sh,
         dsn = "../data/derived_data/vector/sh.shp",
         delete_dsn = TRUE)
#> Deleting source `../data/derived_data/vector/sh.shp' using driver `ESRI Shapefile'
#> Writing layer `sh' to data source `../data/derived_data/vector/sh.shp' using driver `ESRI Shapefile'
#> features:       2
#> fields:         23
#> geometry type:  Multi Polygon
st_write(obj = sh_nomarsh,
         dsn = "../data/derived_data/vector/sh_nomarsh.shp",
         delete_dsn = TRUE)
#> Deleting source `../data/derived_data/vector/sh_nomarsh.shp' using driver `ESRI Shapefile'
#> Writing layer `sh_nomarsh' to data source `../data/derived_data/vector/sh_nomarsh.shp' using driver `ESRI Shapefile'
#> features:       2
#> fields:         0
#> geometry type:  Polygon

Modify EU-DEM using gdalwarp

Using the following command the elevation model is cropped and reprojected

file.remove("../data/derived_data/raster/dem.tif")
#> [1] TRUE
system(paste(gdalwarp,
             "-s_srs EPSG:3035",
             "-t_srs EPSG:25832",
             "-r cubicspline",
             "-cutline ../data/derived_data/vector/sh_nomarsh.shp",
             "-crop_to_cutline",
             "../data/raw_data/raster/eu_dem_v11_E40N30/eu_dem_v11_E40N30.tif",
             "../data/derived_data/raster/dem.tif"))

library(raster)
plot(raster("../data/derived_data/raster/dem.tif"))