From e7efeeede6987ddc883d00aa36d2a840859977ff Mon Sep 17 00:00:00 2001 From: Jelle Stuurman Date: Tue, 16 Apr 2024 22:37:16 +0200 Subject: [PATCH 1/4] added min/max coordinates --- data-raw/DATASET.R | 12 ++++++++++-- 1 file changed, 10 insertions(+), 2 deletions(-) diff --git a/data-raw/DATASET.R b/data-raw/DATASET.R index 27ab2bc..c041b3a 100644 --- a/data-raw/DATASET.R +++ b/data-raw/DATASET.R @@ -1,4 +1,5 @@ ## code to prepare `DATASET` dataset goes here +#'@title create global variables #'@noRd define_globals <- function() { pdok_baseurl <- "https://service.pdok.nl" @@ -11,6 +12,13 @@ define_globals <- function() { proj_rd <- "epsg:28992" proj_wgs <- "epsg:4326" + RD_min <- list("X" = 482.06, "Y" = 306602.42) + RD_max <- list("X" = 284182.97, "Y" = 637049.52) + + WGS_min <- list("LON" = 3.2, "LAT" = 50.75) + WGS_max <- list("LON" = 7.22, "LAT" = 53.7) + + AHN_versions <- c("AHN5", "AHN4", "AHN3", "AHN2", "AHN1") AHN_05res_versions <- c("AHN4", "AHN3", "AHN2") latest_pdok_version <- "AHN4" @@ -31,7 +39,7 @@ define_globals <- function() { if ("AHN4" %in% operational_ahns_in_package) { ahn_DSM05_bladIndex <- sf::st_read("data-raw/AHN_bladIndexes.gpkg", layer = "AHN_DSM05_bladIndex_rd") ahn_DTM05_bladIndex <- sf::st_read("data-raw/AHN_bladIndexes.gpkg", layer = "AHN_DTM05_bladIndex_rd") - ahn4_bladIndex <- print("See 'ahn_DSM05_bladIndex' or 'ahn_DTM05_bladIndex' variables.") + ahn4_bladIndex <- "See 'ahn_DSM05_bladIndex' and 'ahn_DTM05_bladIndex' variables." } else { ahn_DSM05_bladIndex <- NULL ahn_DTM05_bladIndex <- NULL @@ -52,7 +60,7 @@ define_globals <- function() { } else { ahn1_bladIndex <- NULL } - usethis::use_data(pdok_baseurl, pdok_wcs_url, ahn_DSM05_bladIndex, ahn_DTM05_bladIndex, ahn1_bladIndex, ahn2_bladIndex, ahn3_bladIndex, ahn4_bladIndex, ahn5_bladIndex, epsg_rd, epsg_wgs, proj_rd, proj_wgs, list_ahn_letters, AHN_versions, AHN_05res_versions, latest_pdok_version, operational_ahns_in_package, pdok_versions, non_pdok_versions, default.output.dir, default.sheets.dir, overwrite = TRUE, internal = TRUE) + usethis::use_data(pdok_baseurl, pdok_wcs_url, ahn_DSM05_bladIndex, ahn_DTM05_bladIndex, ahn1_bladIndex, ahn2_bladIndex, ahn3_bladIndex, ahn4_bladIndex, ahn5_bladIndex, epsg_rd, epsg_wgs, proj_rd, proj_wgs, RD_min, RD_max, WGS_min, WGS_max, list_ahn_letters, AHN_versions, AHN_05res_versions, latest_pdok_version, operational_ahns_in_package, pdok_versions, non_pdok_versions, default.output.dir, default.sheets.dir, overwrite = TRUE, internal = TRUE) } define_globals() From 36acfd9c7962b32e9150fb0f9b8850c4bc8f82b8 Mon Sep 17 00:00:00 2001 From: Jelle Stuurman Date: Tue, 16 Apr 2024 22:41:17 +0200 Subject: [PATCH 2/4] improved code to always get rectified grid or give warnings --- R/ahn_area.R | 27 ++++++------ R/ahn_point.R | 18 ++++---- R/create_area.R | 42 ++++++++++++------- R/create_bbox_polygon.R | 19 +++++++-- R/find_ahn_sheets.R | 15 ++++--- R/get_rectified_coordinates.R | 30 +++++++++++++ ..._grid.R => get_rectified_grid_for_point.R} | 42 ++++++------------- R/merge_and_crop_ahn.R | 4 +- 8 files changed, 117 insertions(+), 80 deletions(-) create mode 100644 R/get_rectified_coordinates.R rename R/{get_rectified_grid.R => get_rectified_grid_for_point.R} (50%) diff --git a/R/ahn_area.R b/R/ahn_area.R index 465540a..e059036 100644 --- a/R/ahn_area.R +++ b/R/ahn_area.R @@ -1,13 +1,13 @@ #'@title AHN raster area #'@description Get AHN raster area of a specified area. #' -#'AHN Area is retrieved through a circle (x, y and radius), BBOX (x, y and radius, OR bbox coordinates), or own geometry (polygon). +#'AHN Area is retrieved through a circle (x, y and radius), BBOX (x, y and radius, OR BBOX coordinates), or own geometry (polygon). #' #'AHN data is obtained from the AHN3 (default), AHN2 or AHN 1 from the available resolutions. #' #'Available for the Digital Surface Model (DSM) and Digital Terrain Model (DTM). #' -#'You can download the AHN data using the WCS method (default, \code{sheets.method = FALSE}) returning a GeoTIFF format. The sheets method (\code{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. +#'You can download the AHN data using the WCS method (default, \code{sheets.method = FALSE}) returning a GeoTIFF format. The sheets method (\code{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 elevations 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). @@ -21,8 +21,8 @@ #'@param sheets Only required when only full AHN sheets need to be downloaded. To do so, fill the list() object with strings of the sheet numbers. #'@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 sheets.method Default FALSE. When TRUE, it downloads AHN area through the available GeoTIFF AHN sheets available on [PDOK](https://www.pdok.nl/atom-downloadservices/-/article/actueel-hoogtebestand-nederland-ahn). -#'@param sheets.dir Default is the 'AHN_sheets' directory in the working directory or \code{output.dir} (when defined). In this sheets directory all the AHN sheets are downloaded. If sheets already exists in this directory, they will be reused and no redownload will take place. Within this sheets directory, a directory structure will be created automatically: if `sheets.dir` is set to 'myFolder', all sheets will be downloaded in their apropriate AHN version and DEM folder. Eg. '\code{output.dir}/myFolder/AHN_sheets/AHN4/DSM'. You want to use this parameter if you want to have a central location where all the AHN sheets are stored so that they can be accessed elsewhere for other purposes. It is recommended to always use the their original file name after download. Due to the size of these GeoTiff sheets, it is not allowed to download this is the RAM memory (tempdir). -#'@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. +#'@param sheets.dir Default is the 'AHN_sheets' directory in the working directory or \code{output.dir} (when defined). In this sheets directory all the AHN sheets are downloaded. If sheets already exists in this directory, they will be reused and no redownload will take place. Within this sheets directory, a directory structure will be created automatically: if `sheets.dir` is set to 'myFolder', all sheets will be downloaded in their appropriate AHN version and DEM folder. Eg. '\code{output.dir}/myFolder/AHN_sheets/AHN4/DSM'. You want to use this parameter if you want to have a central location where all the AHN sheets are stored so that they can be accessed elsewhere for other purposes. It is recommended to always use the their original file name after download. Due to the size of these GeoTiff sheets, it is not allowed to download this is the RAM memory (tempdir). +#'@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 extraction will be followed. #'@author Jelle Stuurman #'@return GeoTIFF file of AHN area #'@export @@ -43,7 +43,7 @@ ahn_area <- function(name = "AHNarea", output.dir, X, Y, radius, bbox, polygon, output.dir <- tempdir() } else { if (dir.exists(output.dir) == FALSE) { - dir.create(output.dir) + dir.create(output.dir, recursive = TRUE) print(paste(output.dir, "directory was not found and was created.", sep = " ")) } #print(output.dir) @@ -69,7 +69,7 @@ ahn_area <- function(name = "AHNarea", output.dir, X, Y, radius, bbox, polygon, if (length(sheets) == 0) { #get and create area - ahn_area <- create_area(X = X, Y = Y, radius = radius, bbox = bbox, polygon = polygon, LONLAT = LONLAT, type = "raster") + ahn_area <- create_area(name = name_trim, X = X, Y = Y, resolution = my_resolution, radius = radius, bbox = bbox, polygon = polygon, LONLAT = LONLAT, type = "raster") } else { if (!is.list(sheets)) { stop("No list is provided as the input for the 'sheets' parameter. Please use a list() with the needed sheet nrs.") @@ -84,9 +84,10 @@ ahn_area <- function(name = "AHNarea", output.dir, X, Y, radius, bbox, polygon, if (sheets.method == FALSE) { #use WCS method and get data (fast) wcs_source <- bladIndex.sf$wcs_url[1] - wcs_url <- create_wcs_url(bbox = ahn_area$bbox, type = "area", AHN = AHN, resolution = my_resolution, dem = dem, interpolate = interpolate, wcs = wcs_source) + wcs_url <- create_wcs_url(bbox = ahn_area$rectified_bbox, type = "area", AHN = AHN, resolution = my_resolution, dem = dem, interpolate = interpolate, wcs = wcs_source, LONLAT = LONLAT) raster_data <- download_wcs_raster(wcsUrl = wcs_url, name = name_trim, AHN = AHN, dem = tolower(dem), radius = radius, resolution = my_resolution, interpolate = interpolate, output.dir = output.dir, type = "raster") ahn_data <- terra::mask(x = raster_data$data, mask = ahn_area$area, filename = raster_data$fileName, overwrite = TRUE) + unlink(raster_data$data) output <- ahn_data if (output.dir == tempdir()) { base::unlink(ahn_data) @@ -114,13 +115,13 @@ ahn_area <- function(name = "AHNarea", output.dir, X, Y, radius, bbox, polygon, if (dir.exists(ahn_sheets_directory) == FALSE) { dir.create(ahn_sheets_directory, showWarnings = FALSE) } - print(sprintf("The AHN sheets are loaded from or downloaded into: %s. If no AHN sheet(s) are found in this directory or if no correct name of AHN sheet(s) are found, the AHN sheet(s) will be downloaded.", ahn_sheets_directory)) + print(sprintf("The AHN sheets are loaded from or downloaded into: %s. If no AHN sheet(s) are found in this directory or if no correct name of AHN sheet(s) are found, the AHN sheet(s) will be downloaded here.", ahn_sheets_directory)) #download AHN sheets and get data (slow) - sheets_df <- data.frame(kaartBlad = character(), filePath = character(), dwnld = logical(), stringsAsFactors = FALSE) + sheets_df <- data.frame(kaartBlad = character(), filePath = character(), downloaded = logical(), stringsAsFactors = FALSE) if (length(sheets) == 0) { #get ahn area from ahn sheets - bladnrs <- find_ahn_sheets(name = name_trim, area = ahn_area, type = "area", bladIndex = bladIndex.sf) + bladnrs <- find_ahn_sheets(name = name_trim, area = ahn_area$area, type = "area", bladIndex = bladIndex.sf) #print(bladnrs) for (b in bladnrs) { url <- bladIndex.sf$atom_url[bladIndex.sf$kaartbladNr == b] @@ -128,7 +129,7 @@ ahn_area <- function(name = "AHNarea", output.dir, X, Y, radius, bbox, polygon, 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) } - ahn_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) + ahn_data <- merge_and_crop_ahn(name = name_trim, sheets = sheets_df, area = ahn_area$area, AHN = AHN, dem = tolower(dem), resolution = my_resolution$res_name, output.dir = output.dir) if (LONLAT == TRUE) { warning("The input geometry was provided using Longitude and Latitude coordinates. The output is exported as a raster using the the RD New (epsg 28992) coordinate system.") @@ -148,13 +149,13 @@ ahn_area <- function(name = "AHNarea", output.dir, X, Y, radius, bbox, polygon, sheets_df <- rbind(sheets_df, path_sheet) } if (LONLAT == TRUE) { - warning("The input geometry was provided using Longitude and Latitude coordinates. The downloaded sheet(s) are geotif files with the RD New cordinate system.") + warning("The input geometry was provided using Longitude and Latitude coordinates. The downloaded sheet(s) are geotiff files with the RD New coordinate system.") } bladen <- data.frame(filePath = character(), stringsAsFactors = FALSE) for (d in seq_along(sheets_df$kaartBlad)) { bladen[d, "filePath"] <- sheets_df[d, "filePath"] if (sheets_df[d, "downloadedNow"] == TRUE) { - print(paste("The AHN sheet", sheets_df[d, "kaartBlad"], "in", sheets_df[d, "filePath"], "is downloaded succesfully.", sep = " ")) + print(paste("The AHN sheet", sheets_df[d, "kaartBlad"], "in", sheets_df[d, "filePath"], "is downloaded successfully.", sep = " ")) } } output <- bladen$filePath diff --git a/R/ahn_point.R b/R/ahn_point.R index abbf872..9384a07 100644 --- a/R/ahn_point.R +++ b/R/ahn_point.R @@ -5,7 +5,7 @@ #' #'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 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. +#'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 elevations 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. @@ -42,15 +42,15 @@ ahn_point <- function(name = "AHNelevation", X, Y, AHN = "AHN", dem, resolution, dem <- get_dem(AHN = AHN, resolution = my_resolution, dem = dem, interpolate = interpolate) #get and create a point - my_point <- get_rectified_grid(name = name_trim, X = X, Y = Y, LONLAT = LONLAT, resolution = my_resolution$res) + my_point <- get_rectified_grid_for_point(name = name_trim, X = X, Y = Y, LONLAT = LONLAT, resolution = my_resolution$res) #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] - 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") + wcs_url <- create_wcs_url(type = "point", bbox = my_point$rectified_bbox, AHN = AHN, dem = dem, resolution = my_resolution, interpolate = interpolate, wcs = wcs_source, LONLAT = LONLAT) + raster_data <- download_wcs_raster(wcsUrl = wcs_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) } else if (sheets.method == TRUE) { #download AHN sheets and get data (slow) @@ -81,23 +81,23 @@ ahn_point <- function(name = "AHNelevation", X, Y, AHN = "AHN", dem, resolution, 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(radius = "", bbox = my_point$bbox, LONLAT = LONLAT, type = "point") + ahn_area <- create_area(radius = "", bbox = my_point$rectified_bbox, resolution = my_resolution, LONLAT = LONLAT, type = "point") #download AHN sheets and get data (slow) - bladnrs <- find_ahn_sheets(name, area = ahn_area, type = "point", bladIndex = bladIndex.sf) + bladnrs <- find_ahn_sheets(name, area = ahn_area$area, type = "point", bladIndex = bladIndex.sf) #get AHN area - sheets_df <- data.frame(kaartBlad = character(), filePath = character(), dwnld = logical(), stringsAsFactors = FALSE) + sheets_df <- data.frame(kaartBlad = character(), filePath = character(), downloaded = 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) + raster_data <- merge_and_crop_ahn(name = name_trim, sheets = sheets_df, area = ahn_area$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 { - stop("Ncorrect value is provided for the sheets.method parameter. Please set it to 'TRUE' or 'FALSE'.") + stop("No correct value is provided for the sheets.method parameter. Please set it to 'TRUE' or 'FALSE'.") } if (is.na(my_elevation)) { diff --git a/R/create_area.R b/R/create_area.R index 08cabee..b41e732 100644 --- a/R/create_area.R +++ b/R/create_area.R @@ -1,10 +1,10 @@ #'@inheritParams ahn_area #'@noRd -create_area <- function(X, Y, radius, bbox, polygon, LONLAT = FALSE, type) { +create_area <- function(name, X, Y, radius, resolution, bbox, polygon, LONLAT = FALSE, type) { if (type == "point") { my_bbox <- bbox - my_bbox_area.sf <- create_bbox_polygon(my_bbox) - my_area.sf <- my_bbox_area.sf + my_bbox_area.sf <- create_bbox_polygon(name = name, bbox = my_bbox, resolution = resolution$res) + my_area.sf <- my_bbox_area.sf$rectified_area } else if (type == "raster" || type == "pc") { if (missing(X) == TRUE && missing(Y) == TRUE && missing(polygon) == TRUE && radius == "") { #create BBOX using BBOX coordinates @@ -13,18 +13,19 @@ create_area <- function(X, Y, radius, bbox, polygon, LONLAT = FALSE, type) { stop("4 coordinates are required: XMIN, YMIN, XMAX, YMAX.") } if (LONLAT == TRUE) { + #create spatial points in RD my_min <- create_spatialpoint(X = bbox$xmin, Y = bbox$ymin, LONLAT = LONLAT) my_max <- create_spatialpoint(X = bbox$xmax, Y = bbox$ymax, LONLAT = LONLAT) min_coords <- sf::st_coordinates(my_min) max_coords <- sf::st_coordinates(my_max) - my_bbox <- data.frame("xmin" = min_coords[1, "X"], "ymin" = min_coords[1, "Y"], "xmax" = max_coords[1, "X"], "ymax" = max_coords[1, "Y"]) + my_bbox <- list("xmin" = min_coords[1, "X"], "ymin" = min_coords[1, "Y"], "xmax" = max_coords[1, "X"], "ymax" = max_coords[1, "Y"]) } else { - my_bbox <- data.frame("xmin" = bbox[1], "ymin" = bbox[2], "xmax" = bbox[3], "ymax" = bbox[4]) + my_bbox <- list("xmin" = bbox[1], "ymin" = bbox[2], "xmax" = bbox[3], "ymax" = bbox[4]) } - my_bbox_area.sf <- create_bbox_polygon(my_bbox) - my_area.sf <- my_bbox_area.sf + my_bbox_area.sf <- create_bbox_polygon(name = name, bbox = my_bbox, resolution = resolution$res) + my_area.sf <- my_bbox_area.sf$rectified_area } else if (radius != "") { #create circle through buffer around a X,Y point if (radius == 0) { @@ -42,21 +43,30 @@ create_area <- function(X, Y, radius, bbox, polygon, LONLAT = FALSE, type) { sf::st_crs(my_area.sf, epsg_rd) my_bbox <- sf::st_bbox(my_area.sf, crs = epsg_rd) - my_bbox <- data.frame("xmin" = floor(my_bbox$xmin), "ymin" = floor(my_bbox$ymin), "xmax" = ceiling(my_bbox$xmax), "ymax" = ceiling(my_bbox$ymax)) - my_bbox_area.sf <- create_bbox_polygon(my_bbox) + my_bbox_area.sf <- create_bbox_polygon(name = name, bbox = my_bbox, resolution = resolution$res) + + # recreate circle based on rectified grid + new_radius <- (my_bbox_area.sf$rectified_bbox$xmax - my_bbox_area.sf$rectified_bbox$xmin) / 2 + if (new_radius != radius) { + warning(paste("The provided X, Y coordinates and radius do not return an area that matches the AHN rectified grid. The radius has therefore been adjusted with ", (new_radius - radius), "meters.", sep = " ")) + } + my_area.sf <- sf::st_buffer(my_point, dist = radius) + sf::st_crs(my_area.sf, epsg_rd) + if (bbox == TRUE) { #create bbox based on only radius print("Creating bbox from X,Y point and radius input.") - my_area.sf <- my_bbox_area.sf + my_area.sf <- my_bbox_area.sf$rectified_area } } else if (missing(X) == TRUE && missing(Y) == TRUE && bbox == FALSE && radius == "") { #load shape through polygon to create area - print(polygon) + warning("Using a polygon as area will most likely not return a raster that aligns with the rectified grid. If this raster will be used for elevation analysis, wrong values may be calculated.") + my_area.sf <- sf::st_as_sf(polygon) geometry_point <- sf::st_is(my_area.sf, "POINT") geometry_polygon <- sf::st_is(my_area.sf, "POLYGON") if (geometry_point == "POINT" || geometry_polygon == FALSE) { - stop("Geometry type is not a polygon. Please make sure the input given for the 'polygon' paramter is a polygon in the correct format.") + stop("Geometry type is not a polygon. Please make sure the input given for the 'polygon' parameter is a polygon in the correct format.") } print("Creating area from custom geometry.") if (LONLAT == TRUE) { @@ -65,9 +75,10 @@ create_area <- function(X, Y, radius, bbox, polygon, LONLAT = FALSE, type) { if (nrow(my_area.sf) != 1) { stop("The selected polygon has no or more than one feature. Add/reduce to one feature or use loop functionalities.") } + + #get bbox of polygon area my_bbox <- sf::st_bbox(my_area.sf, crs = epsg_rd) - my_bbox <- data.frame("xmin" = floor(my_bbox$xmin), "ymin" = floor(my_bbox$ymin), "xmax" = ceiling(my_bbox$xmax), "ymax" = ceiling(my_bbox$ymax)) - my_bbox_area.sf <- create_bbox_polygon(my_bbox) + my_bbox_area.sf <- create_bbox_polygon(name = name, bbox = my_bbox, resolution = resolution$res) } if (!exists("my_area.sf")) { stop("Too many or little parameters have been defined. Please add or remove them.") @@ -76,5 +87,6 @@ create_area <- function(X, Y, radius, bbox, polygon, LONLAT = FALSE, type) { #sf::st_write(obj = my_area.sf, dsn = "C:/ROutput/Utrecht/My_area.shp", append=FALSE) } } - return(list("area" = my_area.sf, "bbox_area" = my_bbox_area.sf, "bbox" = my_bbox)) + + return(list("area" = my_area.sf, "rectified_bbox" = my_bbox_area.sf$rectified_bbox)) } diff --git a/R/create_bbox_polygon.R b/R/create_bbox_polygon.R index 93db2f7..9574e3c 100644 --- a/R/create_bbox_polygon.R +++ b/R/create_bbox_polygon.R @@ -1,9 +1,20 @@ #'@inheritParams ahn_area #'@noRd -create_bbox_polygon <- function(bbox) { - polygon_list <- list(rbind(c(bbox$xmin, bbox$ymin), c(bbox$xmin, bbox$ymax), c(bbox$xmax, bbox$ymax), c(bbox$xmax, bbox$ymin), c(bbox$xmin, bbox$ymin))) +create_bbox_polygon <- function(name, bbox, resolution) { + #create BBOX on rectified grid by correcting the BBOX coordinates + new_bbox <- list( + "xmin" = get_rectified_coordinates(x = bbox$xmin, resolution = resolution, rounding = "down"), + "ymin" = get_rectified_coordinates(x = bbox$ymin, resolution = resolution, rounding = "down"), + "xmax" = get_rectified_coordinates(x = bbox$xmax, resolution = resolution, rounding = "up"), + "ymax" = get_rectified_coordinates(x = bbox$ymax, resolution = resolution, rounding = "up") + ) + + if (new_bbox$xmin != bbox$xmin || new_bbox$xmax != bbox$xmax || new_bbox$ymin != bbox$ymin || new_bbox$ymax != bbox$ymax) { + warning(paste0("Found bbox coordinates for '", name, "' did not align to the rectified grid of the AHN. The BBOX coordinates are therefore adjusted to align the rectified grid of the AHN to avoid resampling.")) + } + polygon_list <- list(rbind(c(new_bbox$xmin, new_bbox$ymin), c(new_bbox$xmin, new_bbox$ymax), c(new_bbox$xmax, new_bbox$ymax), c(new_bbox$xmax, new_bbox$ymin), c(new_bbox$xmin, new_bbox$ymin))) polygon_sfc <- sf::st_polygon(polygon_list) polygon_geom <- sf::st_sfc(polygon_sfc, crs = 28992) my_area.sf <- sf::st_sf(geometry = polygon_geom) - return(my_area.sf) -} \ No newline at end of file + return(list("rectified_area" = my_area.sf, "rectified_bbox" = new_bbox)) +} diff --git a/R/find_ahn_sheets.R b/R/find_ahn_sheets.R index 772c208..82b2ade 100644 --- a/R/find_ahn_sheets.R +++ b/R/find_ahn_sheets.R @@ -1,15 +1,14 @@ #'@inheritParams ahn_area #'@noRd find_ahn_sheets <- function(name, area, type = "", bladIndex) { - shape_area <- area$area - shape_area <- sf::st_transform(shape_area, sf::st_crs(bladIndex)) + sf::st_transform(area, sf::st_crs(bladIndex)) sf::st_agr(bladIndex) <- "constant" - sf::st_agr(shape_area) <- "constant" + sf::st_agr(area) <- "constant" if (type == "point" || type == "area") { #download raster sheets for area or point intersection - bladnrsIntersect.sf <- sf::st_intersection(bladIndex, sf::st_buffer(shape_area, 0)) + bladnrsIntersect.sf <- sf::st_intersection(bladIndex, sf::st_buffer(area, 0)) bladnrs <- bladnrsIntersect.sf$kaartbladNr if (length(bladnrs) == 0) { stop("No intersection found between the area and the AHN sheets.") @@ -18,7 +17,7 @@ find_ahn_sheets <- function(name, area, type = "", bladIndex) { #if (length(bladnrs) == 4) { #stop("The selected point is exactly on the intersect of 4 AHN sheets. Pease adjust the X and Y coordinates by at least 1 meter.") #} else if (length(bladnrs) == 2) { - #stop("The selected point is exactly on the intersect of 2 AHN sheets. Pease adjust the X OR Y coordinates by at least 1 meter. If changinig either coordinate doe not work, change both.") + #stop("The selected point is exactly on the intersect of 2 AHN sheets. Pease adjust the X OR Y coordinates by at least 1 meter. If changing either coordinate doe not work, change both.") #} } else if (type == "area") { geom_types <- sf::st_geometry_type(bladnrsIntersect.sf, by_geometry = TRUE) @@ -30,18 +29,18 @@ find_ahn_sheets <- function(name, area, type = "", bladIndex) { output <- bladnrs } else if (type == "pc") { # #download point clouds sheets - # bladnrsIntersect.sf <- sf::st_crop(bladIndex, sf::st_buffer(shape_area, 0)) + # bladnrsIntersect.sf <- sf::st_crop(bladIndex, sf::st_buffer(area, 0)) # bladnrs <- bladnrsIntersect.sf$bladnr # bboxes <- c() # for (f in bladnrs){ # bladnr <- bladnrsIntersect.sf$bladnr == bladnrs[f] # singlebladNr.sf <- bladnrsIntersect.sf[bladnr, ] # sf::st_agr(singlebladNr.sf) <- "constant" - # singlebladNr.sf <- sf::st_crop(singlebladNr.sf, sf::st_buffer(shape_area, 0)) + # singlebladNr.sf <- sf::st_crop(singlebladNr.sf, sf::st_buffer(area, 0)) # my_bbox <- sf::st_bbox(singlebladNr.sf) # bboxes <- cbind(bboxes, my_bbox) # } - # ahn_data <- download_pointCloud(name = name_trim, output.dir = output.dir, AHN = AHN, bladnrs = bladnrs, area = shape_area, radius = radius, bboxes = bboxes, gefilterd = gefilterd, sheets.location = sheets.location, sheets.keep = sheets.keep) + # ahn_data <- download_pointCloud(name = name_trim, output.dir = output.dir, AHN = AHN, bladnrs = bladnrs, area = area, radius = radius, bboxes = bboxes, gefilterd = gefilterd, sheets.location = sheets.location, sheets.keep = sheets.keep) # output <- ahn_data } return(output) diff --git a/R/get_rectified_coordinates.R b/R/get_rectified_coordinates.R new file mode 100644 index 0000000..c17f54d --- /dev/null +++ b/R/get_rectified_coordinates.R @@ -0,0 +1,30 @@ +#'@title get coordinates from rectified grid +#'@noRd +get_rectified_coordinates <- function(x, resolution, rounding = "nearest") { + + if (resolution == 0.5) { + # don't round if it is a integer + # don't round if decimal digit == 0.5 + # round down if first decimal digit is < 5 (0.5-) + # round up if first decimal digit is > 5 (0.5+) + has_decimal <- x != floor(x) + if (has_decimal == TRUE) { + rounded_value <- round(x, 1) # First round to 1 decimal digit + # Check if the rounded value is exactly 0.5 + if (abs(rounded_value - floor(rounded_value)) == 0.5) { + z <- rounded_value # If it's exactly 0.5, retain the rounded value + } else { + if (rounding == "nearest") { + z <- round(x) # Otherwise, round using the round() function + } else if (rounding == "up") { + z <- ceiling(x) + } else if (rounding == "down") { + z <- floor(x) + } + } + } else { + z <- x # Return x without rounding if it's an integer + } + } + return(z) +} diff --git a/R/get_rectified_grid.R b/R/get_rectified_grid_for_point.R similarity index 50% rename from R/get_rectified_grid.R rename to R/get_rectified_grid_for_point.R index ab1f507..a43bd0c 100644 --- a/R/get_rectified_grid.R +++ b/R/get_rectified_grid_for_point.R @@ -1,38 +1,23 @@ +#'@title get rectified grid for point #'@inheritParams ahn_area #'@noRd -get_rectified_grid <- function(name = "", X, Y, LONLAT = FALSE, resolution) { - #RD new coordinaten systeem +get_rectified_grid_for_point <- function(name = "", X, Y, LONLAT = FALSE, resolution) { + #RD new coordinate system my_point <- create_spatialpoint(X = X, Y = Y, LONLAT = LONLAT) coords <- sf::st_coordinates(my_point) - if (coords[1, "X"] < 482.06 || coords[1, "X"] > 284182.97) { + if (coords[1, "X"] < RD_min$X || coords[1, "X"] > RD_max$X) { stop("X coordinate out of range.") } - if (coords[1, "Y"] < 284182.97 || coords[1, "Y"] > 637049.52) { + if (coords[1, "Y"] < RD_min$Y || coords[1, "Y"] > RD_max$Y) { stop("Y coordinate out of range.") } - ## round number to whole number unless if decimal is 0.5 - #round down if decimal digit is lower than 5 (0.5-) - #round up if decimal digit greater than 5 (0.5+) - #don't round if decimal digit = 0.5 + #round to nearest x,y pixel bbox value to make middle pixel that aligns with rectified grid + xround <- get_rectified_coordinates(coords[1, "X"], resolution = resolution) + yround <- get_rectified_coordinates(coords[1, "Y"], resolution = resolution) - rounding <- function(x) { - rounded_value <- round(x, 1) # First round to 1 decimal digit - # Check if the rounded value is exactly 0.5 - if (abs(rounded_value - floor(rounded_value)) == 0.5) { - z <- rounded_value # If it's exactly 0.5, retain the rounded value - } else { - z <- round(x) # Otherwise, round using the round() function - } - return(z) - } - - #round to nearest x,y pixel bbox value to make middle pixel - xround <- rounding(coords[1, "X"]) - yround <- rounding(coords[1, "Y"]) - - ## create 8 pixel around middle pixel + ## get 8 rectified pixels around middle rectified pixel if (resolution == 0.5) { #x coordinate if (coords[1, "X"] - xround > 0) { @@ -63,12 +48,11 @@ get_rectified_grid <- function(name = "", X, Y, LONLAT = FALSE, resolution) { my_ymin <- coords[1, "Y"] - (1 * resolution) my_ymax <- coords[1, "Y"] + (2 * resolution) } + } else if (resolution == 5 || resolution == 100) { + stop("No support at the moment to use this resolutions.") } else { stop("No correct WCS resolution is provided. Please try again.") } - - bbox <- data.frame("xmin" = my_xmin, "xmax" = my_xmax, "ymin" = my_ymin, "ymax" = my_ymax) - return(list("name" = name, "point" = my_point, "bbox" = bbox)) + bbox <- list("xmin" = my_xmin[[1]], "xmax" = my_xmax[[1]], "ymin" = my_ymin[[1]], "ymax" = my_ymax[[1]]) + return(list("name" = name, "point" = my_point, "rectified_bbox" = bbox)) } - -get_rectified_grid(X = 150000.6, Y = 450000.3, resolution = 0.5) \ No newline at end of file diff --git a/R/merge_and_crop_ahn.R b/R/merge_and_crop_ahn.R index 1e7a2da..cfef9f1 100644 --- a/R/merge_and_crop_ahn.R +++ b/R/merge_and_crop_ahn.R @@ -5,7 +5,7 @@ merge_and_crop_ahn <- function(name, sheets, AHN, dem, resolution, area, output. cropped_sheets <- list() for (i in seq_along(sheets$filePath)) { raster <- terra::rast(sheets$filePath[i]) - cropped_sheets[[i]] <- terra::crop(raster, area$area) + cropped_sheets[[i]] <- terra::crop(raster, area) } # Merge raster merged_raster <- cropped_sheets[[1]] @@ -25,7 +25,7 @@ merge_and_crop_ahn <- function(name, sheets, AHN, dem, resolution, area, output. message(paste("Masked raster at ", ahn_masked_raster_filepath, " already exists and will be overwritten.", sep = "")) } # masked raster - ahn_masked_raster <- terra::mask(x = merged_raster, mask = area$area, filename = ahn_masked_raster_filepath, overwrite = TRUE) + ahn_masked_raster <- terra::mask(x = merged_raster, mask = area, filename = ahn_masked_raster_filepath, overwrite = TRUE) if (output.dir != tempdir()) { # Write the cropped raster to a new tif file From 64a3276d1fdc362c81b645c0bac808e4bfffd01e Mon Sep 17 00:00:00 2001 From: Jelle Stuurman Date: Tue, 16 Apr 2024 22:42:26 +0200 Subject: [PATCH 3/4] improved documentation and fixed spelling errors --- R/ahn_sheets_info.R | 4 ++-- R/check_ahn_version.R | 4 ++-- R/create_wcs_url.R | 8 +++----- R/download_ahn_sheet.R | 6 +++--- R/download_wcs_raster.R | 5 ++--- R/extract_elevation.R | 6 +++--- 6 files changed, 15 insertions(+), 18 deletions(-) diff --git a/R/ahn_sheets_info.R b/R/ahn_sheets_info.R index 91a418d..18b2596 100644 --- a/R/ahn_sheets_info.R +++ b/R/ahn_sheets_info.R @@ -5,7 +5,7 @@ #'@author Jelle Stuurman ahn_sheets_info <- function(AHN, dem, resolution, output.dir) { if (missing(AHN) == TRUE) { - stop("No value provided for the AHN prameter. Please enter a valid AHN parameter.") + stop("No value provided for the AHN parameter. Please enter a valid AHN parameter.") } check_ahn_version(AHN) @@ -15,7 +15,7 @@ ahn_sheets_info <- function(AHN, dem, resolution, output.dir) { } if (missing(dem) == TRUE) { - stop("No value provided for the dem prameter. Please enter a valid dem parameter.") + stop("No value provided for the dem parameter. Please enter a valid dem parameter.") } ahn_bi <- get_bladindex(AHN = AHN, dem = dem, resolution = my_resolution$res) diff --git a/R/check_ahn_version.R b/R/check_ahn_version.R index ca99c75..b5f3f2e 100644 --- a/R/check_ahn_version.R +++ b/R/check_ahn_version.R @@ -1,4 +1,4 @@ -#'@description check the version of AHN that is called and whehter it is available within this package. +#'@description check the version of AHN that is called and whether it is available within this package. #'@inheritParams ahn_area #'@noRd check_ahn_version <- function(AHN) { @@ -18,7 +18,7 @@ check_ahn_version <- function(AHN) { 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.") + stop(" The AHN version found is currently not supported at the moment with this package.") } else { AHN <- toupper(this_ahn) if (note == TRUE) { diff --git a/R/create_wcs_url.R b/R/create_wcs_url.R index bb20cdf..f34dae7 100644 --- a/R/create_wcs_url.R +++ b/R/create_wcs_url.R @@ -1,10 +1,9 @@ #'@inheritParams ahn_area #'@noRd -create_wcs_url <- function(type, bbox, AHN = "AHN", resolution = list(res = 0.5, res_name = "05m"), dem = "DSM", interpolate = TRUE, wcs) { +create_wcs_url <- function(type, bbox, AHN = "AHN", resolution = list(res = 0.5, res_name = "05m"), dem = "DSM", interpolate = TRUE, wcs, LONLAT = FALSE) { #get BBOX extent of buffer area my_bbox <- paste(bbox$xmin, bbox$ymin, bbox$xmax, bbox$ymax, sep = ",") bbox_url <- paste0("BBOX=", my_bbox) - #create image pixel dimensions if (type == "point") { my_width <- 3 @@ -13,7 +12,7 @@ create_wcs_url <- function(type, bbox, AHN = "AHN", resolution = list(res = 0.5, my_width <- bbox$xmax - bbox$xmin my_height <- bbox$ymax - bbox$ymin } - + #print(my_height) if (AHN %in% pdok_versions) { dimensions_url <- paste0("WIDTH=", toString(my_width), "&HEIGHT=", toString(my_height)) #name of layer @@ -44,6 +43,5 @@ create_wcs_url <- function(type, bbox, AHN = "AHN", resolution = list(res = 0.5, #generate URL wcsUrl <- paste(wcs, name_layer_url, bbox_url, crs_url, imgFormat_url, dimensions_url, sep = "&") - return(wcsUrl) -} \ No newline at end of file +} diff --git a/R/download_ahn_sheet.R b/R/download_ahn_sheet.R index aed8f1c..08c7565 100644 --- a/R/download_ahn_sheet.R +++ b/R/download_ahn_sheet.R @@ -29,16 +29,16 @@ download_ahn_sheets <- function(name, AHN, dem, url, output.dir, sheets.dir, she ahn_dem_file_path <- paste(ahn_dem_directory, ahn_dem_raster_filename, sep = "/") #check if sheet exists - dwnld <- FALSE + downloaded <- FALSE if (!file.exists(ahn_dem_file_path)) { print(paste0("Downloading ", AHN, " ", dem, " sheet ", ahn_dem_raster_filename, "...")) utils::download.file(url = url, destfile = ahn_dem_file_path, mode = "wb", quiet = FALSE) - dwnld <- TRUE + downloaded <- TRUE } else { message(paste("Corresponding dem sheet at", ahn_dem_file_path, "already exists and will be used.", sep = " ")) if (sheets.keep == FALSE) { warning(paste("Through the parameter `sheets.keep` is it provided that the sheet needs to be removed. The", AHN, dem, kaartNadNr, "sheet already existed on the filepath location and therefore it will not be removed. Only existing sheets that do not exist through the `sheets.dir` wil be removed when `sheets.keep` is set to 'FALSE'.", sep = " ")) } } - return(list("kaartBlad" = kaartNadNr, "filePath" = ahn_dem_file_path, "downloadedNow" = dwnld)) + return(list("kaartBlad" = kaartNadNr, "filePath" = ahn_dem_file_path, "downloadedNow" = downloaded)) } diff --git a/R/download_wcs_raster.R b/R/download_wcs_raster.R index 5f53681..ba1c516 100644 --- a/R/download_wcs_raster.R +++ b/R/download_wcs_raster.R @@ -1,6 +1,6 @@ #'@inheritParams ahn_area #'@noRd -download_wcs_raster <- function(wcsUrl, name = "elevation", area, AHN = "AHN", dem = "DSM", resolution, radius, interpolate, output.dir, type = "raster") { +download_wcs_raster <- function(wcsUrl, name = "elevation", AHN = "AHN", dem = "DSM", resolution, radius, interpolate, output.dir, type = "raster") { ahn_letter <- get_ahn_letter(AHN = AHN, dem = dem, resolution = resolution$res, interpolate = interpolate, method = "raster") #define radius @@ -12,7 +12,7 @@ download_wcs_raster <- function(wcsUrl, name = "elevation", area, AHN = "AHN", d overwriteText <- paste0("(", radius, "m)") } - #set or get outut directory + #set or get output directory if (missing(output.dir) == TRUE || output.dir == tempdir()) { output.dir <- tempdir() } else if (output.dir == default.output.dir) { @@ -38,7 +38,6 @@ download_wcs_raster <- function(wcsUrl, name = "elevation", area, AHN = "AHN", d utils::download.file(url = wcsUrl, destfile = image_name, mode = "wb", quiet = FALSE) print("Download raster image succeeded.") ahn_raster <- terra::rast(image_name) - #print(ahn_raster$extent) #ahn_raster <- terra::project(image_name, crs = CRS("+init:epsg:28992"), overwrite = TRUE) terra::NAflag(ahn_raster) <- -32768.0 diff --git a/R/extract_elevation.R b/R/extract_elevation.R index 2c6d819..9042b3d 100644 --- a/R/extract_elevation.R +++ b/R/extract_elevation.R @@ -1,7 +1,7 @@ -#'@title extracte AHN elevation (m) at specified point location. -#'@description Exctract the elevation value (in meters) from the AHN. +#'@title extract AHN elevation (m) at specified point location. +#'@description Extract the elevation value (in meters) from the AHN. #'@param ras Spatial raster on which the elevation will be extracted -#'@param point Spatial point from where the elebation will be extraced +#'@param point Spatial point from where the elevation will be extracted #'@inheritParams ahn_point #'@noRd extract_elevation <- function(ras, point, extract.method = "simple") { From cb5d215f5940e2ba9f0c2818948a95850d662e01 Mon Sep 17 00:00:00 2001 From: Jelle Stuurman Date: Tue, 16 Apr 2024 22:42:47 +0200 Subject: [PATCH 4/4] release of version 0.98 --- .Rbuildignore | 3 +++ .gitignore | 3 +-- DESCRIPTION | 2 +- R/sysdata.rda | Bin 27010 -> 27137 bytes README.Rmd | 6 +++--- README.html | 11 ++++++----- README.md | 9 +++++---- docs/index.html | 6 +++--- docs/pkgdown.yml | 2 +- docs/reference/ahn_area.html | 16 ++++++++-------- docs/reference/ahn_point.html | 10 +++++----- man/ahn_area.Rd | 8 ++++---- man/ahn_point.Rd | 6 +++--- 13 files changed, 43 insertions(+), 39 deletions(-) diff --git a/.Rbuildignore b/.Rbuildignore index 12b73ac..9b855fa 100644 --- a/.Rbuildignore +++ b/.Rbuildignore @@ -18,3 +18,6 @@ ^\.github$ ^Packaging$ ^\.Rbuildignore$ +^data-raw$ +^builds$ + diff --git a/.gitignore b/.gitignore index 3115703..1539de2 100644 --- a/.gitignore +++ b/.gitignore @@ -6,7 +6,6 @@ inst/doc tests/ scratch/ AHN_output/ -Rcode my_output/ README_cache Packaging/ @@ -20,4 +19,4 @@ deprecated/ docs .Rproj.user builds/ -.gitignore \ No newline at end of file +.gitignore diff --git a/DESCRIPTION b/DESCRIPTION index cf46b24..96283f2 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -2,7 +2,7 @@ Type: Package Package: rAHNextract Title: Extract elevation points or areas in the Netherlands from the AHN Version: 0.98 -Date: 11-04-2024 +Date: 16-04-2024 Author: Jelle Stuurman Maintainer: Jelle Stuurman Description: Extract elevation points or raster areas in the Netherlands diff --git a/R/sysdata.rda b/R/sysdata.rda index 68f176b4da678f2b9fcf52de180a8b46c5ab057f..4f8aa016c4bef1e7b602e4e7c4f2f93e6d3f4859 100644 GIT binary patch delta 8832 zcmYj$2{@GB_xOmELQKe#vJ5gI49ZrPlCi}wWGq=HGt*FH%QBQLA}8AWV=MJE!hMsWcCo{( zaew~8!OowJKRbVZv#RSJ*Bb6Sc`h3Up`mCD;`g6Pty?E@jnApW%1rTon)v_-I*9TQ!?1sG#V#zs&!J_SQOOIT?ydP2zsJsKe z3sBw@2e9o70`4Eey;+S6DG6A>2L=TO$231jqoBKJ5n&0w=a2Q5-T=H<{aMh;^IyrL z$1)%$=vM`MvMFr@L;WQaK)09Lf6VUx$sCd$GI^$=3l)wbiK=<^N_h0@d2uGHawh1p zpET3sj?9rIamOU8siKm3aKh#)BBE*>FDNN8x<(`vK@Y(CDBG+5>H1e;%h6Yd9?fYxm;5fdhTgAsL$Ua3-?L_5p-ro_^1c-Nl!b zgSv~(OZ0F*5>PAflI4n4laE%h>E|aaL?>57CzlFEsX-qrCmwB3f#(X5-Bd+iC%1_L z`)o~GE1pO-yv!y>xhaRYx-{%elvewGtr>Iso$J%WbI*Zk_b$UH6Z6RL1kC4y$ft3S zUpPVf&z9rZ?l@2*+nm&2x6iAmBhgHZ9r#j?Oh~uyR1$*8nqc`&9JAa)Z&H}I75n}ksS{bQ z#w?nl`tq`~PZ(l1lx{N1Sl@D_&%hgs;KQ7wRh>vph*~GczGy(q_lShc6{*5BB~d2k zuTdfPnbnKWeYtsSF3vRH*6FaNVi0q*N;O|2=?3#k#!R7lI+|JSH{$IvZ6^2Rx{~AQ z&qeig9I1?YTT9a_C43>Xaw^pw(M)^Z>H8$^tzr}XB1jUiylLJX}|5l-$L|CHsORl)6EgGwpNeVrh-h**SI zqIj1dDTqCoKnx&qy+%$N3M(%<`BTNHgY{%8vz>Bs81;~Efkb{-f2YvUUA*kkU_H(p zpiah{0mQ{>*h<~mWSb>N#PU>-68O}b<5a(Sa*uUpyXFnKNI8E0j^{>d zZJO351ZWR|VQi6_qwQ?d6&)r|G#Bu=I;RjNowDhFQX89>a=vT*TF=e#^th`6m>x?) z?r{&Jd*}L_`BNdnsls0$mEx)jf(4q-w-Nel*2($$T9#(CAVtHsFA8HSSA@m{#(0K5 z2|cv$FgWjOTu9MbZxxUPjsedZ^vwv%okp6^jRku<$!E=k&f5&f&Ry-fe{zIl)%{#M zEaz3~^mT)jxZLvuUyZgQu50|a94?~j$!aHi{KRI=XY!KG2UdmEJY#T*9sERl@0sUQ zA7V9mf}2S~Y19J#E$gpJ^e(-mvOg4 z-?Uk5uFF)k@{v)^VrvpEo7ts(U0lUtTx&CNL}PiR2O{=BW2H4ehHGA|ku$gWRSnOp z?ah2eqep6}SYStzCYQ9>wI)VwPgx^p_Woyo)YtYOq0+&2_oYZurbR}vpn7GU!xfpsx4hno80(Z&P@Rs=j2>HzM4VY`J$|#E8C$bS zAo@z15G30VaN;Yx5#aKCPQFQrV%QKE0uf<-(lis zV2a&F!nsgM>J2~#pv)rtYPsgE>ffSu`MEoxT|R!_e+P8MjAv16Ny`VZ%-=61V1Bbm z-(Slw2bZ=2T&H{DjRe>~7mja*&2U}tpPwCbn0zpDX`(nP?Xt%F*f|;LE+S}g35nHb zr_VmN-uMxZOTx*81qEM)-J5Ki{QRkLt^8KZ-aBkz$(7h-Q0@EtjBoEI24}_&ZisgM z2`rOKjY>UaIuL|uEPoP}jey8sEteG%iI<1)k(-hF&~)yGYkULQHiYup>AKH5ALKTF z{T}Eff|j?nOz`)A``ymHJ7hy(yc_VDclM)w_Qy=_;|}(H><|;(1L*M|!&{zEjX;$~ zi>zw*O2T4!tG`CRyCCSD52kPS{*unbiINdou#V3M3H+iZ{`=C(VW|K(Fqdd)Ix-jy ztb6SO^>hQ?O3E(9rAZ~L#kE=Tc)f1OkRUh7$nt}J4B8OO1mm-BLfx`uA}`!Y?uoYu zHAvCse98+?SMVe_^PI(Du3?3!R~3q&$fF2qnT@GE_Q}mXM&>JLVy=_aO=^LxCDhF- zzWcl)SVX|I^|dNW@_BN*l|l++jLloF%aWeaZp&ui3TvG9I)0_JNzu9PVVrD>T-P z%cK=`s^>ve5|F0Y?h6y9V1s#2n}koV1MXVsRwW#IK(l))tycuw?=y>k@$gj6O5L|D zCS%)d`?B@`uKLBn;?rAiaAKGX^9*P zD@TWEmiyJrkG?bA705!{!Z(DP-wc9(MWk8}>eW$jK<6zkr9Gf?*V_1pJ?n#o1eq8& zhGzm&6@t2Qkv!&b3|Ely~K)nkTx%+A6(QbVH<0Uv~Z}s|y$&GMaH?9 z(hvGQq^Big;xzWqlpK{-#--FL2oF97y%Uwj@{@*XjP=>C$b3$c35==NR;K0y{ z)!miyhBaRzX#8>*qxY9DIkUMd@Z4y1x{+!=xm|7791e$7)ezxu@D`P)!tju$~~!)|e6jru>DUp$rE6Kj5{naCQ#VpWW&RJ-|{i%RhLIo>bvTkTlS6 z1A35v&-689+r2or4>7o#6LL&I$*sTy6_Bb(oDqyK`Em5okLQd|LZE8K$_eHitQzKkrdzu#Z|hyRAoBE{RtzW}Bb|1I1& z!T%T$@-kS^2$nvA0^y+=|L6rJd-blg3Fa`f_Cx_qq;Lw30RRC2HkrawC#7Nn)yZvM ztFDvDAY2m0;O4iyiEw`=b8bvdsKy=3UQ>@`&+Xn(T6{lwItODy&6K0FIiZnrOA97mh4+zOzNgdI8=0DOL(dSm&tn^+HOsOJ* z;c%;J7=@CD9R*Cn0j#mIPzrDaCrV+BEbu=F_@7gz|8>plw;b1OMeg5+qaH5RCM-AC z-xKCZzxH|XJ1hJ6%O~TZfqxZcOk&yC4jBU?BDCT6;EaE3tj7(9JJDHoS==2(h_Fx8p-DtCC;R1rM5}56Bm0q+Zkz4`8+(cNmO`L|Pd@f=$m@NKew&5pM`#t(aXM z#tB959j6sn45Lp%S z^ePlY*}kTT@Nhn|bU`OtQEUZq*%gAxuAG+bqAYw(`AbjgJ%0jza;TnmI1*)TLNaJY zL}v>bGzcAjh^Ba2>w78cyCCrWP|x#VVr1kSWpy4y8&y%G0uk;vZYrYB;(F8)jX$A8 z$!_T;ZqM`L-AJ;JUn`T0p1^W#-j3Sd$q+813HR$e94^AE=+Ewm4d9uiVm`Pu!lRAZv07%PYg|Vb@`HrE#=Q{6U4?f34pBNROR7`CUXZ+B zv5p&h1!J&T#RLK)-9a#zG*UVEO&ZLMF{jyem)l3>U8_h0(u1PQTUC-dcl-8j%`OH5 z&og^5on4tbm}p~KAesI*{Gs8B`f}9Ibw05tG0P!tu*xKM#=(6KNj7$>UuEX|KDxnkqz@J9OOf%``-_txXvq2Vkkc&7IRy{QHBu^M{f#R@%<289+Eam^d34_@wbe zUFEkme2`xRz%L@`74pQTm9J6knbUKc^)9C;Mc9>tTRym#Jo(!RS94S>8T<((nF030 z_$RoNm?{arOtcj#1=&>8iMC&Bhn!(!FuHu9%F5jJZLHBxr{9_8pVaj|#1|0}kqqp8 z+X#Uzx(n9|cL?E@k$6sx z-Ric$aB=_QKlY0cwu@=h;P8&psI)g+N)5nU;z~c~K3j3IvGeh5N|;>g;8u3#b{D1< zn0zdBmZfC~(^xkeq|HXA>-1zjUhfj8i$`SQdkd%?PWOg$cx|0+K2Do;u1756*+!3q za~cbe%xOX%-^=^ir7;PFw#oWJGr3JQ>)$J23jn;vBU7Oa?qrq9E>1^mr34#mo8Zkf zW$y-Kq(jTg5kLH?n6LiKt;L$nt*tZObx#okxTj5D|ANFs+|e`25OxrEJp}LB2uf$* zB94F}T>ehaMT9X2p6gsJu32=#P(xy9V*`k;m9DO?mGOC=?MU`2?)uho_ux!|;^f-o zE+-TaB=YEQL#(V_E)_G=jt@RH^#KvLb53bsSU~k{CYT0+Dqq391_Ko*8yL603 zb$v%x76t~md&!vR;^4`0yxAnq&0cI~8-rLJ@GnjpZ-f8s505-* zOy|aX%z^MnAW2#}mGhec^L`Uy{pfYrlu_H(^HA~O$7tG0XFpR9zFf2>4T7M@+a zeqgy(0AH0C6#?N{*10SBH@*mH1J87@j9dyhWX4v!9;j4a_ANA=l2N!ibu zzopDa^P-^|LpaHhINz(#s}^k!98V;2DlJ`B|Be6i7J%)V5HpZB$(&F1D zWAi|!Q@_tQg-g;`A_?_fvO%|XYV?uxL&?*=5sJzDIN_5Pj_z@`+<3I(V0)C5()FX2 zN3Th!^XPf@sLKzDWnVYv>$PAzZFQB5|7Z$&J31zFQE93>%KmPDk9c1I8VwboF_MB| z7FL}WM3u6d6&=4(jeCLKTk{`fMrx9w0etwyA@@%&Y}Nw}{e|f3{q&C#01l!f8fnAe z-;t5SMQc|Mg-sX8JqbJ*`{{&~l<5kn^h;(5*$TAYjP!r!OU^pnX+Z{6zB7JzmJyZqaJG2e zv;W5nRiLwWEf~5s1J(YlGOm8+@bNKii*FmrIQ+($DN0n@`EPvqx!`VS(68~5{)Mfj z@x0|*Z?flYlkHNm_-g`@ezs~;Noboy{RZ_jgb=AT_PbG%zSm(CuheO{qQ>M7C{XO2 z*&`4L1_gOUB}wq3@cgZlrg(Xe=s|fu+hm($bVow+qoXlAT;!;09C2~RkLgQLbM4>lDq_Q;h>ny0sd>^tp!=_?4AByVTZmjAi{rU~+h*s9I-_#y#&tf^)i=cQ=+ z4CDNT8Y-Jq!z^(wL<}>N%*Tbg9MDx$?Md8uD#PEl-}VMhwTYC_@K+`_3#?a;tB9PG zzWk%{cuOe8HlfTFvviOn^n7Y+X}8P5;W8w>Tp`pBl3q(lUw42*byUur2Fc+jbU)y? zUMO}__m5D|BozfdyW$zW872ogg4T_xci;lsO|h{ZJv+t*40~#rzWJwj)2E&SoqSrW z8Ln0!huwH@4D&O=huLqwlpmVA`SZmI%{5*z?Gp6LVYy2%*XGjY%ea}#$180JX<6aD zys?eTJ|Gw?7xZabRQAtjnVmF$JltVYGFfY|I+jQn8g&P$+55+3>{rsYoQO2 z1u_F=pqD@PWC<+X7wcTio!R&byK(yW#mK?CV|en53}j(o<_;WoBXFkh-t{9Fj0Lh* zi#NtI!h4o1+CHE0$s3f%U0H&!%)yy0Rm$IC>=V(lbFD7hv*dO z(8i`EujlN$>-Ao%7a$!d#zIJ}Yx(W4?KJzYE!M9^CCrJ-^T;~$y`pz5xAuNhboyyu>C|A+JO6Q3*ZzcC7FsX8V zXFu_<2mVNJES-dFC#P#K-4KLXptr;{56h*OHIBmpy|+&U%q>~mZrY&ldbeaW`cXf3 zR0zDFg};R7%Vv;nX)n9iIz3|szd!rc1pnX)C2v24{mLoLU)Pg+sXHa?OPabRC95gy z9Xiy&qWwQ3ql-JC0~JHzoIT_h--9M_Yv#XFHM8zFudgrR#{hW z3yoEi_9{)LFwk*18NY1|4Hnl4wn!b`Y*9Q&hgq=jZG_+VmRAoF_dKCmzg`m7d|B^@ zPI83#-7MZp@IdkG==XOIdJhqvx@qX1|G2%Q=vbqmrC+u5?|}aJC-{${*7bKcQvwq?`<@&Uyt6I28n{rJPAH`v~E1 zDP>&`H$rEFN-i%+%5-xn>B|h{^cUEAB9SWT$+)f;McU1}^legb&WFAgosyZ{TYRZ9 z=?RR+v#GyIYD^AN&OUhdM>Xx0{vWBrsKQ7&4y9vBf(RxY1AkppOR4qtY?G9Pu*(`M zw;Ur6+w9cRu$B1%3mFT`wWwio@AqEw+`3*=9x-o>p2125SveuS7fWX4gfDu{Y1lCw zykgjUVOfL43mCai88O%GuUJ;6`cu6NGzc||RE$Gy#fXm+rnu53!zTn)(TmjuQ+p~> zE1?Lg+On>X1M7vsffhg^1-ux>YN0B2YHe-~_vDb6g%?y>5d97iWeTOcy3=oWKROTLRTS%s<;#c$$@C%*T~Qo}iL` z^da9d)u;p?2o?P+@>CZ1vUR8RiP5$F%XWS7FPqNoHBk_x{$Gn2KOB9;(mVnH;J_P; z*1Y($p+T8gMCA8A3}4NFy>}lvI3~Sw?m}zYJ7~JQYS+5q#wVw@Ik*F>vj--7jo95n z#MC_rAjhdg#!|+dEhrjgiAdk1(fB(Wd^eTc< zhvz#A;5eiQkT?~H3-k~4B*zS+zN?K5L*qDk9w}A~vb(zKpq*{-; z23DDG82S8Wheo{uwp93I#2~lG-(Pj(l>txQJ0#xP>vt9%LFloxGRmcX5w9(VlIwQ2W_8cwF1_W(apP1Zah~N{Oo0C~U`4 zOUnC39P+62O(uQR4!|aN`rgN;7S^+k1OQs8(eE!Z)47lB)-6O4izFOLqhPxqLSLnxd_5AL06N+CssXu{9v;6p~?SHig> zZ~d0Ho9Ek6?XQZMO1kt9d6oW^5r6~+!6Q~GCpG>l>cT)CUM|=q)a05h{%IoXX(Vjp n+x?fv&qgE79-7YIUkfA+L0cn6grauaP@}Xcj~%173Q;SuRn?|y zh9*X9mzq_IdvE{g{eJ)7*YCVeKKVT7oX_Wc*7M{%??g?}J)NYhbvL{uFphf3V4e5t z-}d7;RRE0DTD(*S`&!5dNAlzSK-WPa&@LljE8U~${$ng&2n0HZn>*`bXiEm2_cu}K zdZymy(WW^XlX0O71Y+2Q^V-S{sM%;N^n(}&PBlUN;Wv(CK)bbcgK%5n{oA9*{!^Ku zZ413d7J`!Cu{cAiG)8)ff?Qb;n_(nvCtAaT=n7T;-(+FK5lkV8S814E$8*t1JM~Lp&dJ;7fJ0zr9swN>4nB(;*y_o#472e|3c7=Modv0Iq?)dq8d^HP)%+yn=|8P{8a}A z@0R1ps0Oo1M7KJG`HWVBq)8QXG-70AG7LO0C>jg#mGTfAA=eR95Fh*@VFHjrNJVs9 zzzEqDfvUsJToy7O6fHEw{?(k9i&Fi7R>OHnXhN7zVlcxm9EQ!NaU&nlf-6Rh_eBY} zA`(3WaO6hX06I%n$Z7~Ivz&xhBcQw|%|_Li;wGwPzjIayxOY^kRX5l{++41Aqds*i zlj}b?5C&;Mz8AQhsfDU0qjIWe#6Bxt8>-P5S$3$hBhakHPG$x~$17j@2KQOb$-|)< zR2)fU%IZZ1hBrp?KS%A*$|KX&;ffL^a4GW)vbYBY6oi+vCGo z8aKi}pitFlhtXtPcTr9YV{!4)xrYBA_>;9H`Ci8<9D)Nl*swnMM>!r+d?MZFr6I}*>t6ig}yWmimAU;{1w!zzM) zf&u1_m^xT>jyriy5s5SMWt*4H1i! zmvTxncog|h-A3s76x$ziDQvWI^=ZnXw$Ef4nttzFqNf!xwotbY9CA0MdF}U}_rrPE;G{Gt8nY@R~ z)u!1kpUZ0XjdiUr?|5gBXT0I!=|3tv{KBc`i0Iy@eMkIld|f^lpDbW4U2lCY`ytS& zSUyzA6F+HA%Fnn)>9tL}rN?;JO(_{e2uvbRi@g6>t+N|s;ZHf1t-5uxy_wKCS+6%- zSbk&tURpsJID_+Qq*j9dbp1Ut~nxCwSN*Bxg?1>`rjvWcupE-pTCd!%!WDx;Fi|tQ|fiFOr?2 zFXG~0s-i8fIX*4I*&`&*JX}ZyXs^PDQ%(?X>?`JZUONyi9*_@PAG%J)pvS!KGiufgxZfCj)~r~ z2RD`3Z>R|HMQDj!ukA6etwn@S{198Ni+#L3Z`1s}eLdvDG`4zlzuoX_UP<%zuz`yI z@<;Ar`8Cvd+_!7InYR=im4uX1SL1B1D*;)_s)TGg+1hK(LPx$`KFt;bX|W!`dz)N$ zOumrv3wKI=`BsNjIPc-(tpz68)A}C-=24N~}pA za{4&DGn?-+F9czc%y@U`YeNiB$R|RtI5)U>uU&OTBxE)iIy%)TOp316`SO5I{(j=6 zpiK>HFT;oUTlr@eL-#R3&xtqA7zETm`;b@@@5i`m*v<~w3zUlWA}0Sxit?0UXD6fj z$|THQf35G@JAT%=Iknz8EK9im^w_qluT_ps+DdIqjl6#Iqna*OlNHp0`yFuAIZS{HDRy^EDG{bl8tt~*7A;r6{=9uBFddp zUa7BH&(uZ#?I~U|XIIthRP@pppAG643iW->?)(|nDAih_9OZE}intdc?Au-Qg$Ke) z_w`a8(j+E{2%bH++-t`FLn-DGH7GQL|G%euO-=$s4Knm&uT@k?&* zaC>LhD3cQW{p4CC@J;9AE@tU4LOra&($|un%bv?J?5jxT%tfe;H%c6qU0>h8V?$#S zN1moBtH>mn5ZA_lZ((=aeippUpFLBBw(0mj$8tG|2 zu(I|rcP`1Ot7u6?2G@Ej6i#Vac!^0%zp##x-H?1v0o)lgCHV>T^%F&@`fZ0&c7Q+v zTL#zQEuEiyOPvd{K27@fY`9r7`LX9-R@j`W>r?b822!lbu~OFjCpIPeNdeyRQ#|>& z7_>P5jANx#+yq;9S$kPAYj}h#Zpf6aWHuV$!DR&G!dK+rVcKX^A?=81Y|IhYOrZi* zV51f0dg$I~?jyQmu6Efw--nrvpa2XC7|+>A<$9Ip_Qu*j{kx}L(&7=8vAR3HGvB03 zdH~;MqsrQQH8dgh_M%mXsOyyT(5F{yAES13s>l&nJtLeNjXnj7x(7YF*nviT(Ajz& zVY&A)TKkvX-VMdmOWZv>*NvoK2+uv(W7SUX^eqwL&grDwg`X5q38BJGSil=)SfnR%*7Bi?YxbVtBY`>1_IjGy2nqBI{6WW)#f+^3YUAld$>$qxBZR%i zp5Y5WeKT{Lyp&3eEX&4S(^hpLlYxg!?ss#)JFQoh@6Wrx=gz$J<%Vnl@(UdL_e}d= z^=TAVi`vpn&&$_$oz}M(g3Dk!VK!#TeivFNhrO2ew2-w2XMYa;($ZBw@&=l$lV*m4nO6eJy^JbOOB}6VKTWCilEMx;;)g6;L01GpQc|tp{%O2triyJ-&M=Rd+#aXq;~P{RHup5k2LUm>4q|Ll!;7 z92~6pW{iuclzu$G63v{Vn`R8Sbup(%0o+PBDe?3HrMMKRvJH%d84E4Fuh<1O8{{@j zFmy2+gapd_37yYd_Y{~e&VMW&mXt6bxFNhT1(T#`t{S>{DU?g9s9d|=;au~t;Uvoc zHyH|iqrDx`JfJywC;0z*0@Oiie$&>afPN3RC(aFMYk|G65y)0pnN(T{F zz9k{yV%9TqnNgV)<;unu`IPa!BESSQcHuWsVH@PQ(q{lM;Zb28slUvrU!r1$y9K?M zz<1wNwery5$VJX6-Y)v5gp@wO14XVcR=E6M{x~Z?>7T>!*3zeX6 zWf}fOrFlVXbY$f3l>P3JG>!tIqw$^~o^*N8e})@D0LVYz?5o1lFTo1s9)I3j2~z(= zeo?=m%2dMp(B<3Gm^MZRFy4>>O#J@D>N4P{x-9>Rmo8`v(@5LlUp+nmOp5sXluxEn zNJD)XJkCi+4RNo{VoN(a1%UZNCtxgVl6m{Sa#kD+Xr(CZHGYw=eAA1n$|66N1}8j zhzlmyrb+XJM1kq1gIJ_!DvL-|-<|FMI7y?l8umQ*zbnV>;DS2!^$Vj42IEdfLoIub zTS};QdtfCm?>}+-AEgxu2&NVJ^)MtvR~<_YrA^RBk*HAXpenAOrt)^D5Oh=1(9Z?- zkmDZ&6bao>v-76&;D6$@)wH|hPDn_I7tLR%`T6G3qNm`fo1JYtdBK>~t{# zx-=;++HNpOOyTcM=%AQd+SmD5AN92rcd<$#?QszMMioi0ri3#`m*+!6riG!-K z^b?IhBLI^#>B*$-_xPh_W| z6SSpx0d}Tnt9#x_*7(E{EF8f)pcpOV&Q9xA)4wT8q1ojqyK_KzajiBJ38nUGikSQ60QpRJh7$5u~C{c=SsIkfvV-g>w*yBYiBg$72&7#BsYc>Mj zGhRP6aoD?@OUgKY(7W=5q9?N|hM?yzA;Tuf4lm-$_8iwQ!%xHrNSHERz*Lsv(gdLt zBG-Lv21t-pl#}A)B10mwsN7SRkSD-6xqR*o!e0Yssjy0Qf}W>L;RG0C6Bqx^73^KC z*LPWo;KRb>b-`}1^mz_H?T18(sfUSD$#1Xfiu7UhtT#ZqC>gNuz-?h^oM14fR*+Rm zkm`Ew1tS=K6U?Ivw~dP*qpgVSG5#m3o4iOV^x4df&2a}M7C6Y*~Pm~T( z+&t6nZT|I3CG%r<;pd-B-PGmZb=DsLdiHtm=DbelW1WItp5U~4A|cUjC9V5JD!?eI zT}3vp%`ba2P-`%6T?^XdLMF?5PU&OVz0+UGF!PJ@b6AJGACIMv$T1HZH{c1z9uF># z2Ii&y)Pjz}l$zq)Z>QJ5az+CIL)xQ{Op<`u#gmI6o@m?Bv5t_8b6wZl?}Pf)Y67eG zwMIkpMh2m5U2BHl!J#Pl?JOI!k#((6|E8#HfzzdRq>*>|E%wU*5{PR$q_4p%SNHqajfSEpuZdOpj5Y8OuS@%cGCF3RiHq*hzbv z2Ur~KsW2#OWcTNVSt|A|s;FoULi_fcvKNdOT>Z499KW}2ZTE^q!^I8Fh6b^*KMlHbafub=Fv0DHG40!&(#u|Z$9+J?PDzPGitwKP$om}B*`{r1#$ zhP;(g$#>_qx*)!=-=(3FwB#_V@yMP&m9z+Dh#KiTFx0DcD!{eCn+Ub9@xNa94G_daGA?YZXHBJla=djwUFEV)B3vD`o#2RX;GMSWVH= z=)F-yD2^lYld4MT){JRIpyz%x5Gkff528u(i|fj^z+H z)D5>V=-?;}>X*$=c0Z7GzNAHN-690l)a3Su=jy=KcdqLK2MY_KE=SMO0z37w`6%qe z29;lu36lFjZ_WiMU;=|CGlK7F?PxAN_zVJp(%42uoExyUZ*n^j9XfB_*IUP6;lk?f z>yNp(>|?T?xT*l{=WBx8tH{(2fXi77*t#isW=edev~{XRIzojz%{M$j25c-RTLC-j zn{xBBrm@XqQ@$mENPx`*nL3muhiT#EMLjx0O-#Gl|kV+?7_;N^1 zSc@BJT*+e6Pk&lk8UiaVy-M&*Kh1#GJ63_gP&k*r*M^Ki3!U+8xA4?@g3Yn_a4lz~qVYgT8~UCcaMZ%~iMO!dn-!!01=tfR?$F z@9S@Mn%3oB$7D$cDpWgQkfZZU&T!)U=I?$vPfZDoSwR;NYlE$FA>iN-whf_=#?F0DFGevU$1s9%#!)^GStCf}c88naepHA2MUb<+Pd} z({_{sY=rE*ALzJH#Bx=vCIuFnSyO3aUg00))U-%!3_a5q(>H5~>Q7YQG{QE&YhHxrSrX*-$LKB9(=EiO&+5dIw2;__OiGRBK}q!G#NeNu3JBj*Wy41~kS&{vJyRtC zZpe+GxIc0tG6@1M-7`Hu%jK2&`#S#fnB z1G?QAGR$h>Wn=mA75wZNo$oroANKL-gHH<)$s0C-u~I_li`1`=HAp*f=<532gZmRd zJBE$8U)t}UE3D_kzLd5Canf^3gw9;tdlLIno%)XRr&V!GR;p`dl;nI(M{<92L?P@+ z7j(@)x?4Zg>G)Q+x+f2A>DpV~pAo_r_7}OdMT~A->rVS%(YpTRHW^P>SJz=UPbVrKRuytA+;sb~ zk^dYXcAy?)cp5yPuN`QAY-C#YQhG70^O^gW2@r<6mVXmAj|tt%2yusLchvA7BYx1Q zC4Z}w-kFC@S~S7h7g`G}9?8PiEe8+sFM{iB)5t#4zK71Ic-Zb=&*!Y)4Nsku#(p|n zXm~=|$S*0-!CrExwZX1bEsnO;y}bUu()iL-p4?XowrR`L(A`t*mb3{PqGmBow6U0} zD2@!gdUN>(*y8m2qa>|-Kb^5p4P7T1_eX!iVcJb<&ng}^ICuUX?2K;s!FqZVTYI}( zkBlL}glK<%PQqN|iIUH$-VO;?G$02A;-ZtHW)K5|u!!@^I+ujU!YvJyE2X~%^Dy|i zXc|a|O`JO0?GJ8hwQd4lZ=Aokcf3M2O^!Yk(J=M6zpItDQvKw{|B6({KToY^r-4E zdDn7@#JvRnt!ZND{GZON6$){8p5+s&oQ;BpJAyb88lI`YE}cI*)<=Ahi?Ry8NtNHMDhK-r?yPJ3H0$3DEj|`C4pV-45 zT1?NqdFfLrJ+<*jYIq=|KHH$cWNU zue;UJ`&>6ZnNEoP4OGQcMM*O%uqN}=up<2uy)+zWpfNG??78fdzJL=IMdAUC2gA?b zx$h!x!#S?U*wo4>E~jfRW(P2W#OP?S3y8B(j5l`LKI+^wg^_i_AkwCa5%PdqM!{x6 zSK3HlRk+M9w6#_7xdI#D*(hd4FSu_;DY4aL+S8x}N$#4T&CwfdGJWH_{i!K<%a?+b zBO*;(vWYwME!p|Gt9Q z1j(<_T=is>D}RFuyUF)%9b+1Vh;4LUko#vNbffAq;fU%$qHaAtnk%9mP@sVNzc9<1 zWVfTMDZdTMm1ATiV|v@W(+)o~M)3Zj23O}wWq!)RsvMK(XG*|yOe2~el4pO>zjBeq za+Lhz!k?IWG~}#}c`eoT8#TCsLb*NO9W8w8ZSpu?I?LwmTm8q_p$i+sW3U`o<&P;= z|Lnb}uEb)h?)*h;TCc)~(cp!8_aY^JZ7+YFejsAy<>?Biz>_~HWKy=brA^&&1{WN4l3}ipF@N8c0aGalk z&%RF}^Il9&df!c!nvXP?ZWaIV_FOuS>xb8+?Wv0K7!2LnCM^A>Mby(jtqit@%#E); zeqOQtPEqcd_3@I3@e2tIX@!QuB=yoz>)$!_DnBp>M7*qr4UBqB|2dq{|9dn2B!ksl tQdAN=l>z++x%5M1;cmJcm-9zU^R=+Q?K3DryTASjFwnG0k?#Hi|9_AInw9_n diff --git a/README.Rmd b/README.Rmd index 6dc1f44..c2a0b75 100644 --- a/README.Rmd +++ b/README.Rmd @@ -25,13 +25,13 @@ knitr::opts_chunk$set( ) ``` -# rAHNextract 0.98 dev +# rAHNextract [![Lifecycle: maturing](https://img.shields.io/badge/lifecycle-maturing-blue.svg)](https://www.tidyverse.org/lifecycle/#maturing)![GitHub release (latest by date including pre-releases)](https://img.shields.io/github/v/release/Jellest/rAHNextract?include_prereleases) -Updated: 11-04-2024 +Updated: 16-04-2024 This R package automatically extracts elevation points or areas of the Netherlands from the Actueel Hoogtebestand Nederland (AHN) datasets collected by [AHN](https://www.ahn.nl/). Individual elevations, and elevation raster areas can be extracted from the AHN4 using the geo and atom web services that are made available by [PDOK](https://www.pdok.nl/introductie/-/article/actueel-hoogtebestand-nederland-ahn). PDOK only makes the most recent AHN available. Currently this is version AHN4. The next release will also support the older versions of AHN which are made available elsewhere. @@ -66,7 +66,7 @@ For both the points and raster areas, the WCS method is set to default. The shee ## 3. Method of elevation extraction -Extraction of the elevation is done based on the raster data it receives from the sources. Determining the elevation can be tricky if no correct resampling is applied. The script always ensures a correct resampling takes place by adjusting the raster cells that need to be downloaded whenever necessary. This is called a 'rectified grid.' To avoid issues, the extraction calculation is therefore always done in the script and not through a OGC method online. Please refer to [this](https://geoforum.nl/t/wcs-van-ahn3-werkt-een-bounding-box-als-clip/6122/4) why it is important to always use rectified grids. +Extraction of the elevation is done based on the raster data it receives from the sources. Determining the elevation can be tricky if no correct resampling is applied. The script always ensures a correct resampling takes place by adjusting the raster cells that need to be downloaded whenever necessary. This is called a 'rectified grid.' To avoid issues, the interpolation on the rectified grid raster to determine the elevation is always done in the script and not through a OGC method online. Please refer to [this](https://geoforum.nl/t/wcs-van-ahn3-werkt-een-bounding-box-als-clip/6122/4) why it is important to always use rectified grids. Please refer to [this](https://gisgeography.com/raster-resampling/) to read between the different extraction methods allowed: 'simple' (nearest) or 'bilinear'. In this package the default for the `extraction.method` parameter is set to 'simple'. diff --git a/README.html b/README.html index adf65a7..4e79ccb 100644 --- a/README.html +++ b/README.html @@ -601,13 +601,13 @@ -

