Skip to content

replaced loop in windfetch() with a parallel backend loop from foreac… #2

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
anguswg-ucsb wants to merge 3 commits into
base: main
Choose a base branch
from
Open

replaced loop in windfetch() with a parallel backend loop from foreac… #2

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
6d81716
replaced loop in windfetch() with a parallel backend loop from foreac…
anguswg-ucsb Sep 16, 2022
027769d
as_sf as a parameter
anguswg-ucsb Sep 16, 2022
3a0e3b7
error searching
anguswg-ucsb Sep 17, 2022
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
5 changes: 4 additions & 1 deletion DESCRIPTION
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -11,7 +11,10 @@ Imports:
magrittr,
sf,
utils,
rlang
rlang,
foreach,
doParallel,
parallel
URL: https://github.com/blasee/windfetch
BugReports: https://github.com/blasee/windfetch/issues
Description: Wind fetch is the unobstructed length of water over which wind can
Expand Down
9 changes: 9 additions & 0 deletions NAMESPACE
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -6,13 +6,17 @@ exportMethods(as_sf)
exportMethods(crs_transform)
exportMethods(plot)
exportMethods(summary)
importFrom(doParallel,registerDoParallel)
importFrom(dplyr,as_tibble)
importFrom(dplyr,bind_cols)
importFrom(dplyr,group_by)
importFrom(dplyr,left_join)
importFrom(dplyr,mutate)
importFrom(dplyr,relocate)
importFrom(dplyr,select)
importFrom(dplyr,summarise)
importFrom(dplyr,tibble)
importFrom(foreach,foreach)
importFrom(magrittr,"%>%")
importFrom(methods,"S3Part<-")
importFrom(methods,as)
Expand All @@ -25,13 +29,17 @@ importFrom(methods,setMethod)
importFrom(methods,setValidity)
importFrom(methods,show)
importFrom(methods,validObject)
importFrom(parallel,detectCores)
importFrom(parallel,makeCluster)
importFrom(parallel,stopCluster)
importFrom(rlang,.data)
importFrom(sf,"st_geometry<-")
importFrom(sf,st_as_sf)
importFrom(sf,st_buffer)
importFrom(sf,st_cast)
importFrom(sf,st_coordinates)
importFrom(sf,st_crs)
importFrom(sf,st_drop_geometry)
importFrom(sf,st_intersection)
importFrom(sf,st_intersects)
importFrom(sf,st_is)
Expand All @@ -44,6 +52,7 @@ importFrom(sf,st_point)
importFrom(sf,st_sf)
importFrom(sf,st_sfc)
importFrom(sf,st_transform)
importFrom(stats,setNames)
importFrom(utils,head)
importFrom(utils,setTxtProgressBar)
importFrom(utils,txtProgressBar)
7 changes: 4 additions & 3 deletions R/internals.R
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -3,14 +3,15 @@
#' @importFrom sf st_intersection st_cast st_nearest_points st_nearest_feature

