Skip to content

Commit

Permalink
Breaking changes
Browse files Browse the repository at this point in the history 
-AHN1-AHN3 is currently not supported to be used. Please change it to 'AHN' to always use the latest version made available by PDOK. Currently that is AHN4.
-The default extraction method to receive the AHN elevation is changed to 'simple' (nearest) instead of bilinear. Please see this [article](https://gisgeography.com/raster-resampling/) that explains the difference.
-In this release downloading point clouds is not supported.
-It is now possible to retrieve AHN elevation data at exact intersection points of 2 or more AHN sheets.

Other changes
-AHN4 is the latest full AHN version made available by PDOK. This release now supports this AHN version.
-in `ahn_area()` introduced the ability to only download AHN sheets using the `sheets` parameter.
-Introduced the function `ahn_sheets_info()` to retrieve information about the AHN sheets and which URLs are used for download.
-Removed the ability to use older AHN versions (AHN1-AHN3). It is in the planning to reintroduce this in a future release.
-Removed the ability to download point cloud data from the AHN. This will possibly be reintroduced in a future release.
-Adjusted URLs so that AHN datasets hosted by PDOK can be retrieved again.
-The code now looks at source URLs from data frames instead of hard coded references in the code.
-Efficient checks are built in to see if the AHN version called is currently supported.
-Reorganised code and removed deprecated/unsupported functions.
-Removed all dependencies of the deprecated 'rgdal' package. 'Terra' package is now used instead of 'raster'.
-Updated all the dependent package to the the latest versions.
-improved R syntax throughout the whole code.
  • Loading branch information
@Jellest
Jellest committed Apr 10, 2024
1 parent 74b234c commit 003547c
Show file tree
Hide file tree
Showing 30 changed files with 658 additions and 1,082 deletions.
190 changes: 129 additions & 61 deletions R/ahn_area.R
View file

Large diffs are not rendered by default.

12 changes: 0 additions & 12 deletions R/ahn_letters.R
View file

This file was deleted.

143 changes: 82 additions & 61 deletions R/ahn_point.R
View file
Original file line number Diff line number Diff line change
@@ -1,97 +1,118 @@
#'@title AHN elevation point
#'@description Get elevation of specific point location.
#'
#'Requires the X and Y coordinates as input. AHN data is obtained from the AHN3 (default), AHN2 or AHN 1 from the available resolutions. Default resolution is always the highest resolution (smallest number).
#'Requires the X and Y coordinates as input. AHN data is obtained from the AHN in the available resolutions. Default resolution is always the highest resolution (smallest number).
#'
#'AHN3 and AHN2 DEM (Digital Elevation Models) are available for the Digital Surface Model (DSM) and Digital Terrain Model (DTM). AHN1 only has the DTM. THE DTM of the AHN2 has an interpolated and a non-interpolate version.
#'Available for the Digital Surface Model (DSM) and Digital Terrain Model (DTM).
#'
#'You can download the AHN data using the WCS method (default, sheets.method = FALSE) returning a GeoTIFF float32 format. The sheets method (sheets.method = TRUE) returns a regular raster GeoTIFF output file that is extracted from the PDOK sheets. The WCS method is recommended if only a few AHN elevation points need to be extracted. The sheets method always requires more data to be downloaded to the client but may be more efficient if many elevatiomns need to be retrieved from a small area. Choosing your method depends on speed and your desired output format. See documentation for all available parameters.
#'@param name Optional. Give a name of the specified area. This name will be used in the output file names.
#'@param X Required. X coordinate in RD New or WGS84 (LON).
#'@param Y Required. Y coordinate in RD New or WGS84 (LAT).
#'@param AHN Default 'AHN3'. Set to 'AHN1', 'AHN2', or 'AHN3'.
#'@param dem Default 'DSM'. Choose type of Digital Elevation Model. 'DSM' or 'DTM'. AHN1 only has 'DTM'.
#'@param resolution Default 0.5 meters for AHN2/AHN3, 5 meters for AHN1. Choose resolution of AHN in meters. AHN3 and AHN2 both have 0.5 and 5 meters. AHN1 has 5 and 100 m.
#'@param interpolate Default TRUE. Only applicable for AHN2 DTM. If true, it gets the interpolated version of the AHN2.
#'@param output.dir Optional but unnecessary. Set location of output raster files. Leaving blank (default) will make all output point files be temporary files. This output directory excludes the location of the AHN sheets which is depicted with the \code{sheets.location} parameter.
#'@param LONLAT Optional. Default FALSE. Set to TRUE if X and Y are in Longitude and Latitude format. Output will be in RD New format.
#'@param extract.method Default 'bilinear'. Choose 'bilinear or 'simple'. Intersection is done using [\code{extract()}](https://www.rdocumentation.org/packages/raster/versions/3.1-5/topics/extract) function from the \code{raster} package.
#'@param decimals Default 2. Decide number of decimal places of output elevations.
#'@param sheets.method Default FALSE. FALSE downloads AHN area through the faster WCS method. Output is 32float GeoTIFF file.TRUE downloads AHN area through the available GeoTIFF AHN sheets available on [PDOK](http://esrinl-content.maps.arcgis.com/apps/Embed/index.html?appid=a3dfa5a818174aa787392e461c80f781).
#'@param sheets.location Default is the 'AHN_sheets' directory in the working directory. Set directory where all the AHN sheets are loaded when pre-existing sheets will be used or when new sheets will be stored. When loading existing files, always use the correct directory structure and capitalization within the selected directory. Example directory structure when this parameter is set to e.g. 'myFolder': 'myFolder/AHN_sheets/AHN3/DSM' or 'myFolder/AHN_sheets/AHN2/DTM'. Only use extracted files in their original name after download.
#'@param sheets.keep Default TRUE. Only applicable if \code{sheets.method} is set to TRUE and sheets were downloaded. Set to FALSE if you want to delete the downloaded sheet. It is recommended to keep the sheets if ahn elevation extarction will be followed.
#'@author Jelle Stuurman
#'You can download the AHN data using the WCS method (default, sheets.method = FALSE) returning a GeoTIFF format. The sheets method (sheets.method = TRUE) returns a regular raster GeoTIFF output file that is extracted from the PDOK sheets. The WCS method is recommended if only a few AHN elevation points need to be extracted. The sheets method always requires more data to be downloaded to the client but may be more efficient if many elevatiomns need to be retrieved from a small area. Choosing your method depends on speed and your desired output format. See documentation for all available parameters.
#'@inheritParams ahn_area
#'@param decimals Default 2. Maximum Determines the number of decimal places in the output elevation.
#'@param extract.method Default 'simple'. Choose between 'simple' (nearest) and 'bilinear'. Intersection is done using [\code{extract()}](https://www.rdocumentation.org/packages/terra/versions/1.7-71/topics/extract) function from the \code{terra} package. See [this](https://gisgeography.com/raster-resampling/) article that explains the difference between the two methods.
#'@return AHN elevation in meters.
#'@export
ahn_point <- function(name = "AHNelevation", X, Y, AHN = "AHN3", dem = "DSM", resolution, interpolate = TRUE, output.dir, LONLAT = FALSE, extract.method = "bilinear", decimals = 2, sheets.method = FALSE, sheets.location, sheets.keep = TRUE){
loadNamespace("raster")
name_trim <- trim_name(name)
ahn_point <- function(name = "AHNelevation", X, Y, AHN = "AHN", dem, resolution, output.dir, LONLAT = FALSE, extract.method = "simple", decimals = 2, sheets.method = FALSE, sheets.dir, sheets.keep = TRUE) {
interpolate <- TRUE
loadNamespace("terra")
# check for selected AHN
AHN <- check_ahn_version(AHN)

if (missing(dem) == TRUE) {
dem <- ""
}
name_trim <- gsub(" ", "_", name)

#set tmp folder if applicable or create output and directory
if(missing(output.dir) == TRUE){
if (missing(output.dir) == TRUE) {
output.dir <- tempdir()
} else {
if(!dir.exists(output.dir) == TRUE){
if (dir.exists(output.dir) == FALSE) {
dir.create(output.dir)
print(paste0("'",output.dir, "' directory was not found and was created."))
print(paste(output.dir, "directory was not found and was created.", sep = " "))
}
#print(output.dir)
}

#selected AHN
if(tolower(AHN) != "ahn1" && tolower(AHN) != "ahn2" && tolower(AHN) != "ahn3"){
stop("No correct AHN is provided. Please select 'AHN3' or 'AHN2'")
} else {
AHN <- toupper(AHN)
#get resolution
if (missing(resolution) == TRUE) {
resolution <- ""
}
my_resolution <- get_resolution(AHN = AHN, resolution = resolution)

#set missing resolution
if(missing(resolution) == TRUE){
if(AHN == "AHN3" || AHN == "AHN2"){
resolution = 0.5
} else if(AHN =="AHN1"){
resolution = 5
}
}
#get dem type
dem <- get_dem(AHN = AHN, resolution = my_resolution, dem = dem, interpolate = interpolate)

#get and create a point
my_point <- generate_ahn_point(name = name_trim, X = X, Y = Y, LONLAT = LONLAT, resolution = resolution)
my_point <- generate_ahn_point(name = name_trim, X = X, Y = Y, LONLAT = LONLAT, resolution = my_resolution$res)

#get AHN elevation
if(sheets.method == TRUE){
#get AHN data
bladIndex.sf <- get_bladindex(AHN = AHN, dem = dem, resolution = my_resolution$res)
if (sheets.method == FALSE) {
##get elevation through WCS method (fast)
wcs_source <- bladIndex.sf$wcs_url[1]
print(bladIndex.sf$wcs_url[1])
wcs_url <- create_wcs_url(type = "point", bbox = my_point$bbox, AHN = AHN, dem = dem, resolution = my_resolution, interpolate = interpolate, wcs = wcs_source)
raster_data <- download_wcs_raster(wcsUrl = wcs_url, name = name_trim, area = ahn_area, AHN = AHN, dem = tolower(dem), resolution = my_resolution, output.dir = output.dir, interpolate = interpolate, type = "point")
my_elevation <- extract_elevation(raster_data$data, my_point$point, extract.method = extract.method)
} else if (sheets.method == TRUE) {
#download AHN sheets and get data (slow)

#set AHN sheets location
if(missing(sheets.location) == TRUE){
sheets.location <- getwd()
print(paste0("The AHN sheets are loaded from or downloaded in: ", sheets.location, "/", default.sheets.dir, "/", AHN, "/", dem))
if (missing(sheets.dir) == TRUE) {
#sheets directory
if (output.dir == tempdir()) {
sheets.dir <- getwd()
} else {
sheets.dir <- output.dir
}
}
if (sheets.dir == tempdir()) {
stop("Due to the size of these AHN sheets, this script does not allow to put these in the RAM memory. Please adjust the location of the sheets.dir.")
}
if (grepl(default.sheets.dir, sheets.dir, fixed = TRUE) == TRUE) {
ahn_sheets_directory <- sheets.dir
} else {
ahn_sheet_directory <- paste(sheets.location, AHN, dem, sep="/")
print(sprintf("The AHN sheets are loaded from or downloaded in: %s. If no AHN sheet in the correct directory or if no correct name of AHN sheet is found, sheet will be downloaded. For first use it is recommended to use the default output directory.", ahn_sheet_directory))
if (dir.exists(sheets.dir) == FALSE) {
dir.create(sheets.dir, showWarnings = FALSE)
}
ahn_sheets_directory <- paste(sheets.dir, default.sheets.dir, sep = "/")
}
if (dir.exists(ahn_sheets_directory) == FALSE) {
dir.create(ahn_sheets_directory, showWarnings = FALSE)
}
print(sprintf("The AHN sheets are loaded from or downloaded in: %s. If no AHN sheet in the correct directory or if no correct name of AHN sheet is found, sheet will be downloaded. For first use it is recommended to use the default output directory.", ahn_sheets_directory))

#create area
ahn_area <- create_area(bbox = my_point$bbox, LONLAT = LONLAT, type = "point")
ahn_area <- create_area(radius = "", bbox = my_point$bbox, LONLAT = LONLAT, type = "point")

#download AHN sheets and get data (slow)
bladnrs <- find_ahn_sheets(name, area = ahn_area, type = "point", bladIndex = bladIndex.sf)
#get AHN area
raster_data <- get_ahn_sheets(name = name_trim, area = ahn_area, type = "point", AHN = AHN, dem = dem, resolution = resolution, radius = "", interpolate = interpolate, output.dir = output.dir, sheets.keep = sheets.keep, sheets.location = sheets.location)
sheets_df <- data.frame(kaartBlad = character(), filePath = character(), dwnld = logical(), stringsAsFactors = FALSE)
for (b in bladnrs) {
url <- bladIndex.sf$atom_url[bladIndex.sf$kaartbladNr == b]
path_sheet <- download_ahn_sheets(name = name_trim, AHN = AHN, dem = dem, url = url, output.dir = output.dir, sheets.dir = ahn_sheets_directory, sheets.keep = sheets.keep)
sheets_df <- rbind(sheets_df, path_sheet)
}
raster_data <- merge_and_crop_ahn(name = name_trim, sheets = sheets_df, area = ahn_area, AHN = AHN, dem = tolower(dem), resolution = my_resolution$res_name, output.dir = output.dir)

#get elevation
my_elevation <- extract_elevation(raster_data$data, my_point$point, extract.method = extract.method)
} else {
#retrieve data through WCS (fast)
my_resolution <- get_resolution(AHN = AHN, resolution = resolution)

my_url <- create_wcs_url(type = "point", bbox = my_point$bbox, AHN = AHN, dem = dem, resolution = my_resolution, interpolate = interpolate)
raster_data <- download_wcs_raster(wcsUrl = my_url, name = name_trim, AHN = AHN, dem = tolower(dem), resolution = my_resolution, output.dir = output.dir, interpolate = interpolate, type = "point")
my_elevation <- extract_elevation(raster_data$data, my_point$point, extract.method = extract.method)
stop("Ncorrect value is provided for the sheets.method parameter. Please set it to 'TRUE' or 'FALSE'.")
}

my_elevation <- format(round(my_elevation, decimals), nsmall = decimals)
print(paste("Elevation of ", name , ": ", my_elevation, " m.", sep=""))
my_elevation <- as.numeric(my_elevation)
# if(output.dir == tempdir()){
# unlink(raster_data$fileName)
# }
return (my_elevation)
if (is.na(my_elevation)) {
print(paste("No elevation is available for this point in the ", AHN, " ", dem, ".", sep = ""))
} else {
my_elevation <- format(round(my_elevation, decimals), nsmall = decimals)
print(paste("Elevation of ", name, ": ", my_elevation, " m.", sep = ""))
my_elevation <- as.numeric(my_elevation)
}
if (output.dir == tempdir()) {
base::unlink(raster_data$fileName)
}
if (sheets.keep == FALSE) {
delete_ahn_sheet(sheets_df)
}
return(my_elevation)
}
27 changes: 27 additions & 0 deletions R/ahn_sheets_info.R
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,27 @@
#'@title AHN sheets info
#'@description get polygon features that desribe the AHN sheets
#'@inheritParams ahn_area
#'@export
#'@author Jelle Stuurman
ahn_sheets_info <- function(AHN, dem, output.dir) {
if (missing(AHN) == TRUE) {
stop("No value provided for the AHN prameter. Please enter a valid AHN parameter.")
}
check_ahn_version(AHN)

if (missing(resolution)) {
resolution <- ""
my_resolution <- get_resolution(AHN = AHN, resolution = resolution)
}

if (missing(dem) == TRUE) {
stop("No value provided for the dem prameter. Please enter a valid dem parameter.")
}

ahn_bi <- get_bladindex(AHN = AHN, dem = dem, resolution = my_resolution$res)

if (missing(output.dir) == FALSE) {
sf::st_write(obj = ahn_bi, dsn = paste(output.dir, "/", paste0(AHN, "_", "bladindex.gpkg"), sep = ""), append = FALSE)
}
return(ahn_bi)
}
36 changes: 36 additions & 0 deletions R/check_ahn_version.R
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,36 @@
#'@description check the version of AHN that is called and whehter it is available within this package.
#'@inheritParams ahn_area
#'@noRd
check_ahn_version <- function(AHN) {
AHN <- toupper(AHN)
note <- FALSE
if (AHN == "AHN" || AHN == latest_pdok_version) {
if (AHN == "AHN") {
note <- TRUE
}
xml_doc <- xml2::read_xml(pdok_wcs_url)

# Extract the abstract content
abstract <- xml2::xml_text(xml2::xml_find_first(x = xml_doc, xpath = ".//ows:Abstract", ns = xml2::xml_ns(xml_doc)))

# Determine AHN version
ahn_version <- base::regmatches(abstract, base::regexpr("(?<=versie\\s)\\d+", abstract, perl = TRUE))
this_ahn <- paste0("AHN", ahn_version)

if (!toupper(this_ahn) %in% pdok_versions) {
stop(" The AHN version found is currently not supoorted at the moment with this package.")
} else {
AHN <- toupper(this_ahn)
if (note == TRUE) {
message(paste("The AHN version found through PDOK is:", AHN))
}
return(AHN)
}
} else if (AHN %in% pdok_versions) {
return(AHN)
} else if (AHN %in% operational_ahns_in_package) {
stop("Currently there is no support to use other AHN versions than the AHN made available by PDOK.")
} else {
stop(" The provided AHN is not found. Please try another name.")
}
}
Loading

0 comments on commit 003547c

Please sign in to comment.