rAHNextract 0.98 dev

+

rAHNextract

-

Lifecycle: maturingGitHub release (latest by date including pre-releases)

+

Lifecycle: maturingGitHub release (latest by date including pre-releases)

-

Updated: 11-04-2024

+

Updated: 16-04-2024

This R package automatically extracts elevation points or areas of the Netherlands from the Actueel Hoogtebestand Nederland (AHN) datasets collected by AHN. Individual @@ -690,8 +690,9 @@

3. Method of elevation correct resampling is applied. The script always ensures a correct resampling takes place by adjusting the raster cells that need to be downloaded whenever necessary. This is called a ‘rectified grid.’ To -avoid issues, the extraction calculation is therefore always done in the -script and not through a OGC method online. Please refer to this +avoid issues, the interpolation on the rectified grid raster to +determine the elevation is always done in the script and not through a +OGC method online. Please refer to this why it is important to always use rectified grids.

Please refer to this to read between the different extraction methods allowed: ‘simple’ (nearest) or diff --git a/README.md b/README.md index a68d871..bcdd62f 100644 --- a/README.md +++ b/README.md @@ -1,5 +1,5 @@ -# rAHNextract 0.98 dev +# rAHNextract @@ -9,7 +9,7 @@ release (latest by date including pre-releases)](https://img.shields.io/github/v/release/Jellest/rAHNextract?include_prereleases) -Updated: 11-04-2024 +Updated: 16-04-2024 This R package automatically extracts elevation points or areas of the Netherlands from the Actueel Hoogtebestand Nederland (AHN) datasets @@ -74,8 +74,9 @@ from the sources. Determining the elevation can be tricky if no correct resampling is applied. The script always ensures a correct resampling takes place by adjusting the raster cells that need to be downloaded whenever necessary. This is called a ‘rectified grid.’ To avoid issues, -the extraction calculation is therefore always done in the script and -not through a OGC method online. Please refer to +the interpolation on the rectified grid raster to determine the +elevation is always done in the script and not through a OGC method +online. Please refer to [this](https://geoforum.nl/t/wcs-van-ahn3-werkt-een-bounding-box-als-clip/6122/4) why it is important to always use rectified grids. diff --git a/docs/index.html b/docs/index.html index 6973510..4b56705 100644 --- a/docs/index.html +++ b/docs/index.html @@ -65,11 +65,11 @@