return_fetch_vector = function(vector, origin, polygon) {

fetch_on_land = suppressWarnings(st_intersection(vector, polygon))

if (nrow(fetch_on_land) == 0) {
fetch = vector
} else {
fetch_ends = suppressWarnings(st_cast(fetch_on_land, "POINT"))
fetch_vector = st_nearest_points(origin, fetch_ends)
which_closest = st_nearest_feature(origin, fetch_ends)
fetch_ends = suppressWarnings(st_cast(fetch_on_land, "POINT"))
fetch_vector = st_nearest_points(origin, fetch_ends)
which_closest = st_nearest_feature(origin, fetch_ends)
fetch = fetch_vector[which_closest]
}
fetch
Expand Down
195 changes: 151 additions & 44 deletions R/windfetch.R
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -37,6 +37,8 @@
#' quadrant (default 9).
#' @param quiet logical. Suppress diagnostic messages? (Default \code{FALSE}).
#' @param progress_bar logical. Show a text progress bar? (Default \code{TRUE})
#' @param as_sf logical. Should an SF object be returned on only the fetch distances. Default is TRUE, returns SF object
#' @param ncores numeric. The number of cores on ones computer to use for parallel processing backend. Default is 2.
#'
#' @return Returns a \code{\link{WindFetch}} object.
#'
Expand All @@ -53,16 +55,25 @@
#'
#' @importFrom utils head txtProgressBar setTxtProgressBar
#' @importFrom methods new is
#' @importFrom sf st_is st_is_longlat st_crs st_transform st_intersects
#' st_buffer st_coordinates st_sf st_sfc st_linestring st_point
#' @importFrom sf st_is st_is_longlat st_crs st_transform st_intersects st_as_sf
#' st_buffer st_coordinates st_sf st_sfc st_linestring st_point st_drop_geometry
#' st_length st_geometry<- st_as_sf
#' @importFrom dplyr left_join
#' @importFrom dplyr left_join mutate relocate
#' @importFrom stats setNames
#' @importFrom parallel detectCores makeCluster stopCluster
#' @importFrom doParallel registerDoParallel
#' @importFrom foreach foreach
#'
#' @include windfetch_class.R
#' @examples
#'
#' # Create the polygon layer ----------------------------------------
#'
#' library(sf)
#' library(doParallel)
#' library(parallel)
#' library(foreach)
#' library(dplyr)
#'
#' # Read in North Carolina shape file
#' nc = st_read(system.file("shape/nc.shp", package = "sf"), quiet = TRUE)
Expand Down Expand Up @@ -133,9 +144,16 @@
#' st_write(nc_fetch_sf, "nc_fetch.shp", driver = "ESRI Shapefile")
#' }
#' @export
windfetch = function(polygon_layer, site_layer, max_dist = 300, n_directions = 9,
quiet = FALSE, progress_bar = TRUE) {

windfetch <- function(
polygon_layer,
site_layer,
max_dist = 300,
n_directions = 9,
quiet = FALSE,
progress_bar = TRUE,
Copy link
Owner

@blasee blasee Sep 26, 2022

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The progress_bar argument is no longer required. I think maybe we should only have a progress bar when 1 core is selected. Thoughts?

as_sf = TRUE,
Copy link
Owner

@blasee blasee Sep 26, 2022

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

as_sf is misleading as the object returned is actually a windFetch object, which doesn't contain the sf class. I think we get rid of this argument, unless you have a good reason why the windfetch function shouldn't act as a constructor function for the windFetch class.

ncores = 2
Copy link
Owner

@blasee blasee Sep 26, 2022

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think a good default is a single core.

) {
if (!any(is(polygon_layer, "sf"), is(polygon_layer, "sfc")))
stop(paste("polygon_layer must be either an sf or sfc object.\nSee",
"`?read_sf` for details on how to read in an",
Expand All @@ -154,8 +172,11 @@ windfetch = function(polygon_layer, site_layer, max_dist = 300, n_directions = 9
if (!is.numeric(max_dist) || length(max_dist) != 1)
stop("max_dist must be a single number.", call. = FALSE)

if (!is.numeric(n_directions) || length(n_directions) != 1)
if (!is.numeric(n_directions) || length(n_directions) != 1) {

stop("n_directions must be a single integer.", call. = FALSE)

}
n_directions = round(n_directions)

if (n_directions < 1 || n_directions > 90)
Expand All @@ -168,32 +189,42 @@ windfetch = function(polygon_layer, site_layer, max_dist = 300, n_directions = 9

which_proj = c(!st_is_longlat(polygon_layer), !st_is_longlat(site_layer))

if (all(!which_proj))
if (all(!which_proj)) {

stop("The polygon_layer or site_layer argument must be projected to calculate fetch.",
call. = FALSE)

if (all(which_proj) && !(st_crs(site_layer) == st_crs(polygon_layer))){
}

if (all(which_proj) && !(st_crs(site_layer) == st_crs(polygon_layer))) {
warning(paste("The CRS for the polygon_layer and site_layer arguments",
"differ; transforming site_layer CRS to match."),
call. = FALSE)

site_layer = st_transform(site_layer, st_crs(polygon_layer))

}

if (!which_proj[1]){
if (!quiet)
if (!quiet) {

message("projecting polygon_layer onto the site_layer CRS")

}
polygon_layer = st_transform(polygon_layer, st_crs(site_layer))
}

if (!which_proj[2]){
if (!quiet)

if (!quiet) {
message("projecting site_layer onto the polygon_layer CRS")
}
site_layer = st_transform(site_layer, st_crs(polygon_layer))
}

if (!quiet)
if (!quiet) {
message("checking site locations are not on land")
}

sites_on_land = st_intersects(site_layer, polygon_layer, sparse = FALSE)[, 1]

Expand All @@ -202,18 +233,27 @@ windfetch = function(polygon_layer, site_layer, max_dist = 300, n_directions = 9
call. = FALSE)
site_layer = site_layer[!sites_on_land, ]
}


# convert dataframes to SF
if (!is.data.frame(site_layer)) {

site_layer = st_as_sf(site_layer)

}

# name sites
if (any(grepl("^[Nn]ames{0,1}$", names(site_layer)))){
name_col = grep("^[Nn]ames{0,1}$", names(site_layer))
message("name section1")
Copy link
Owner

@blasee blasee Sep 26, 2022

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This is an uninformative message.

name_col = grep("^[Nn]ames{0,1}$", names(site_layer))

site_names = as.character(site_layer[[name_col]])

} else {
message("name section2")
Copy link
Owner

@blasee blasee Sep 26, 2022

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Another uninformative message.

site_names = paste("Site", seq_len(dim(site_layer)[[1]]))
}

}
message("name section3")
Copy link
Owner

@blasee blasee Sep 26, 2022

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I don't think we need any of these extra messages.

site_layer$site_names = site_names

# Convert max_dist to appropriate units.
Expand All @@ -223,12 +263,18 @@ windfetch = function(polygon_layer, site_layer, max_dist = 300, n_directions = 9
# Double check if metres are the correct units, and warn otherwise.
proj_unit = st_crs(polygon_layer)$units

if (proj_unit != "m")
if (proj_unit != "m") {

warning(paste("The CRS unit is not metres; ensure that the max_dist argument",
"has been scaled appropriately."), call. = FALSE)

}

message("dir section1")
Copy link
Owner

@blasee blasee Sep 26, 2022

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Remove.

directions = head(seq(0, 360, by = 360 / (n_directions * 4)), -1)


message("dir section2")
Copy link
Owner

@blasee blasee Sep 26, 2022

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Remove.

# Return the quadrant the directions belong to
dirs = as.numeric(directions)
dirs_bin = findInterval(dirs, seq(45, 315, by = 90))
Expand All @@ -237,60 +283,121 @@ windfetch = function(polygon_layer, site_layer, max_dist = 300, n_directions = 9
quadrant[dirs_bin == 2] = "South"
quadrant[dirs_bin == 3] = "West"

message("dir section3")
Copy link
Owner

@blasee blasee Sep 26, 2022

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Remove.


# Rearrange sequence order to start at 90 degrees
directions = unlist(split(directions, directions < 90), use.names = FALSE)

fetch_ends = st_buffer(site_layer, max_dist, n_directions)
fetch_ends_df = as.data.frame(st_coordinates(fetch_ends))

message("name section4")
Copy link
Owner

@blasee blasee Sep 26, 2022

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Remove :-)


fetch_ends_df$site_names = site_names[fetch_ends_df[, 4]]

fetch_locs_df = as.data.frame(st_coordinates(site_layer))

colnames(fetch_locs_df) = c("X0", "Y0")

message("name section5")
Copy link
Owner

@blasee blasee Sep 26, 2022

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Remove.


fetch_locs_df$site_names = site_names

fetch_df = unique(left_join(fetch_ends_df, fetch_locs_df, by = "site_names"))
fetch_df$directions = c(head(seq(90, 360, by = 360 / (n_directions * 4)), -1),
head(seq(0, 90, by = 360 / (n_directions * 4)), -1))
fetch_df = fetch_df[with(fetch_df, order(site_names, directions)), ]


message("nested section")
Copy link
Owner

@blasee blasee Sep 26, 2022

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Remove.

fetch_df = st_sf(fetch_df[, c("site_names", "directions")],
geom = st_sfc(sapply(apply(fetch_df, 1, function(x) {
X0 = as.numeric(x['X0'])
Y0 = as.numeric(x['Y0'])
X = as.numeric(x['X'])
Y = as.numeric(x['Y'])
st_sfc(st_linestring(matrix(c(X0, X, Y0, Y), 2, 2), dim = "XY"))
}), st_sfc), crs = st_crs(polygon_layer)),
origin = st_sfc(sapply(apply(fetch_df, 1, function(x) {
geom = st_sfc(
sapply(
apply(fetch_df, 1, function(x) {
X0 = as.numeric(x['X0'])
Y0 = as.numeric(x['Y0'])
X = as.numeric(x['X'])
Y = as.numeric(x['Y'])
st_sfc(
st_linestring(matrix(c(X0, X, Y0, Y), 2, 2), dim = "XY")
)
}),
st_sfc
),
crs = st_crs(polygon_layer)
),
origin = st_sfc(
sapply(apply(fetch_df, 1, function(x) {
st_sfc(st_point(c(as.numeric(x['X0']), as.numeric(x['Y0']))))
}), st_sfc), crs = st_crs(polygon_layer)))

message("Subset section")
Copy link
Owner

@blasee blasee Sep 26, 2022

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Remove.


poly_subset = subset(polygon_layer, lengths(st_intersects(polygon_layer, fetch_df)) > 0)
message("parralel")
Copy link
Owner

@blasee blasee Sep 26, 2022

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Remove.

# make a cluster
cl <- parallel::makeCluster(ncores) #not to overload your computer
Copy link
Owner

@blasee blasee Sep 26, 2022

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think we should include an error message if ncores > parallel::detectCores()


# register parallel backend
doParallel::registerDoParallel(cl)

# run a loop using a parallel backend. ~2x faster
fetch_loop <- foreach::foreach(i = 1:nrow(fetch_df),
.combine = "rbind",
.export = c('return_fetch_vector'),
.packages = c("dplyr", "sf", "stats")) %dopar% {

fetch_calc <- return_fetch_vector(
vector = fetch_df[i, "geom"],
origin = fetch_df$origin[i],
polygon = poly_subset
) %>%
as.data.frame() %>%
stats::setNames(c("geometry")) %>%
dplyr::mutate(
site_names = fetch_df$site_names[i],
directions = fetch_df$directions[i]
) %>%
dplyr::relocate(site_names, directions)

fetch_calc
}

# stop extra workers
doParallel::stopImplicitCluster()
parallel::stopCluster(cl)

# convert linestring dataframe to sf dataframe
final_fetch <-
fetch_loop %>%
sf::st_as_sf() %>%
sf::st_transform(sf::st_crs(polygon_layer))

message("quadrant add")
Copy link
Owner

@blasee blasee Sep 26, 2022

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Remove.


# add quadrant factor column
final_fetch$quadrant <- factor(quadrant, levels = c("North", "East", "South", "West"))

# add length of fetch linestring
final_fetch$fetch <- sf::st_length(final_fetch$geometry)

# if no SF object is desired, just return fetch distances
if(as_sf == TRUE) {
message("new part")
Copy link
Owner

@blasee blasee Sep 26, 2022

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Remove.

# make the "WindFetch" object
return(new("WindFetch", final_fetch, names = site_names, max_dist = max_dist / 1000))

if (progress_bar)
pb = txtProgressBar(max = nrow(fetch_df), style = 3)
} else {

for (i in 1:nrow(fetch_df)) {
fetch_df$fetch[i] = as.data.frame(
return_fetch_vector(fetch_df[i, "geom"],
fetch_df$origin[i],
poly_subset))
if (progress_bar)
setTxtProgressBar(pb, i)
}
cat("\n")
message("no new part")
Copy link
Owner

@blasee blasee Sep 26, 2022

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Remove.

# return only distances, no sf object
final_fetch <-
final_fetch %>%
sf::st_drop_geometry()

fetch_df$fetch = st_sfc(lapply(fetch_df$fetch, `[[`, 1),
crs = st_crs(polygon_layer))
return(final_fetch)

}

st_geometry(fetch_df) = fetch_df$fetch
fetch_df = fetch_df[, 1:2]
fetch_df$quadrant = factor(quadrant,
levels = c("North", "East", "South", "West"))
fetch_df$fetch = st_length(fetch_df$geom)

new("WindFetch", fetch_df, names = site_names, max_dist = max_dist / 1000)
}

Loading