R/sf_utils.R
In ggRaver/Lslide: Functions for object-oriented image analysis in landslide science
Defines functions
st_erase
Documented in
st_erase
#' Erase one geometry from another
#'
#' This function erase one geometry from another. The projection must be identical.
#'
#' @param x object of class sf
#' @param y object of class sf
#' @importFrom magrittr "%>%"
#' @return
#' Geometry of class sfc
#'
#'
#' @keywords simple feature, erase
#'
#'
#' @export
#'
st_erase = function(x, y, precision = 0, do.subset = TRUE)
{
if(do.subset)
{
inter <- sf::st_intersects(x = y, y = x) %>% unlist(.) %>% unique(.)
x.remain <- x[-inter,]
x <- x[inter,]
}
if(precision != 0)
{
x <- x %>% sf::st_set_precision(x = ., precision = precision) %>% lwgeom::st_make_valid(.)
y <- y %>% sf::st_combine(.) %>% sf::st_union(.) %>%
sf::st_set_precision(x = ., precision = precision) %>% lwgeom::st_make_valid(.)
}
out <- sf::st_difference(x = x, y = y) # erase y from x
if(do.subset)
{
out <- rbind(out, x.remain)
}
return(out)
}
#' Erase one geometry from another using Saga GIS
#'
#' This function erase one geometry from another. The projection must be identical.
#'
#' @param x object of class sf. First element: Should be either line or polygon
#' @param y object of class sf. Second element: Always polygon.
#' @param method method of erase. Either "1": , or "2": Line-Polygon Intersection. Default: "1"
#' @param split Set to "1", if multi-part polygons should be splitted to single-part polygons. Default: "0"
#' @param attributes attributes inherited to intersection result. [0] polygon, [1] line, [2] line and polygon. Default: "1"
#' @param env.rsaga SAGA GIS environemnt. Default: RSAGA::rsaga.env()
#' @param check.geom If set to TRUE then geometry is checked with sf::st_is_valid(). If there are invalid geometries, geometries are repaired using lwgeom::st_make_valid(). Default: TRUE
#' @importFrom magrittr "%>%"
#' @return
#' Geometry of class sfc
#'
#'
#' @keywords simple feature, erase
#'
#'
#' @export
#'
rsaga_erase = function(x, y, method = "1", split = "0", attributes = "1", env.rsaga = RSAGA::rsaga.env(), check.geom = TRUE, quiet = TRUE)
{
path.x <- file.path(tempdir(), "tmp_x.shp")
path.y <- file.path(tempdir(), "tmp_y.shp")
path.result <- file.path(tempdir(), "tmp_result.shp")
sf::st_write(obj = x, dsn = path.x, delete_layer = TRUE, quiet = quiet)
sf::st_write(obj = y, dsn = path.y, delete_layer = TRUE, quiet = quiet)
if(method == "1")
{
# RSAGA::rsaga.get.usage(lib = "shapes_polygons", module = 15, env = env.rsaga)
RSAGA::rsaga.geoprocessor(lib = "shapes_polygons", module = 15, env = env.rsaga, show.output.on.console = !quiet, param = list(
A = path.x, B = path.y, RESULT = path.result, SPLIT = split))
}
if(method == "2")
{
# RSAGA::rsaga.get.usage(lib = "shapes_lines", module = 3, env = env.rsaga)
# ATTRIBUTES: [1] line
RSAGA::rsaga.geoprocessor(lib = "shapes_lines", module = 3, env = env.rsaga, show.output.on.console = !quiet, param = list(
LINES = path.x, POLYGONS = path.y, ATTRIBUTES = attributes, DIFFERENCE = path.result))
}
## read data
out <- sf::st_read(dsn = path.result, quiet = quiet)
## check validity
if(check.geom && !all(sf::st_is_valid(out)))
{
warning('Some invalid geometries by "rsaga_erase". Try to correct geomeries using lwgeom::st_make_valid()!')
out <- lwgeom::st_make_valid(x = out)
if(method == "1")
{
out <- suppressWarnings(out %>% sf::st_collection_extract(x = ., type = c("POLYGON")))
}
if(method == "2")
{
out <- suppressWarnings(out %>% sf::st_collection_extract(x = ., type = c("LINESTRING")))
}
}
return(out)
}
#' Return geometry from bounding box
#'
#' This function calculate the geometry based on a bounding box.
#'
#' @param x object of class sf
#' @param extent vector containing numeric values, in the following order: xmin, xmax, ymax, ymin
#' @importFrom magrittr "%>%"
#' @return
#' Geometry of class sfc
#'
#'
#' @keywords simple feature, geometry of bounding box
#'
#'
#' @export
#'
st_bbox_geom = function(x, extent = NULL)
{
if(!is.null(x) & !is.null(extent))
{
stop("Input conflict: either x or extent!")
}
if(is.null(extent))
{
out <- sf::st_bbox(x) %>% sf::st_as_sfc(.)
} else {
names(extent) <- c("xmin", "xmax", "ymax", "ymin")
out <- sf::st_bbox(extent) %>% sf::st_as_sfc(.)
}
return(out)
} # end of function st_bbox_geom
#' Plot sf geometry
#'
#' This function simply plots the geometry of a sf object.
#'
#' @param x object of class sf
#' @return
#' plot
#'
#'
#' @keywords simple feature, plot
#'
#'
#' @export
#'
st_plot = function(x) plot(sf::st_geometry(x))
#' Dissolve geometry
#'
#' This function dissovle a geometry based on a field. Optionnally, field statistics can be computed.
#'
#' @param x object of class sf
#' @param by field for dissolving. Default: NULL
#' @param ... Optional: field statistcs
#' @importFrom magrittr "%>%"
#' @return
#' dissolved geometry of class sf
#'
#'
#' @keywords simple feature, dissolve
#'
#'
#' @export
#'
st_dissolve = function(x, by = NULL, ...) x %>% dplyr::group_by(.dots = by) %>% dplyr::summarise(...)
#' Row bind simple features
#'
#' This function combines simple feature with different columns using dplr::bind_rows().
#'
#' @param x list of objects of class sf
#' @param geom vector containing the geometries of x. Order must be the same as in list.
#' @importFrom magrittr "%>%"
#' @return
#' sf object with combined data.frames
#'
#'
#' @keywords simple feature, rbind
#'
#'
#' @export
#'
st_rbind = function(x, geom) x %>% lapply(X = ., FUN = function(x) sf::st_set_geometry(x = x, value = NULL)) %>%
dplyr::bind_rows() %>% sf::st_set_geometry(x = ., value = geom)
#' Calculate the perimeter of a polygon
#'
#' This function calculates the perimeter of a polygon.
#'
#' @param x object of class sf
#' @importFrom magrittr "%>%"
#' @return
#' Vector with perimeter values
#'
#'
#' @keywords simple feature, perimeter
#'
#'
#' @export
#'
st_perimeter = function(x, ...) suppressWarnings(x %>% sf::st_cast(x = ., to = "MULTILINESTRING", ...) %>% sf::st_length(.) %>% as.numeric(.)) # units::drop_units(.)
#' Calculate the shape indices of a polygon
#'
#' This function calculates the shape index of a polygon.
#'
#' @param x object of class sf
#' @importFrom magrittr "%>%"
#' @return
#' Vector with shape index and interior edge ratio values
#'
#'
#' @keywords simple feature, shape index, interior edge ratio
#'
#'
#' @export
#'
st_shape_indices = function(x){
perimeter <- x %>% st_perimeter(.)
area <- x %>% sf::st_area(.) %>% as.numeric(.) # units::drop_units(.)
shape_index <- perimeter / (2 * sqrt(pi * area))
interior_edge_ratio <- perimeter / area
return(list(shape_index = shape_index, interior_edge_ratio = interior_edge_ratio))
}
#' Calculate integration index
#'
#' This function calculates an index defining the intergration of new land use area into conisting land use area
#' of the same land use class conisdering two different time steps.
#'
#' @param geom.old object of class sf representing a land use class of a previous time step
#' @param geom.new object of class sf representing a land use class of a following time step
#' @param geom.boundary polygon of class sf representing subregions, e.g. administrative boundaries
#' @param tol tolerance value for overlapping area m square
#' @param precision precision for process. See sf::st_set_precision(). Default: 0
#' @param env.rsaga environment of SAGA GIS. If st_erase fails then SAGA GIS erase is used. Default: NULL, but in function call if not set: RSAGA::rsaga.env()
#' @param use.saga use SAGA GIS for erase process. Default: FALSE
#' @param return.geom If set to TRUE, intermediate geometries are returned as well. Default: FALSE
#' @param quiet show output on console. Default: FALSE
#' @importFrom magrittr "%>%"
#' @return
#' Vector with integration index (completly integrated: 2/3 < R < 1; good integrated: 1/3 < R < 2/3; low integrated: 0 < R < 1/3; not integrated: 0)
#'
#'
#' @keywords simple feature, shape index, interior edge ratio
#'
#'
#' @export
#'
st_integration_index = function(geom.old, geom.new, geom.boundary = NULL, tol = 0.1, precision = 0,
ignr.overlap = FALSE, env.rsaga = NULL, use.saga = FALSE, return.geom = FALSE, quiet = FALSE){
# get start time of process
process.time.start <- proc.time()
## check input
if(missing(geom.old) || missing(geom.new)){ stop('Input is missing!')}
## check and set precision
if(!quiet) cat("... check and set st_precision \n")
if(sf::st_precision(geom.old) != precision)
{
geom.old <- geom.old %>% sf::st_set_precision(x = ., precision = precision) %>%
lwgeom::st_make_valid(.) %>%
sf::st_collection_extract(x = ., type = "POLYGON")
}
if(sf::st_precision(geom.new) != precision)
{
geom.new <- geom.new %>% sf::st_set_precision(x = ., precision = precision) %>%
lwgeom::st_make_valid(.) %>%
sf::st_collection_extract(x = ., type = "POLYGON")
}
## check validity of geometries
if(!all(sf::st_is_valid(geom.old))){ stop('Input of "geom.old" contains not valid geometries. Please try lwgeom::st_make_valid().')}
if(!all(sf::st_is_valid(geom.new))){ stop('Input of "geom.new" contains not valid geometries. Please try lwgeom::st_make_valid().')}
if(!is.null(geom.boundary) && !all(sf::st_is_valid(geom.boundary))){ stop('Input of "geom.boundary" contains not valid geometries. Please try lwgeom::st_make_valid().')}
# if(!quiet) cat("... union input geometries \n")
# geom.old <- sf::st_union(x = geom.old) %>% sf::st_cast(., "POLYGON") %>% sf::st_set_precision(x = ., precision = precision)
# geom.new <- sf::st_union(x = geom.new) %>% sf::st_cast(., "POLYGON") %>% sf::st_set_precision(x = ., precision = precision)
# # # # START CALCULATION OF INTEGRATION INDEX
## check for overlapping polygon
## common area
if(!quiet) cat("... intersection of input geometries \n")
inter <- suppressWarnings(sf::st_intersection(x = geom.old, y = geom.new) %>%
sf::st_collection_extract(x = ., type = c("POLYGON")))# %>%
# sf::st_set_precision(x = ., precision = precision)
## new area
if(!quiet) cat('... erase intersection from "geom.new" (this can take a while!) \n')
if(use.saga)
{
if(is.null(env.rsaga))
{
env.rsaga <- RSAGA::rsaga.env()
}
erase <- Lslide::rsaga_erase(x = geom.new, y = inter, method = "1",
split = "1", env.rsaga = env.rsaga) %>%
.[which(x = as.numeric(sf::st_area(.)) >= tol),]
} else {
erase <- tryCatch({
suppressWarnings(Lslide::st_erase(x = geom.new, y = inter, precision = precision) %>%
sf::st_collection_extract(x = ., type = c("POLYGON")) %>%
sf::st_cast(x = ., to = "POLYGON") %>%
.[which(x = as.numeric(sf::st_area(.)) >= tol),])
}, error = function(e){
warning(paste('SAGA GIS is used due to error in st_erase():', e))
if(is.null(env.rsaga))
{
env.rsaga <- RSAGA::rsaga.env()
}
tmp.erase <- Lslide::rsaga_erase(x = geom.new, y = inter, method = "1",
split = "1", env.rsaga = env.rsaga) %>%
.[which(x = as.numeric(sf::st_area(.)) >= tol),]
return(tmp.erase)
})
} # end of use.saga
if(!quiet) cat('... conversion to lines \n')
line.erase <- sf::st_cast(x = erase, to = "MULTILINESTRING") %>% sf::st_set_precision(x = ., precision = precision)
line.inter <- sf::st_cast(x = inter, to = "MULTILINESTRING") %>% sf::st_set_precision(x = ., precision = precision)
## common border
if(!quiet) cat('... find border lines by intersection \n')
unique.border <- suppressWarnings(sf::st_intersection(x = line.inter, y = line.erase) %>%
sf::st_collection_extract(x = ., type = c("LINESTRING")))
## check boundary constraints
if(!is.null(geom.boundary))
{
geom.boundary$ID_BOUNDS <- 1:nrow(geom.boundary)
geom.boundary <- geom.boundary[, c("ID_BOUNDS", "geometry")]
if(!quiet) cat('... intersection with boundaries \n')
unique.border <- suppressWarnings(sf::st_intersection(x = geom.boundary, y = unique.border) %>%
sf::st_collection_extract(x = ., type = c("LINESTRING")) %>%
st_cast(., "MULTILINESTRING") %>% sf::st_cast(., "LINESTRING")) # cast is necessairy to split multi-object
dt.unique.border <- unique.border %>%
sf::st_set_geometry(x = ., value = NULL) %>%
data.table::as.data.table(.)
erase <- suppressWarnings(sf::st_intersection(x = geom.boundary, y = erase) %>%
sf::st_collection_extract(x = ., type = c("POLYGON")) %>%
st_cast(., "MULTIPOLYGON") %>% sf::st_cast(., "POLYGON"))
dt.erase <- erase %>%
sf::st_set_geometry(x = ., value = NULL) %>%
data.table::as.data.table(.)
if(!quiet) cat('... get statistics and calculate index \n')
dt.unique.border$L <- sf::st_length(x = unique.border) %>% as.numeric() # units::drop_units(x = .)
dt.erase$P <- Lslide::st_perimeter(x = erase)
dt.result.UB <- dt.unique.border[,list(L = sum(L, na.rm = TRUE)), by = ID_BOUNDS]
dt.result.E <- dt.erase[,list(P = sum(P, na.rm = TRUE)), by = ID_BOUNDS]
dt.result <- merge(x = dt.result.E, y = dt.result.UB, by = "ID_BOUNDS") %>%
dplyr::mutate(.data = ., R = L/P)
} else {
if(!quiet) cat('... get statistics and calculate index \n')
dt.result <- data.table::data.table(P = sum(Lslide::st_perimeter(x = erase), na.rm = TRUE),
L = sum((sf::st_length(x = unique.border) %>% as.numeric(.)), na.rm = TRUE)) %>% # units::drop_units(x = .))) %>%
dplyr::mutate(.data = ., R = L/P)
} # end of if-else boundary check
process.time.run <- proc.time() - process.time.start
if(quiet == FALSE) cat("------ Run of st_integration_index: " , round(process.time.run["elapsed"][[1]]/60, digits = 4), " Minutes \n")
if(return.geom)
{
return(list(integration_index = dt.result, unique_border = unique.border, new_area = erase))
} else {
return(dt.result)
}
} # end of function st_integration_index
#' Calculate the effective mesh size
#'
#' This function calculates the effective mesh size.
#'
#' @param geom.frag polygon of class sf representing the fragmentation geometry
#' @param geom.boundary polygon of class sf representing subregions, e.g. administrative boundaries
#' @param total.area Numeric value representing size for area. Only to use if geom.boundary is not present. Value must match with the conversion constant c (default hectare). Default: NULL
#' @param conv constant to convert original square m output. Default: 10000 to convert to hectare. If set to 1, than meter square is the result.
#' @param do.preProcessing If TRUE (default), the input of geom.frag is, first, dissolved to single part feature, and second, splitted to multi-parts. By this step it is assured, that polygon connected to each other are summarized
#' @param return.geom If set to TRUE, intermediate geometries are returned as well. Default: FALSE
#' @param quiet If set to FALSE, actual state is printed to console. Default: TRUE.
#' @importFrom magrittr "%>%"
#' @return
#' Vector with shape index and interior edge ratio values
#'
#'
#' @keywords simple feature, mesh, effective mesh size
#'
#'
#' @export
#'
st_mesh = function(geom.frag, geom.boundary = NULL, total.area = NULL, conv = 10000, do.preProcessing = TRUE, return.geom = FALSE, quiet = TRUE)
{
# get start time of process
process.time.start <- proc.time()
## check input
if(missing(geom.frag)){ stop('Input of "geom.frag" is missing!')}
if(!missing(geom.frag) && !("sf" %in% class(geom.frag))){ stop('Input of "geom.frag" is not of class "sf"!')}
if(!is.null(geom.boundary) && !("sf" %in% class(geom.boundary))){ stop('Input of "geom.boundary" is not of class "sf"!')}
if(is.null(geom.boundary) && (is.null(total.area) || !is.numeric(total.area))){
stop('If input of "geom.boundary" is null, then input of "total.area" must be a numeric value representing the size of the study area.')
}
## check validity of geometries
if(!all(sf::st_is_valid(geom.frag))){ stop('Input of "geom.frag" contains not valid geometries. Please try lwgeom::st_make_valid().')}
if(!is.null(geom.boundary) && !all(sf::st_is_valid(geom.boundary))){ stop('Input of "geom.boundary" contains not valid geometries. Please try lwgeom::st_make_valid().')}
if(do.preProcessing)
{
if(!quiet) cat("... union geometries to a single geometry with resolved boundaries \n")
geom.frag <- geom.frag %>% sf::st_union(.)
if(!quiet) cat("... split multi-parts to single-parts polygon \n")
geom.frag <- geom.frag %>% st_cast(., "POLYGON") %>% sf::st_sf(ID_FRAG = 1:length(.), geometry = .)
} else {
geom.frag$ID_FRAG <- 1:nrow(geom.frag) ## add unique IDs
}
if(!is.null(geom.boundary)){ geom.boundary$ID_BOUNDS <- 1:nrow(geom.boundary) } # ## add unique IDs
## add Area in m_sq
geom.frag$A_FRAG <- sf::st_area(geom.frag) %>% (function(x, conv = conv) as.numeric(x)/conv) # as.numeric() %>% "/" (conv) # units::drop_units(.)
if(!is.null(geom.boundary)){ geom.boundary$A_BOUNDS <- sf::st_area(geom.boundary) %>% (function(x, conv = conv) as.numeric(x)/conv)}
## subset data
geom.frag <- geom.frag[, c("ID_FRAG", "A_FRAG", "geometry")]
if(!is.null(geom.boundary)){geom.boundary <- geom.boundary[, c("ID_BOUNDS", "A_BOUNDS", "geometry")]}
# # # CALCULATE MESH SIZE INDICES
## start calculation
if(!is.null(geom.boundary))
{
## get intersection
if(!quiet) cat("... intersection to boundary \n")
inter <- suppressWarnings(sf::st_intersection(x = geom.boundary, y = geom.frag))
inter <- suppressWarnings(inter %>% st_cast(., "MULTIPOLYGON") %>% sf::st_cast(., "POLYGON")) # cast is necessairy to split multi-polygons
inter$A_FRAG_INTER <- sf::st_area(inter) %>% (function(x, conv = conv) as.numeric(x)/conv) # overwrite area of fragments, # units::drop_units(.)
## calculation of mesh indices
df.inter <- sf::st_set_geometry(x = inter, value = NULL) %>% data.table::as.data.table(.)
df.inter.multi <- df.inter[, list(A_BOUNDS = unique(A_BOUNDS), # splitted multi-parts are joined together
A_FRAG = unique(A_FRAG),
A_FRAG_INTER = sum(A_FRAG_INTER, na.rm = TRUE)), by = list(ID_BOUNDS, ID_FRAG)]
if(!quiet) cat("... get statistics \n")
## cutting-out (CUT) procedure
mesh.CUT <- df.inter[, list(Fg = unique(A_BOUNDS),
Fi = sum(A_FRAG_INTER^2, na.rm = TRUE),
Fi_count = length(!is.na(A_FRAG_INTER)),
Fi_sum = sum(A_FRAG_INTER, na.rm = TRUE),
Fi_min = min(A_FRAG_INTER, na.rm = TRUE),
Fi_max = max(A_FRAG_INTER, na.rm = TRUE),
Fi_ave = mean(A_FRAG_INTER, na.rm = TRUE),
Ci = sum((A_FRAG_INTER/A_BOUNDS)^2, na.rm = TRUE)), by = ID_BOUNDS] %>%
dplyr::mutate(., Di = 1-Ci) %>%
dplyr::mutate(., Si = 1/Ci) %>%
dplyr::mutate(., mEff_CUT = Fi/Fg)
# View(mesh.CUT)
## cross-boundary connections (CBC) procedure
mesh.CBC <- df.inter.multi[, list(Fg = unique(A_BOUNDS),
FiErg_count = length(!is.na(A_FRAG)),
FiErg_sum = sum(A_FRAG, na.rm = TRUE),
FiErg_min = min(A_FRAG, na.rm = TRUE),
FiErg_max = max(A_FRAG, na.rm = TRUE),
FiErg_ave = mean(A_FRAG, na.rm = TRUE),
CiErg = sum((A_FRAG/A_BOUNDS)^2, na.rm = TRUE),
Fi_CBC1 = sum(A_FRAG_INTER*A_FRAG, na.rm = TRUE),
Fi_CBC2 = sum(2*A_FRAG_INTER*A_FRAG-A_FRAG_INTER^2, na.rm = TRUE)), by = ID_BOUNDS] %>%
dplyr::mutate(., DiErg = 1-CiErg) %>%
dplyr::mutate(., SiErg = 1/CiErg) %>%
dplyr::mutate(., mEff_CBC1 = Fi_CBC1/Fg) %>%
dplyr::mutate(., mEff_CBC2 = Fi_CBC2/Fg)
# View(mesh.CBC)
df.result <- merge(x = mesh.CUT, y = mesh.CBC[, -2], by = "ID_BOUNDS")
} else {
if(!quiet) cat("... get statistics \n")
df.result <- sf::st_set_geometry(x = geom.frag, value = NULL) %>% data.table::as.data.table(.)
df.result <- df.result[, list(Fg = total.area,
Fi = sum(A_FRAG^2, na.rm = TRUE),
Fi_count = length(!is.na(A_FRAG)),
Fi_sum = sum(A_FRAG, na.rm = TRUE),
Fi_min = min(A_FRAG, na.rm = TRUE),
Fi_max = max(A_FRAG, na.rm = TRUE),
Fi_ave = mean(A_FRAG, na.rm = TRUE),
Ci = sum((A_FRAG/total.area)^2, na.rm = TRUE)),] %>%
dplyr::mutate(., Di = 1-Ci) %>%
dplyr::mutate(., Si = 1/Ci) %>%
dplyr::mutate(., mEff_CUT = Fi/total.area)
}
process.time.run <- proc.time() - process.time.start
if(quiet == FALSE) cat("------ Run of st_mesh: " , round(process.time.run["elapsed"][[1]]/60, digits = 4), " Minutes \n")
if(return.geom)
{
if(is.null(geom.boundary))
{
return(list(mesh = df.result, geom.frag = geom.frag))
} else {
return(list(mesh = df.result, geom.frag = geom.frag, geom.inter = inter))
}
} else {
return(df.result)
}
} # end of function st_mesh
#' Calculate urban sprawl
#'
#' This function calculates the urban sprawl of urban area.
#'
#' @param geom.urban polygon of class sf representing the fragmentation geometry
#' @param geom.boundary polygon of class sf representing subregions, e.g. administrative boundaries
#' @param dist vector containing distance between lines in x and y direction. Default: c(100, 100) [m]
#' @param trans transformation function {x-1+1/(x+trans.k)} with x as free line and trans.k a constant
#' @param trans.k constant in km for transformation function trans. Default: 1
#' @param tol tolerance value for intersection with erased lines. Buffering procedure is used.Default: 0.1 [m]
#' @param precision precision for process. See sf::st_set_precision(). Default: 0
#' @param extent Numeric value representing extent for area. Format of vector: c(xmin, xmax, ymax, ymin) . Default: NULL
#' @param force.extent If TRUE extent is used instead of geom.boundary (if both are present). Default: FALSE
#' @param do.preProcessing If TRUE (default), the input of geom.frag is, first, dissolved to single part feature, and second, splitted to multi-parts. By this step it is assured, that polygon connected to each other are summarized
#' @param return.geom If set to TRUE, intermediate geometries are returned as well. Default: FALSE
#' @param env.rsaga environment of SAGA GIS. If st_erase fails then SAGA GIS erase is used. Default: NULL, but in function call if not set: RSAGA::rsaga.env()
#' @param use.saga use SAGA GIS for erase process. Default: FALSE
#' @param quiet If set to FALSE, actual state is printed to console. Default: TRUE.
#' @importFrom magrittr "%>%"
#' @return
#' strong urban sprawl: 40-50%, less urban sprawl: 80-90%
#'
#'
#' @keywords simple feature, urban sprawl
#'
#'
#' @export
#'
st_urban_sprawl = function(geom.urban, geom.boundary = NULL, dist = c(100, 100), trans = function(x, trans.k){x-1+1/(x+trans.k)}, trans.k = 1, tol = 0.1, extent = NULL, force.extent = FALSE,
precision = 0, do.preProcessing = TRUE, env.rsaga = NULL, use.saga = FALSE, return.geom = FALSE, quiet = TRUE)
{
# get start time of process
process.time.start <- proc.time()
## check input
if(missing(geom.urban)){ stop('Input of "geom.urban" is missing!')}
if(!missing(geom.urban) && !("sf" %in% class(geom.urban))){ stop('Input of "geom.urban" is not of class "sf"!')}
if(!is.null(geom.boundary) && !("sf" %in% class(geom.boundary))){ stop('Input of "geom.boundary" is not of class "sf"!')}
if(is.null(geom.boundary) && is.null(extent)){
warning('If input of "geom.boundary" and "extent" is null. Fish net is created using the extent of "geom.urban"!')
}
if(is.null(geom.boundary) & is.null(extent) & force.extent){
stop('If "force.extent" is TRUE, than extent should be set!"!')
}
if(!is.null(geom.boundary) & !is.null(extent)){
warning('Fish net is created using extent of "geom.boundary". "Extent" is skipped. For "extent" use "force.extent"!')
}
if(sf::st_precision(geom.urban) != precision)
{
geom.urban <- geom.urban %>% sf::st_set_precision(x = ., precision = precision) %>%
lwgeom::st_make_valid(.) %>%
sf::st_collection_extract(x = ., type = "POLYGON")
}
## check validity of geometries
if(!all(sf::st_is_valid(geom.urban))){ stop('Input of "geom.urban" contains not valid geometries. Please try lwgeom::st_make_valid().')}
if(!is.null(geom.boundary) && !all(sf::st_is_valid(geom.boundary))){ stop('Input of "geom.boundary" contains not valid geometries. Please try lwgeom::st_make_valid().')}
## add unique ID and subset data
if(!is.null(geom.boundary)){ geom.boundary$ID_BOUNDS <- 1:nrow(geom.boundary) } # ## add unique IDs
if(!is.null(geom.boundary)){geom.boundary <- geom.boundary[, c("ID_BOUNDS", "geometry")]}
## ... create boundary for fishnet
if(!is.null(geom.boundary) & !force.extent){bbox.fishnet <- Lslide::st_bbox_geom(x = geom.boundary)
} else if(!is.null(geom.boundary) & force.extent){bbox.fishnet <- Lslide::st_bbox_geom(extent = extent)
} else if(is.null(geom.boundary) & !is.null(extent)){bbox.fishnet <- Lslide::st_bbox_geom(extent = extent)
} else {bbox.fishnet <- Lslide::st_bbox_geom(x = geom.urban) }
## create and subset fish net
if(!quiet) cat("... create fishnet \n")
fishnet <- Lslide::st_make_grid_lines(x = bbox.fishnet, cellsize = dist)
if(!is.null(geom.boundary) & !force.extent){
if(!quiet) cat("... subset fishnet to area of interest \n")
fishnet <- suppressWarnings(sf::st_intersection(x = fishnet, y = geom.boundary)) %>%
sf::st_collection_extract(x = ., type = "LINESTRING", warn = FALSE)
} else {
fishnet <- suppressWarnings(sf::st_intersection(x = fishnet, y = bbox.fishnet)) %>%
sf::st_collection_extract(x = ., type = "LINESTRING", warn = FALSE)
}
# fishnet <- fishnet[, c("ID", "geometry")]
## do post-processing of urban geometry
if(do.preProcessing)
{
if(!quiet) cat("... union geometries to a single geometry with resolved boundaries \n")
geom.urban <- geom.urban %>% sf::st_union(.)
if(!quiet) cat("... split multi-parts to single-parts polygon \n")
geom.urban <- geom.urban %>% st_cast(., "POLYGON") %>% sf::st_sf(ID_URBAN = 1:length(.), geometry = .)
} else {
geom.urban$ID_URBAN <- 1:nrow(geom.urban) ## add unique IDs
}
## check intersection
inter.check <- sf::st_intersects(x = geom.urban, y = fishnet)
inter.empty <- which(sapply(inter.check, function(x) length(x) == 0))
if(length(inter.empty) > 0)
{
inter.A <- geom.urban[inter.empty,] %>% sf::st_area(.) %>% as.numeric(.)
inter.num <- length(inter.empty)
warning('Some urban polygons are not intersected by fishnet: ', inter.num,
' | area min: ', min(inter.A, na.rm = TRUE), ' - max: ', max(inter.A, na.rm = TRUE), ' [m_sqr] \n')
}
## erase urban area from fishnet using SAGA GIS
if(!quiet) cat("... erase urban area from fishnet \n")
if(use.saga)
{
if(is.null(env.rsaga))
{
env.rsaga <- RSAGA::rsaga.env()
}
erase <- Lslide::rsaga_erase(x = fishnet, y = geom.urban, method = "2", env.rsaga = env.rsaga)
} else {
erase <- suppressWarnings(Lslide::st_erase(x = fishnet, y = geom.urban, precision = precision) %>%
sf::st_collection_extract(x = ., type = "LINESTRING"))
}
## split to single part
erase.single <- erase %>% sf::st_cast(x = ., to = 'LINESTRING', warn = FALSE)
## selection of lines
if(!quiet) cat("... selection of lines \n")
lines.urban <- sf::st_intersects(x = geom.urban %>%
sf::st_buffer(x = ., dist = tol),
y = erase.single) %>%
unlist(.) %>% unique(.)
if(!is.null(geom.boundary) & !force.extent)
{
lines.boundary <- sf::st_intersects(x = geom.boundary %>%
sf::st_boundary(x = .) %>%
sf::st_buffer(x = ., dist = tol),
y = erase.single) %>%
unlist(.) %>% unique(.)
} else {
lines.boundary <- sf::st_intersects(x = bbox.fishnet %>%
sf::st_boundary(x = .) %>%
sf::st_buffer(x = ., dist = tol),
y = erase.single) %>%
unlist(.) %>% unique(.)
}
# lines.urban %>% length(.)
# lines.boundary %>% length(.)
check.missing <- setdiff(c(1:nrow(erase.single)),
(c(lines.urban, lines.boundary) %>% unique(.)))
if(length(check.missing) > 0)
{
warning("Some lines are neither intersected by boundary nor by urban area: ", length(check.missing), " in total \n")
}
erase.single$Type <- NA
# TPYE 1: Lines between urban area (red lines)
type.urban <- setdiff(x = lines.urban, y = lines.boundary)
erase.single$Type[type.urban] <- 1
# TPYE 2: Lines not touching urban area (blue lines)
type.boundary <- setdiff(x = lines.boundary, y = lines.urban)
erase.single$Type[type.boundary] <- 2
# TYPE 3: Lines between boundary and urban area (green lines)
type.bound.urb <- intersect(x = lines.urban, y = lines.boundary)
erase.single$Type[type.bound.urb] <- 3
# summary(erase.single$Type)
# c(type.bound.urb, type.boundary, type.urban) %>% unique(.) %>% length(.)
# browser()
## group single lines back to multi-lines
if((!is.null(geom.boundary) & !force.extent))
{
erase.final <- st_dissolve(x = erase.single[which(erase.single$Type == 2 | erase.single$Type == 3),],
by = list("Type", "ID_FNET", "ID_BOUNDS"))
erase.final <- rbind(erase.final, erase.single[which(erase.single$Type == 1),])
} else {
erase.final <- st_dissolve(x = erase.single[which(erase.single$Type == 2 | erase.single$Type == 3),],
by = list("Type", "ID_FNET"))
erase.final <- rbind(erase.final, erase.single[which(erase.single$Type == 1),])
}
erase.final <- erase.final %>% sf::st_collection_extract(x = ., type = "LINESTRING")
erase.final$L <- sf::st_length(x = erase.final) %>% (function(x) as.numeric(x)/1000) # as.numeric(.) %>% "/" (1000) # convert from meter to km!
erase.final$L_trans <- ifelse(erase.final$Type == 2, erase.final$L, trans(x = erase.final$L, trans.k = trans.k))
## summarizing final data
if(!quiet) cat("... get statistics: urban sprawl \n")
if((!is.null(geom.boundary) & !force.extent))
{
df.result <- erase.final %>% sf::st_set_geometry(x = ., value = NULL) %>%
data.table::as.data.table(.) %>%
.[,list(L = sum(L, na.rm = TRUE),
L_trans = sum(L_trans, na.rm = TRUE)), by = ID_BOUNDS] %>%
dplyr::mutate(.data = ., FFE = L_trans/L*100)
} else {
df.result <- erase.final %>% sf::st_set_geometry(x = ., value = NULL) %>%
data.table::as.data.table(.) %>%
.[,list(L = sum(L, na.rm = TRUE),
L_trans = sum(L_trans, na.rm = TRUE)),] %>%
dplyr::mutate(.data = ., FFE = L_trans/L*100)
}
process.time.run <- proc.time() - process.time.start
if(quiet == FALSE) cat("------ Run of urban_sprawl: " , round(process.time.run["elapsed"][[1]]/60, digits = 4), " Minutes \n")
if(return.geom)
{
return(list(FFE = df.result, geom = erase.final))
} else {
return(df.result)
}
} # end of function st_urban_sprawl
#' Make multiline fishnet
#'
#' This function creates a fishnet, consisting of multilines in x and y direction.
#'
#' @param x object of class sf
#' @importFrom magrittr "%>%"
#' @return
#' Geometry of class sfc
#'
#'
#' @keywords simple feature, fishnet, multilines
#'
#'
#' @export
#'
st_make_grid_lines = function(x, cellsize)
{
# ... create fishnet using corners (lower left corner is starting point)
fishnet <- sf::st_make_grid(x = x, cellsize = cellsize, what = "corners")
# ... get sf::st_make_grid function parameters (automatically calculated in function call)
offset <- sf::st_bbox(obj = x)[1:2]
nx <- ceiling((sf::st_bbox(obj = x)[3] - offset[1])/cellsize[1]) + 1 # add 1 to fit overall columns
ny <- ceiling((sf::st_bbox(obj = x)[4] - offset[2])/cellsize[2]) + 1 # add 1 to fit overall rows
# ... create lines
linesX <- lapply(1:ny, function(i, nx, fishnet){
if(i == 1)
{
index <- 1:nx
} else {
index <- ((i-1)*nx+1):(i*nx)
}
line <- fishnet[index] %>% st_combine(.) %>%
sf::st_multilinestring(x = ., dim = "XY") %>%
st_sfc(.) %>%
sf::st_sf(ID_FNET = paste0("X", i), geometry = .)
return(line)
}, nx = nx, fishnet = fishnet) %>% do.call(what = rbind, args = .)
seqY <- seq(from = 1, to = length(fishnet), by = nx)
linesY <- lapply(0:(nx-1), FUN = function(i, seqY, fishnet){
if(i == 0){
index <- seqY
} else {
index <- seqY+i
}
line <- fishnet[index] %>% st_combine(.) %>%
sf::st_multilinestring(x = ., dim = "XY") %>%
st_sfc(.) %>%
sf::st_sf(ID_FNET = paste0("Y", i), geometry = .)
return(line)
}, seqY = seqY, fishnet) %>% do.call(what = rbind, args = .)
fishnet.multi <- rbind(linesX, linesY)
sf::st_crs(fishnet.multi) <- sf::st_crs(fishnet)
if(!all(sf::st_is_valid(fishnet.multi)))
{
warning("Fishnet contains invalid geometries!")
}
return(fishnet.multi)
} # end of function st_make_grid_lines
ggRaver/Lslide documentation built on April 8, 2022, 7:14 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions st_erase
Documented in st_erase
#' Erase one geometry from another
#'
#' This function erase one geometry from another. The projection must be identical.
#'
#' @param x object of class sf
#' @param y object of class sf
#' @importFrom magrittr "%>%"
#' @return
#' Geometry of class sfc
#'
#'
#' @keywords simple feature, erase
#'
#'
#' @export
#'
st_erase = function(x, y, precision = 0, do.subset = TRUE)
{
if(do.subset)
{
inter <- sf::st_intersects(x = y, y = x) %>% unlist(.) %>% unique(.)
x.remain <- x[-inter,]
x <- x[inter,]
}
if(precision != 0)
{
x <- x %>% sf::st_set_precision(x = ., precision = precision) %>% lwgeom::st_make_valid(.)
y <- y %>% sf::st_combine(.) %>% sf::st_union(.) %>%
sf::st_set_precision(x = ., precision = precision) %>% lwgeom::st_make_valid(.)
}
out <- sf::st_difference(x = x, y = y) # erase y from x
if(do.subset)
{
out <- rbind(out, x.remain)
}
return(out)
}
#' Erase one geometry from another using Saga GIS
#'
#' This function erase one geometry from another. The projection must be identical.
#'
#' @param x object of class sf. First element: Should be either line or polygon
#' @param y object of class sf. Second element: Always polygon.
#' @param method method of erase. Either "1": , or "2": Line-Polygon Intersection. Default: "1"
#' @param split Set to "1", if multi-part polygons should be splitted to single-part polygons. Default: "0"
#' @param attributes attributes inherited to intersection result. [0] polygon, [1] line, [2] line and polygon. Default: "1"
#' @param env.rsaga SAGA GIS environemnt. Default: RSAGA::rsaga.env()
#' @param check.geom If set to TRUE then geometry is checked with sf::st_is_valid(). If there are invalid geometries, geometries are repaired using lwgeom::st_make_valid(). Default: TRUE
#' @importFrom magrittr "%>%"
#' @return
#' Geometry of class sfc
#'
#'
#' @keywords simple feature, erase
#'
#'
#' @export
#'
rsaga_erase = function(x, y, method = "1", split = "0", attributes = "1", env.rsaga = RSAGA::rsaga.env(), check.geom = TRUE, quiet = TRUE)
{
path.x <- file.path(tempdir(), "tmp_x.shp")
path.y <- file.path(tempdir(), "tmp_y.shp")
path.result <- file.path(tempdir(), "tmp_result.shp")
sf::st_write(obj = x, dsn = path.x, delete_layer = TRUE, quiet = quiet)
sf::st_write(obj = y, dsn = path.y, delete_layer = TRUE, quiet = quiet)
if(method == "1")
{
# RSAGA::rsaga.get.usage(lib = "shapes_polygons", module = 15, env = env.rsaga)
RSAGA::rsaga.geoprocessor(lib = "shapes_polygons", module = 15, env = env.rsaga, show.output.on.console = !quiet, param = list(
A = path.x, B = path.y, RESULT = path.result, SPLIT = split))
}
if(method == "2")
{
# RSAGA::rsaga.get.usage(lib = "shapes_lines", module = 3, env = env.rsaga)
# ATTRIBUTES: [1] line
RSAGA::rsaga.geoprocessor(lib = "shapes_lines", module = 3, env = env.rsaga, show.output.on.console = !quiet, param = list(
LINES = path.x, POLYGONS = path.y, ATTRIBUTES = attributes, DIFFERENCE = path.result))
}
## read data
out <- sf::st_read(dsn = path.result, quiet = quiet)
## check validity
if(check.geom && !all(sf::st_is_valid(out)))
{
warning('Some invalid geometries by "rsaga_erase". Try to correct geomeries using lwgeom::st_make_valid()!')
out <- lwgeom::st_make_valid(x = out)
if(method == "1")
{
out <- suppressWarnings(out %>% sf::st_collection_extract(x = ., type = c("POLYGON")))
}
if(method == "2")
{
out <- suppressWarnings(out %>% sf::st_collection_extract(x = ., type = c("LINESTRING")))
}
}
return(out)
}
#' Return geometry from bounding box
#'
#' This function calculate the geometry based on a bounding box.
#'
#' @param x object of class sf
#' @param extent vector containing numeric values, in the following order: xmin, xmax, ymax, ymin
#' @importFrom magrittr "%>%"
#' @return
#' Geometry of class sfc
#'
#'
#' @keywords simple feature, geometry of bounding box
#'
#'
#' @export
#'
st_bbox_geom = function(x, extent = NULL)
{
if(!is.null(x) & !is.null(extent))
{
stop("Input conflict: either x or extent!")
}
if(is.null(extent))
{
out <- sf::st_bbox(x) %>% sf::st_as_sfc(.)
} else {
names(extent) <- c("xmin", "xmax", "ymax", "ymin")
out <- sf::st_bbox(extent) %>% sf::st_as_sfc(.)
}
return(out)
} # end of function st_bbox_geom
#' Plot sf geometry
#'
#' This function simply plots the geometry of a sf object.
#'
#' @param x object of class sf
#' @return
#' plot
#'
#'
#' @keywords simple feature, plot
#'
#'
#' @export
#'
st_plot = function(x) plot(sf::st_geometry(x))
#' Dissolve geometry
#'
#' This function dissovle a geometry based on a field. Optionnally, field statistics can be computed.
#'
#' @param x object of class sf
#' @param by field for dissolving. Default: NULL
#' @param ... Optional: field statistcs
#' @importFrom magrittr "%>%"
#' @return
#' dissolved geometry of class sf
#'
#'
#' @keywords simple feature, dissolve
#'
#'
#' @export
#'
st_dissolve = function(x, by = NULL, ...) x %>% dplyr::group_by(.dots = by) %>% dplyr::summarise(...)
#' Row bind simple features
#'
#' This function combines simple feature with different columns using dplr::bind_rows().
#'
#' @param x list of objects of class sf
#' @param geom vector containing the geometries of x. Order must be the same as in list.
#' @importFrom magrittr "%>%"
#' @return
#' sf object with combined data.frames
#'
#'
#' @keywords simple feature, rbind
#'
#'
#' @export
#'
st_rbind = function(x, geom) x %>% lapply(X = ., FUN = function(x) sf::st_set_geometry(x = x, value = NULL)) %>%
dplyr::bind_rows() %>% sf::st_set_geometry(x = ., value = geom)
#' Calculate the perimeter of a polygon
#'
#' This function calculates the perimeter of a polygon.
#'
#' @param x object of class sf
#' @importFrom magrittr "%>%"
#' @return
#' Vector with perimeter values
#'
#'
#' @keywords simple feature, perimeter
#'
#'
#' @export
#'
st_perimeter = function(x, ...) suppressWarnings(x %>% sf::st_cast(x = ., to = "MULTILINESTRING", ...) %>% sf::st_length(.) %>% as.numeric(.)) # units::drop_units(.)
#' Calculate the shape indices of a polygon
#'
#' This function calculates the shape index of a polygon.
#'
#' @param x object of class sf
#' @importFrom magrittr "%>%"
#' @return
#' Vector with shape index and interior edge ratio values
#'
#'
#' @keywords simple feature, shape index, interior edge ratio
#'
#'
#' @export
#'
st_shape_indices = function(x){
perimeter <- x %>% st_perimeter(.)
area <- x %>% sf::st_area(.) %>% as.numeric(.) # units::drop_units(.)
shape_index <- perimeter / (2 * sqrt(pi * area))
interior_edge_ratio <- perimeter / area
return(list(shape_index = shape_index, interior_edge_ratio = interior_edge_ratio))
}
#' Calculate integration index
#'
#' This function calculates an index defining the intergration of new land use area into conisting land use area
#' of the same land use class conisdering two different time steps.
#'
#' @param geom.old object of class sf representing a land use class of a previous time step
#' @param geom.new object of class sf representing a land use class of a following time step
#' @param geom.boundary polygon of class sf representing subregions, e.g. administrative boundaries
#' @param tol tolerance value for overlapping area m square
#' @param precision precision for process. See sf::st_set_precision(). Default: 0
#' @param env.rsaga environment of SAGA GIS. If st_erase fails then SAGA GIS erase is used. Default: NULL, but in function call if not set: RSAGA::rsaga.env()
#' @param use.saga use SAGA GIS for erase process. Default: FALSE
#' @param return.geom If set to TRUE, intermediate geometries are returned as well. Default: FALSE
#' @param quiet show output on console. Default: FALSE
#' @importFrom magrittr "%>%"
#' @return
#' Vector with integration index (completly integrated: 2/3 < R < 1; good integrated: 1/3 < R < 2/3; low integrated: 0 < R < 1/3; not integrated: 0)
#'
#'
#' @keywords simple feature, shape index, interior edge ratio
#'
#'
#' @export
#'
st_integration_index = function(geom.old, geom.new, geom.boundary = NULL, tol = 0.1, precision = 0,
ignr.overlap = FALSE, env.rsaga = NULL, use.saga = FALSE, return.geom = FALSE, quiet = FALSE){
# get start time of process
process.time.start <- proc.time()
## check input
if(missing(geom.old) || missing(geom.new)){ stop('Input is missing!')}
## check and set precision
if(!quiet) cat("... check and set st_precision \n")
if(sf::st_precision(geom.old) != precision)
{
geom.old <- geom.old %>% sf::st_set_precision(x = ., precision = precision) %>%
lwgeom::st_make_valid(.) %>%
sf::st_collection_extract(x = ., type = "POLYGON")
}
if(sf::st_precision(geom.new) != precision)
{
geom.new <- geom.new %>% sf::st_set_precision(x = ., precision = precision) %>%
lwgeom::st_make_valid(.) %>%
sf::st_collection_extract(x = ., type = "POLYGON")
}
## check validity of geometries
if(!all(sf::st_is_valid(geom.old))){ stop('Input of "geom.old" contains not valid geometries. Please try lwgeom::st_make_valid().')}
if(!all(sf::st_is_valid(geom.new))){ stop('Input of "geom.new" contains not valid geometries. Please try lwgeom::st_make_valid().')}
if(!is.null(geom.boundary) && !all(sf::st_is_valid(geom.boundary))){ stop('Input of "geom.boundary" contains not valid geometries. Please try lwgeom::st_make_valid().')}
# if(!quiet) cat("... union input geometries \n")
# geom.old <- sf::st_union(x = geom.old) %>% sf::st_cast(., "POLYGON") %>% sf::st_set_precision(x = ., precision = precision)
# geom.new <- sf::st_union(x = geom.new) %>% sf::st_cast(., "POLYGON") %>% sf::st_set_precision(x = ., precision = precision)
# # # # START CALCULATION OF INTEGRATION INDEX
## check for overlapping polygon
## common area
if(!quiet) cat("... intersection of input geometries \n")
inter <- suppressWarnings(sf::st_intersection(x = geom.old, y = geom.new) %>%
sf::st_collection_extract(x = ., type = c("POLYGON")))# %>%
# sf::st_set_precision(x = ., precision = precision)
## new area
if(!quiet) cat('... erase intersection from "geom.new" (this can take a while!) \n')
if(use.saga)
{
if(is.null(env.rsaga))
{
env.rsaga <- RSAGA::rsaga.env()
}
erase <- Lslide::rsaga_erase(x = geom.new, y = inter, method = "1",
split = "1", env.rsaga = env.rsaga) %>%
.[which(x = as.numeric(sf::st_area(.)) >= tol),]
} else {
erase <- tryCatch({
suppressWarnings(Lslide::st_erase(x = geom.new, y = inter, precision = precision) %>%
sf::st_collection_extract(x = ., type = c("POLYGON")) %>%
sf::st_cast(x = ., to = "POLYGON") %>%
.[which(x = as.numeric(sf::st_area(.)) >= tol),])
}, error = function(e){
warning(paste('SAGA GIS is used due to error in st_erase():', e))
if(is.null(env.rsaga))
{
env.rsaga <- RSAGA::rsaga.env()
}
tmp.erase <- Lslide::rsaga_erase(x = geom.new, y = inter, method = "1",
split = "1", env.rsaga = env.rsaga) %>%
.[which(x = as.numeric(sf::st_area(.)) >= tol),]
return(tmp.erase)
})
} # end of use.saga
if(!quiet) cat('... conversion to lines \n')
line.erase <- sf::st_cast(x = erase, to = "MULTILINESTRING") %>% sf::st_set_precision(x = ., precision = precision)
line.inter <- sf::st_cast(x = inter, to = "MULTILINESTRING") %>% sf::st_set_precision(x = ., precision = precision)
## common border
if(!quiet) cat('... find border lines by intersection \n')
unique.border <- suppressWarnings(sf::st_intersection(x = line.inter, y = line.erase) %>%
sf::st_collection_extract(x = ., type = c("LINESTRING")))
## check boundary constraints
if(!is.null(geom.boundary))
{
geom.boundary$ID_BOUNDS <- 1:nrow(geom.boundary)
geom.boundary <- geom.boundary[, c("ID_BOUNDS", "geometry")]
if(!quiet) cat('... intersection with boundaries \n')
unique.border <- suppressWarnings(sf::st_intersection(x = geom.boundary, y = unique.border) %>%
sf::st_collection_extract(x = ., type = c("LINESTRING")) %>%
st_cast(., "MULTILINESTRING") %>% sf::st_cast(., "LINESTRING")) # cast is necessairy to split multi-object
dt.unique.border <- unique.border %>%
sf::st_set_geometry(x = ., value = NULL) %>%
data.table::as.data.table(.)
erase <- suppressWarnings(sf::st_intersection(x = geom.boundary, y = erase) %>%
sf::st_collection_extract(x = ., type = c("POLYGON")) %>%
st_cast(., "MULTIPOLYGON") %>% sf::st_cast(., "POLYGON"))
dt.erase <- erase %>%
sf::st_set_geometry(x = ., value = NULL) %>%
data.table::as.data.table(.)
if(!quiet) cat('... get statistics and calculate index \n')
dt.unique.border$L <- sf::st_length(x = unique.border) %>% as.numeric() # units::drop_units(x = .)
dt.erase$P <- Lslide::st_perimeter(x = erase)
dt.result.UB <- dt.unique.border[,list(L = sum(L, na.rm = TRUE)), by = ID_BOUNDS]
dt.result.E <- dt.erase[,list(P = sum(P, na.rm = TRUE)), by = ID_BOUNDS]
dt.result <- merge(x = dt.result.E, y = dt.result.UB, by = "ID_BOUNDS") %>%
dplyr::mutate(.data = ., R = L/P)
} else {
if(!quiet) cat('... get statistics and calculate index \n')
dt.result <- data.table::data.table(P = sum(Lslide::st_perimeter(x = erase), na.rm = TRUE),
L = sum((sf::st_length(x = unique.border) %>% as.numeric(.)), na.rm = TRUE)) %>% # units::drop_units(x = .))) %>%
dplyr::mutate(.data = ., R = L/P)
} # end of if-else boundary check
process.time.run <- proc.time() - process.time.start
if(quiet == FALSE) cat("------ Run of st_integration_index: " , round(process.time.run["elapsed"][[1]]/60, digits = 4), " Minutes \n")
if(return.geom)
{
return(list(integration_index = dt.result, unique_border = unique.border, new_area = erase))
} else {
return(dt.result)
}
} # end of function st_integration_index
#' Calculate the effective mesh size
#'
#' This function calculates the effective mesh size.
#'
#' @param geom.frag polygon of class sf representing the fragmentation geometry
#' @param geom.boundary polygon of class sf representing subregions, e.g. administrative boundaries
#' @param total.area Numeric value representing size for area. Only to use if geom.boundary is not present. Value must match with the conversion constant c (default hectare). Default: NULL
#' @param conv constant to convert original square m output. Default: 10000 to convert to hectare. If set to 1, than meter square is the result.
#' @param do.preProcessing If TRUE (default), the input of geom.frag is, first, dissolved to single part feature, and second, splitted to multi-parts. By this step it is assured, that polygon connected to each other are summarized
#' @param return.geom If set to TRUE, intermediate geometries are returned as well. Default: FALSE
#' @param quiet If set to FALSE, actual state is printed to console. Default: TRUE.
#' @importFrom magrittr "%>%"
#' @return
#' Vector with shape index and interior edge ratio values
#'
#'
#' @keywords simple feature, mesh, effective mesh size
#'
#'
#' @export
#'
st_mesh = function(geom.frag, geom.boundary = NULL, total.area = NULL, conv = 10000, do.preProcessing = TRUE, return.geom = FALSE, quiet = TRUE)
{
# get start time of process
process.time.start <- proc.time()
## check input
if(missing(geom.frag)){ stop('Input of "geom.frag" is missing!')}
if(!missing(geom.frag) && !("sf" %in% class(geom.frag))){ stop('Input of "geom.frag" is not of class "sf"!')}
if(!is.null(geom.boundary) && !("sf" %in% class(geom.boundary))){ stop('Input of "geom.boundary" is not of class "sf"!')}
if(is.null(geom.boundary) && (is.null(total.area) || !is.numeric(total.area))){
stop('If input of "geom.boundary" is null, then input of "total.area" must be a numeric value representing the size of the study area.')
}
## check validity of geometries
if(!all(sf::st_is_valid(geom.frag))){ stop('Input of "geom.frag" contains not valid geometries. Please try lwgeom::st_make_valid().')}
if(!is.null(geom.boundary) && !all(sf::st_is_valid(geom.boundary))){ stop('Input of "geom.boundary" contains not valid geometries. Please try lwgeom::st_make_valid().')}
if(do.preProcessing)
{
if(!quiet) cat("... union geometries to a single geometry with resolved boundaries \n")
geom.frag <- geom.frag %>% sf::st_union(.)
if(!quiet) cat("... split multi-parts to single-parts polygon \n")
geom.frag <- geom.frag %>% st_cast(., "POLYGON") %>% sf::st_sf(ID_FRAG = 1:length(.), geometry = .)
} else {
geom.frag$ID_FRAG <- 1:nrow(geom.frag) ## add unique IDs
}
if(!is.null(geom.boundary)){ geom.boundary$ID_BOUNDS <- 1:nrow(geom.boundary) } # ## add unique IDs
## add Area in m_sq
geom.frag$A_FRAG <- sf::st_area(geom.frag) %>% (function(x, conv = conv) as.numeric(x)/conv) # as.numeric() %>% "/" (conv) # units::drop_units(.)
if(!is.null(geom.boundary)){ geom.boundary$A_BOUNDS <- sf::st_area(geom.boundary) %>% (function(x, conv = conv) as.numeric(x)/conv)}
## subset data
geom.frag <- geom.frag[, c("ID_FRAG", "A_FRAG", "geometry")]
if(!is.null(geom.boundary)){geom.boundary <- geom.boundary[, c("ID_BOUNDS", "A_BOUNDS", "geometry")]}
# # # CALCULATE MESH SIZE INDICES
## start calculation
if(!is.null(geom.boundary))
{
## get intersection
if(!quiet) cat("... intersection to boundary \n")
inter <- suppressWarnings(sf::st_intersection(x = geom.boundary, y = geom.frag))
inter <- suppressWarnings(inter %>% st_cast(., "MULTIPOLYGON") %>% sf::st_cast(., "POLYGON")) # cast is necessairy to split multi-polygons
inter$A_FRAG_INTER <- sf::st_area(inter) %>% (function(x, conv = conv) as.numeric(x)/conv) # overwrite area of fragments, # units::drop_units(.)
## calculation of mesh indices
df.inter <- sf::st_set_geometry(x = inter, value = NULL) %>% data.table::as.data.table(.)
df.inter.multi <- df.inter[, list(A_BOUNDS = unique(A_BOUNDS), # splitted multi-parts are joined together
A_FRAG = unique(A_FRAG),
A_FRAG_INTER = sum(A_FRAG_INTER, na.rm = TRUE)), by = list(ID_BOUNDS, ID_FRAG)]
if(!quiet) cat("... get statistics \n")
## cutting-out (CUT) procedure
mesh.CUT <- df.inter[, list(Fg = unique(A_BOUNDS),
Fi = sum(A_FRAG_INTER^2, na.rm = TRUE),
Fi_count = length(!is.na(A_FRAG_INTER)),
Fi_sum = sum(A_FRAG_INTER, na.rm = TRUE),
Fi_min = min(A_FRAG_INTER, na.rm = TRUE),
Fi_max = max(A_FRAG_INTER, na.rm = TRUE),
Fi_ave = mean(A_FRAG_INTER, na.rm = TRUE),
Ci = sum((A_FRAG_INTER/A_BOUNDS)^2, na.rm = TRUE)), by = ID_BOUNDS] %>%
dplyr::mutate(., Di = 1-Ci) %>%
dplyr::mutate(., Si = 1/Ci) %>%
dplyr::mutate(., mEff_CUT = Fi/Fg)
# View(mesh.CUT)
## cross-boundary connections (CBC) procedure
mesh.CBC <- df.inter.multi[, list(Fg = unique(A_BOUNDS),
FiErg_count = length(!is.na(A_FRAG)),
FiErg_sum = sum(A_FRAG, na.rm = TRUE),
FiErg_min = min(A_FRAG, na.rm = TRUE),
FiErg_max = max(A_FRAG, na.rm = TRUE),
FiErg_ave = mean(A_FRAG, na.rm = TRUE),
CiErg = sum((A_FRAG/A_BOUNDS)^2, na.rm = TRUE),
Fi_CBC1 = sum(A_FRAG_INTER*A_FRAG, na.rm = TRUE),
Fi_CBC2 = sum(2*A_FRAG_INTER*A_FRAG-A_FRAG_INTER^2, na.rm = TRUE)), by = ID_BOUNDS] %>%
dplyr::mutate(., DiErg = 1-CiErg) %>%
dplyr::mutate(., SiErg = 1/CiErg) %>%
dplyr::mutate(., mEff_CBC1 = Fi_CBC1/Fg) %>%
dplyr::mutate(., mEff_CBC2 = Fi_CBC2/Fg)
# View(mesh.CBC)
df.result <- merge(x = mesh.CUT, y = mesh.CBC[, -2], by = "ID_BOUNDS")
} else {
if(!quiet) cat("... get statistics \n")
df.result <- sf::st_set_geometry(x = geom.frag, value = NULL) %>% data.table::as.data.table(.)
df.result <- df.result[, list(Fg = total.area,
Fi = sum(A_FRAG^2, na.rm = TRUE),
Fi_count = length(!is.na(A_FRAG)),
Fi_sum = sum(A_FRAG, na.rm = TRUE),
Fi_min = min(A_FRAG, na.rm = TRUE),
Fi_max = max(A_FRAG, na.rm = TRUE),
Fi_ave = mean(A_FRAG, na.rm = TRUE),
Ci = sum((A_FRAG/total.area)^2, na.rm = TRUE)),] %>%
dplyr::mutate(., Di = 1-Ci) %>%
dplyr::mutate(., Si = 1/Ci) %>%
dplyr::mutate(., mEff_CUT = Fi/total.area)
}
process.time.run <- proc.time() - process.time.start
if(quiet == FALSE) cat("------ Run of st_mesh: " , round(process.time.run["elapsed"][[1]]/60, digits = 4), " Minutes \n")
if(return.geom)
{
if(is.null(geom.boundary))
{
return(list(mesh = df.result, geom.frag = geom.frag))
} else {
return(list(mesh = df.result, geom.frag = geom.frag, geom.inter = inter))
}
} else {
return(df.result)
}
} # end of function st_mesh
#' Calculate urban sprawl
#'
#' This function calculates the urban sprawl of urban area.
#'
#' @param geom.urban polygon of class sf representing the fragmentation geometry
#' @param geom.boundary polygon of class sf representing subregions, e.g. administrative boundaries
#' @param dist vector containing distance between lines in x and y direction. Default: c(100, 100) [m]
#' @param trans transformation function {x-1+1/(x+trans.k)} with x as free line and trans.k a constant
#' @param trans.k constant in km for transformation function trans. Default: 1
#' @param tol tolerance value for intersection with erased lines. Buffering procedure is used.Default: 0.1 [m]
#' @param precision precision for process. See sf::st_set_precision(). Default: 0
#' @param extent Numeric value representing extent for area. Format of vector: c(xmin, xmax, ymax, ymin) . Default: NULL
#' @param force.extent If TRUE extent is used instead of geom.boundary (if both are present). Default: FALSE
#' @param do.preProcessing If TRUE (default), the input of geom.frag is, first, dissolved to single part feature, and second, splitted to multi-parts. By this step it is assured, that polygon connected to each other are summarized
#' @param return.geom If set to TRUE, intermediate geometries are returned as well. Default: FALSE
#' @param env.rsaga environment of SAGA GIS. If st_erase fails then SAGA GIS erase is used. Default: NULL, but in function call if not set: RSAGA::rsaga.env()
#' @param use.saga use SAGA GIS for erase process. Default: FALSE
#' @param quiet If set to FALSE, actual state is printed to console. Default: TRUE.
#' @importFrom magrittr "%>%"
#' @return
#' strong urban sprawl: 40-50%, less urban sprawl: 80-90%
#'
#'
#' @keywords simple feature, urban sprawl
#'
#'
#' @export
#'
st_urban_sprawl = function(geom.urban, geom.boundary = NULL, dist = c(100, 100), trans = function(x, trans.k){x-1+1/(x+trans.k)}, trans.k = 1, tol = 0.1, extent = NULL, force.extent = FALSE,
precision = 0, do.preProcessing = TRUE, env.rsaga = NULL, use.saga = FALSE, return.geom = FALSE, quiet = TRUE)
{
# get start time of process
process.time.start <- proc.time()
## check input
if(missing(geom.urban)){ stop('Input of "geom.urban" is missing!')}
if(!missing(geom.urban) && !("sf" %in% class(geom.urban))){ stop('Input of "geom.urban" is not of class "sf"!')}
if(!is.null(geom.boundary) && !("sf" %in% class(geom.boundary))){ stop('Input of "geom.boundary" is not of class "sf"!')}
if(is.null(geom.boundary) && is.null(extent)){
warning('If input of "geom.boundary" and "extent" is null. Fish net is created using the extent of "geom.urban"!')
}
if(is.null(geom.boundary) & is.null(extent) & force.extent){
stop('If "force.extent" is TRUE, than extent should be set!"!')
}
if(!is.null(geom.boundary) & !is.null(extent)){
warning('Fish net is created using extent of "geom.boundary". "Extent" is skipped. For "extent" use "force.extent"!')
}
if(sf::st_precision(geom.urban) != precision)
{
geom.urban <- geom.urban %>% sf::st_set_precision(x = ., precision = precision) %>%
lwgeom::st_make_valid(.) %>%
sf::st_collection_extract(x = ., type = "POLYGON")
}
## check validity of geometries
if(!all(sf::st_is_valid(geom.urban))){ stop('Input of "geom.urban" contains not valid geometries. Please try lwgeom::st_make_valid().')}
if(!is.null(geom.boundary) && !all(sf::st_is_valid(geom.boundary))){ stop('Input of "geom.boundary" contains not valid geometries. Please try lwgeom::st_make_valid().')}
## add unique ID and subset data
if(!is.null(geom.boundary)){ geom.boundary$ID_BOUNDS <- 1:nrow(geom.boundary) } # ## add unique IDs
if(!is.null(geom.boundary)){geom.boundary <- geom.boundary[, c("ID_BOUNDS", "geometry")]}
## ... create boundary for fishnet
if(!is.null(geom.boundary) & !force.extent){bbox.fishnet <- Lslide::st_bbox_geom(x = geom.boundary)
} else if(!is.null(geom.boundary) & force.extent){bbox.fishnet <- Lslide::st_bbox_geom(extent = extent)
} else if(is.null(geom.boundary) & !is.null(extent)){bbox.fishnet <- Lslide::st_bbox_geom(extent = extent)
} else {bbox.fishnet <- Lslide::st_bbox_geom(x = geom.urban) }
## create and subset fish net
if(!quiet) cat("... create fishnet \n")
fishnet <- Lslide::st_make_grid_lines(x = bbox.fishnet, cellsize = dist)
if(!is.null(geom.boundary) & !force.extent){
if(!quiet) cat("... subset fishnet to area of interest \n")
fishnet <- suppressWarnings(sf::st_intersection(x = fishnet, y = geom.boundary)) %>%
sf::st_collection_extract(x = ., type = "LINESTRING", warn = FALSE)
} else {
fishnet <- suppressWarnings(sf::st_intersection(x = fishnet, y = bbox.fishnet)) %>%
sf::st_collection_extract(x = ., type = "LINESTRING", warn = FALSE)
}
# fishnet <- fishnet[, c("ID", "geometry")]
## do post-processing of urban geometry
if(do.preProcessing)
{
if(!quiet) cat("... union geometries to a single geometry with resolved boundaries \n")
geom.urban <- geom.urban %>% sf::st_union(.)
if(!quiet) cat("... split multi-parts to single-parts polygon \n")
geom.urban <- geom.urban %>% st_cast(., "POLYGON") %>% sf::st_sf(ID_URBAN = 1:length(.), geometry = .)
} else {
geom.urban$ID_URBAN <- 1:nrow(geom.urban) ## add unique IDs
}
## check intersection
inter.check <- sf::st_intersects(x = geom.urban, y = fishnet)
inter.empty <- which(sapply(inter.check, function(x) length(x) == 0))
if(length(inter.empty) > 0)
{
inter.A <- geom.urban[inter.empty,] %>% sf::st_area(.) %>% as.numeric(.)
inter.num <- length(inter.empty)
warning('Some urban polygons are not intersected by fishnet: ', inter.num,
' | area min: ', min(inter.A, na.rm = TRUE), ' - max: ', max(inter.A, na.rm = TRUE), ' [m_sqr] \n')
}
## erase urban area from fishnet using SAGA GIS
if(!quiet) cat("... erase urban area from fishnet \n")
if(use.saga)
{
if(is.null(env.rsaga))
{
env.rsaga <- RSAGA::rsaga.env()
}
erase <- Lslide::rsaga_erase(x = fishnet, y = geom.urban, method = "2", env.rsaga = env.rsaga)
} else {
erase <- suppressWarnings(Lslide::st_erase(x = fishnet, y = geom.urban, precision = precision) %>%
sf::st_collection_extract(x = ., type = "LINESTRING"))
}
## split to single part
erase.single <- erase %>% sf::st_cast(x = ., to = 'LINESTRING', warn = FALSE)
## selection of lines
if(!quiet) cat("... selection of lines \n")
lines.urban <- sf::st_intersects(x = geom.urban %>%
sf::st_buffer(x = ., dist = tol),
y = erase.single) %>%
unlist(.) %>% unique(.)
if(!is.null(geom.boundary) & !force.extent)
{
lines.boundary <- sf::st_intersects(x = geom.boundary %>%
sf::st_boundary(x = .) %>%
sf::st_buffer(x = ., dist = tol),
y = erase.single) %>%
unlist(.) %>% unique(.)
} else {
lines.boundary <- sf::st_intersects(x = bbox.fishnet %>%
sf::st_boundary(x = .) %>%
sf::st_buffer(x = ., dist = tol),
y = erase.single) %>%
unlist(.) %>% unique(.)
}
# lines.urban %>% length(.)
# lines.boundary %>% length(.)
check.missing <- setdiff(c(1:nrow(erase.single)),
(c(lines.urban, lines.boundary) %>% unique(.)))
if(length(check.missing) > 0)
{
warning("Some lines are neither intersected by boundary nor by urban area: ", length(check.missing), " in total \n")
}
erase.single$Type <- NA
# TPYE 1: Lines between urban area (red lines)
type.urban <- setdiff(x = lines.urban, y = lines.boundary)
erase.single$Type[type.urban] <- 1
# TPYE 2: Lines not touching urban area (blue lines)
type.boundary <- setdiff(x = lines.boundary, y = lines.urban)
erase.single$Type[type.boundary] <- 2
# TYPE 3: Lines between boundary and urban area (green lines)
type.bound.urb <- intersect(x = lines.urban, y = lines.boundary)
erase.single$Type[type.bound.urb] <- 3
# summary(erase.single$Type)
# c(type.bound.urb, type.boundary, type.urban) %>% unique(.) %>% length(.)
# browser()
## group single lines back to multi-lines
if((!is.null(geom.boundary) & !force.extent))
{
erase.final <- st_dissolve(x = erase.single[which(erase.single$Type == 2 | erase.single$Type == 3),],
by = list("Type", "ID_FNET", "ID_BOUNDS"))
erase.final <- rbind(erase.final, erase.single[which(erase.single$Type == 1),])
} else {
erase.final <- st_dissolve(x = erase.single[which(erase.single$Type == 2 | erase.single$Type == 3),],
by = list("Type", "ID_FNET"))
erase.final <- rbind(erase.final, erase.single[which(erase.single$Type == 1),])
}
erase.final <- erase.final %>% sf::st_collection_extract(x = ., type = "LINESTRING")
erase.final$L <- sf::st_length(x = erase.final) %>% (function(x) as.numeric(x)/1000) # as.numeric(.) %>% "/" (1000) # convert from meter to km!
erase.final$L_trans <- ifelse(erase.final$Type == 2, erase.final$L, trans(x = erase.final$L, trans.k = trans.k))
## summarizing final data
if(!quiet) cat("... get statistics: urban sprawl \n")
if((!is.null(geom.boundary) & !force.extent))
{
df.result <- erase.final %>% sf::st_set_geometry(x = ., value = NULL) %>%
data.table::as.data.table(.) %>%
.[,list(L = sum(L, na.rm = TRUE),
L_trans = sum(L_trans, na.rm = TRUE)), by = ID_BOUNDS] %>%
dplyr::mutate(.data = ., FFE = L_trans/L*100)
} else {
df.result <- erase.final %>% sf::st_set_geometry(x = ., value = NULL) %>%
data.table::as.data.table(.) %>%
.[,list(L = sum(L, na.rm = TRUE),
L_trans = sum(L_trans, na.rm = TRUE)),] %>%
dplyr::mutate(.data = ., FFE = L_trans/L*100)
}
process.time.run <- proc.time() - process.time.start
if(quiet == FALSE) cat("------ Run of urban_sprawl: " , round(process.time.run["elapsed"][[1]]/60, digits = 4), " Minutes \n")
if(return.geom)
{
return(list(FFE = df.result, geom = erase.final))
} else {
return(df.result)
}
} # end of function st_urban_sprawl
#' Make multiline fishnet
#'
#' This function creates a fishnet, consisting of multilines in x and y direction.
#'
#' @param x object of class sf
#' @importFrom magrittr "%>%"
#' @return
#' Geometry of class sfc
#'
#'
#' @keywords simple feature, fishnet, multilines
#'
#'
#' @export
#'
st_make_grid_lines = function(x, cellsize)
{
# ... create fishnet using corners (lower left corner is starting point)
fishnet <- sf::st_make_grid(x = x, cellsize = cellsize, what = "corners")
# ... get sf::st_make_grid function parameters (automatically calculated in function call)
offset <- sf::st_bbox(obj = x)[1:2]
nx <- ceiling((sf::st_bbox(obj = x)[3] - offset[1])/cellsize[1]) + 1 # add 1 to fit overall columns
ny <- ceiling((sf::st_bbox(obj = x)[4] - offset[2])/cellsize[2]) + 1 # add 1 to fit overall rows
# ... create lines
linesX <- lapply(1:ny, function(i, nx, fishnet){
if(i == 1)
{
index <- 1:nx
} else {
index <- ((i-1)*nx+1):(i*nx)
}
line <- fishnet[index] %>% st_combine(.) %>%
sf::st_multilinestring(x = ., dim = "XY") %>%
st_sfc(.) %>%
sf::st_sf(ID_FNET = paste0("X", i), geometry = .)
return(line)
}, nx = nx, fishnet = fishnet) %>% do.call(what = rbind, args = .)
seqY <- seq(from = 1, to = length(fishnet), by = nx)
linesY <- lapply(0:(nx-1), FUN = function(i, seqY, fishnet){
if(i == 0){
index <- seqY
} else {
index <- seqY+i
}
line <- fishnet[index] %>% st_combine(.) %>%
sf::st_multilinestring(x = ., dim = "XY") %>%
st_sfc(.) %>%
sf::st_sf(ID_FNET = paste0("Y", i), geometry = .)
return(line)
}, seqY = seqY, fishnet) %>% do.call(what = rbind, args = .)
fishnet.multi <- rbind(linesX, linesY)
sf::st_crs(fishnet.multi) <- sf::st_crs(fishnet)
if(!all(sf::st_is_valid(fishnet.multi)))
{
warning("Fishnet contains invalid geometries!")
}
return(fishnet.multi)
} # end of function st_make_grid_lines
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.