-

rAHNextract 0.98 dev +

rAHNextract

-

Updated: 11-04-2024

+

Updated: 16-04-2024

This R package automatically extracts elevation points or areas of the Netherlands from the Actueel Hoogtebestand Nederland (AHN) datasets collected by AHN. Individual elevations, and elevation raster areas can be extracted from the AHN4 using the geo and atom web services that are made available by PDOK. PDOK only makes the most recent AHN available. Currently this is version AHN4. The next release will also support the older versions of AHN which are made available elsewhere.

@@ -133,7 +133,7 @@

2. WCS vs. Sheets method

3. Method of elevation extraction

-

Extraction of the elevation is done based on the raster data it receives from the sources. Determining the elevation can be tricky if no correct resampling is applied. The script always ensures a correct resampling takes place by adjusting the raster cells that need to be downloaded whenever necessary. This is called a ‘rectified grid.’ To avoid issues, the extraction calculation is therefore always done in the script and not through a OGC method online. Please refer to this why it is important to always use rectified grids.

+

Extraction of the elevation is done based on the raster data it receives from the sources. Determining the elevation can be tricky if no correct resampling is applied. The script always ensures a correct resampling takes place by adjusting the raster cells that need to be downloaded whenever necessary. This is called a ‘rectified grid.’ To avoid issues, the interpolation on the rectified grid raster to determine the elevation is always done in the script and not through a OGC method online. Please refer to this why it is important to always use rectified grids.

Please refer to this to read between the different extraction methods allowed: ‘simple’ (nearest) or ‘bilinear’. In this package the default for the extraction.method parameter is set to ‘simple’.

diff --git a/docs/pkgdown.yml b/docs/pkgdown.yml index 5f8a840..2b7ea28 100644 --- a/docs/pkgdown.yml +++ b/docs/pkgdown.yml @@ -2,5 +2,5 @@ pandoc: 3.1.13 pkgdown: 2.0.7 pkgdown_sha: ~ articles: {} -last_built: 2024-04-11T08:07Z +last_built: 2024-04-16T20:35Z diff --git a/docs/reference/ahn_area.html b/docs/reference/ahn_area.html index de0a022..0870b54 100644 --- a/docs/reference/ahn_area.html +++ b/docs/reference/ahn_area.html @@ -1,13 +1,13 @@ AHN raster area — ahn_area • rAHNextractAHN raster area — ahn_area • rAHNextract @@ -55,10 +55,10 @@

Get AHN raster area of a specified area.

-

AHN Area is retrieved through a circle (x, y and radius), BBOX (x, y and radius, OR bbox coordinates), or own geometry (polygon).

+

AHN Area is retrieved through a circle (x, y and radius), BBOX (x, y and radius, OR BBOX coordinates), or own geometry (polygon).

AHN data is obtained from the AHN3 (default), AHN2 or AHN 1 from the available resolutions.

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 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.

+

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 elevations 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.

@@ -137,11 +137,11 @@

Arguments diff --git a/docs/reference/ahn_point.html b/docs/reference/ahn_point.html index 2b8d3ec..a7a752c 100644 --- a/docs/reference/ahn_point.html +++ b/docs/reference/ahn_point.html @@ -2,10 +2,10 @@ AHN elevation point — ahn_point • rAHNextractAHN elevation point — ahn_point • rAHNextract @@ -55,7 +55,7 @@

Get elevation of specific point location.

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).

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 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.

+

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 elevations 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.

@@ -124,11 +124,11 @@

Arguments diff --git a/man/ahn_area.Rd b/man/ahn_area.Rd index 6606418..f54ffd0 100644 --- a/man/ahn_area.Rd +++ b/man/ahn_area.Rd @@ -49,9 +49,9 @@ ahn_area( \item{sheets.method}{Default FALSE. When TRUE, it downloads AHN area through the available GeoTIFF AHN sheets available on \href{https://www.pdok.nl/atom-downloadservices/-/article/actueel-hoogtebestand-nederland-ahn}{PDOK}.} -\item{sheets.dir}{Default is the 'AHN_sheets' directory in the working directory or \code{output.dir} (when defined). In this sheets directory all the AHN sheets are downloaded. If sheets already exists in this directory, they will be reused and no redownload will take place. Within this sheets directory, a directory structure will be created automatically: if \code{sheets.dir} is set to 'myFolder', all sheets will be downloaded in their apropriate AHN version and DEM folder. Eg. '\code{output.dir}/myFolder/AHN_sheets/AHN4/DSM'. You want to use this parameter if you want to have a central location where all the AHN sheets are stored so that they can be accessed elsewhere for other purposes. It is recommended to always use the their original file name after download. Due to the size of these GeoTiff sheets, it is not allowed to download this is the RAM memory (tempdir).} +\item{sheets.dir}{Default is the 'AHN_sheets' directory in the working directory or \code{output.dir} (when defined). In this sheets directory all the AHN sheets are downloaded. If sheets already exists in this directory, they will be reused and no redownload will take place. Within this sheets directory, a directory structure will be created automatically: if \code{sheets.dir} is set to 'myFolder', all sheets will be downloaded in their appropriate AHN version and DEM folder. Eg. '\code{output.dir}/myFolder/AHN_sheets/AHN4/DSM'. You want to use this parameter if you want to have a central location where all the AHN sheets are stored so that they can be accessed elsewhere for other purposes. It is recommended to always use the their original file name after download. Due to the size of these GeoTiff sheets, it is not allowed to download this is the RAM memory (tempdir).} -\item{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.} +\item{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 extraction will be followed.} } \value{ GeoTIFF file of AHN area @@ -59,13 +59,13 @@ GeoTIFF file of AHN area \description{ Get AHN raster area of a specified area. -AHN Area is retrieved through a circle (x, y and radius), BBOX (x, y and radius, OR bbox coordinates), or own geometry (polygon). +AHN Area is retrieved through a circle (x, y and radius), BBOX (x, y and radius, OR BBOX coordinates), or own geometry (polygon). AHN data is obtained from the AHN3 (default), AHN2 or AHN 1 from the available resolutions. Available for the Digital Surface Model (DSM) and Digital Terrain Model (DTM). -You can download the AHN data using the WCS method (default, \code{sheets.method = FALSE}) returning a GeoTIFF format. The sheets method (\code{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. +You can download the AHN data using the WCS method (default, \code{sheets.method = FALSE}) returning a GeoTIFF format. The sheets method (\code{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 elevations 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. } \author{ Jelle Stuurman diff --git a/man/ahn_point.Rd b/man/ahn_point.Rd index 6671b2d..fd3b011 100644 --- a/man/ahn_point.Rd +++ b/man/ahn_point.Rd @@ -43,9 +43,9 @@ ahn_point( \item{sheets.method}{Default FALSE. When TRUE, it downloads AHN area through the available GeoTIFF AHN sheets available on \href{https://www.pdok.nl/atom-downloadservices/-/article/actueel-hoogtebestand-nederland-ahn}{PDOK}.} -\item{sheets.dir}{Default is the 'AHN_sheets' directory in the working directory or \code{output.dir} (when defined). In this sheets directory all the AHN sheets are downloaded. If sheets already exists in this directory, they will be reused and no redownload will take place. Within this sheets directory, a directory structure will be created automatically: if \code{sheets.dir} is set to 'myFolder', all sheets will be downloaded in their apropriate AHN version and DEM folder. Eg. '\code{output.dir}/myFolder/AHN_sheets/AHN4/DSM'. You want to use this parameter if you want to have a central location where all the AHN sheets are stored so that they can be accessed elsewhere for other purposes. It is recommended to always use the their original file name after download. Due to the size of these GeoTiff sheets, it is not allowed to download this is the RAM memory (tempdir).} +\item{sheets.dir}{Default is the 'AHN_sheets' directory in the working directory or \code{output.dir} (when defined). In this sheets directory all the AHN sheets are downloaded. If sheets already exists in this directory, they will be reused and no redownload will take place. Within this sheets directory, a directory structure will be created automatically: if \code{sheets.dir} is set to 'myFolder', all sheets will be downloaded in their appropriate AHN version and DEM folder. Eg. '\code{output.dir}/myFolder/AHN_sheets/AHN4/DSM'. You want to use this parameter if you want to have a central location where all the AHN sheets are stored so that they can be accessed elsewhere for other purposes. It is recommended to always use the their original file name after download. Due to the size of these GeoTiff sheets, it is not allowed to download this is the RAM memory (tempdir).} -\item{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.} +\item{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 extraction will be followed.} } \value{ AHN elevation in meters. @@ -57,5 +57,5 @@ Requires the X and Y coordinates as input. AHN data is obtained from the AHN in 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 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. +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 elevations 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. }