-
Notifications
You must be signed in to change notification settings - Fork 0
/
fetch.R
307 lines (271 loc) · 13.1 KB
/
fetch.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
#' Calculate Wind Fetch
#'
#' Wind fetch is the unobstructed length of water over which wind can blow, and
#' it is commonly used as a measure of exposure to wind and waves at coastal
#' sites. The \code{fetch} function automatically calculates the wind fetch for
#' marine locations within the boundaries of the specified coastline layer.
#' This allows wind fetch to be calculated anywhere around the globe.
#'
#' The function takes a \code{\link[sp]{SpatialPolygons-class}} object
#' (\code{polygon_layer}) that represents the coastline, surrounding islands,
#' and any other obstructions, and calculates the wind fetch for every specifed
#' direction. This is calculated for all the user-defined sites, that are
#' represented as the point geometries in a
#' \code{\link[sp]{SpatialPoints-class}} object.
#'
#' The directions for which the wind fetch are calculated for each site are
#' determined by the number of directions per quadrant (\code{n_directions}).
#' The default value of 9 calculates 9 fetch vectors per quadrant (90 degrees),
#' or equivalently, one fetch vector every 10 degrees. The first fetch vector is
#' always calculated for the northerly direction (0/360 degrees).
#'
#' @param polygon_layer \code{\link[sp]{SpatialPolygons}}* object where the
#' polygon geometries represent any obstructions to fetch
#' calculations including the coastline, islands and/or
#' exposed reefs.
#' @param site_layer \code{\link[sp]{SpatialPoints}}* object where the point
#' geometries represent the site locations.
#' @param max_dist numeric. Maximum distance in kilometres (default 300). This
#' will need to be scaled manually if the units for the CRS are
#' not 'm'.
#' @param n_directions numeric. The number of fetch vectors to calculate per
#' quadrant (default 9).
#' @param site_names character vector of the site names. If missing, the site
#' names are taken from a column of the data associated with
#' \code{site_layer} matching the regular expression
#' \code{^[Nn]ames{0,1}}. If there is no such column, then
#' default names are created ('Site 1', 'Site 2', ...).
#' @param quiet logical. Suppress diagnostic messages? (Default \code{FALSE}).
#'
#' @return Returns a \code{\link{Fetch}} object.
#'
#' @note At least one of the inputs to the \code{polygon_layer} or
#' \code{site_layer} arguments must be projected. If one of the inputs are
#' not projected, then it will be transformed to have the same projection
#' as the other. If both are projected, but do not have identical
#' coordinate reference systems (CRS) then \code{site_layer} will be
#' transformed to the same CRS as \code{polygon_layer}.
#'
#' @seealso \code{\link[rgdal]{spTransform}} for methods on transforming map
#' projections and datum.
#' @seealso \code{\link[sp]{is.projected}} for checking whether a spatial object
#' is projected.
#' @seealso \code{\link{fetchR}} for an overview of this package with an
#' extensive, reproducible example.
#' @seealso \code{\link{summary,Fetch-method}} for summarising the fetch lengths.
#'
#' @importFrom rgdal CRSargs
#' @importFrom sp SpatialPoints CRS over spTransform coordinates SpatialLinesLengths SpatialLinesDataFrame identicalCRS
#' @importFrom rgeos gBuffer gIntersects gIntersection
#' @importFrom utils head
#' @importFrom methods new is
#' @import sp
#' @examples
#'
#' # Create the polygon layer ----------------------------------------
#' #
#' # This is the layer that represents any obstacles that obstruct wind flow.
#'
#' # Import map data for the Philippines.
#' philippines.df = ggplot2::map_data("world", region = "Philippines")
#'
#' # Create a list for each separate polygon
#' philippines.list = split(philippines.df[, c("long", "lat")],
#' philippines.df$group)
#'
#' library(sp)
#'
#' philippines.Poly = lapply(philippines.list, Polygon)
#' philippines.Polys = list(Polygons(philippines.Poly, ID = "Philippines"))
#'
#' # Include CRS information to make it a SpatialPolygons object
#' philippines.sp = SpatialPolygons(philippines.Polys,
#' proj4string = CRS("+init=epsg:4326"))
#'
#' # Create the points layer ----------------------------------------
#' #
#' # The points layer represents the locations for which the wind fetch needs to
#' # be calculated.
#'
#' # We need to calculate wind fetch for the following 3 sites:
#' sites.df = data.frame(lon = c(124.4824, 125.8473, 124.8416),
#' lat = c(9.167999, 9.751394, 11.478243),
#' site = c("Camiguin Island", "Bucas Grande Island",
#' "Talalora"))
#'
#' # Create the SpatialPoints object
#' sites.sp = SpatialPoints(sites.df[, 1:2], CRS("+init=epsg:4326"))
#'
#' # Map projection -------------------------------------------------
#' #
#' # At least one of the polygon or points layers need to be projected to
#' # calculate wind fetch.
#'
#' # All these locations lie within the Philippines zone 5 / PRS92, that has
#' # WGS84 Bounds: 123.8000, 5.3000, 126.7000, 12.7500
#' # (http://spatialreference.org/ref/epsg/3125/)
#' # This suggests that this is a suitable map projection.
#' philippines.proj = spTransform(philippines.sp, "+init=epsg:3125")
#'
#' # Calculate wind fetch -------------------------------------------
#' #
#' # Calculate wind fetch at all the 3 locations for every 10 degrees on the
#' # compass rose, with a maximum distance for any fetch vector of 300 km.
#' my_fetch = fetch(philippines.proj, sites.sp, site_names = sites.df$site)
#' my_fetch
#'
#' # Return only the summary data frame
#' summary(my_fetch)
#'
#' # Transform the fetch vectors back to the original CRS
#' my_fetch_latlon = spTransform(my_fetch, proj4string(philippines.sp))
#'
#' # Return the raw data in the original, lat/lon coordinates
#' my_fetch_latlon.df = as(my_fetch_latlon, "data.frame")
#' my_fetch_latlon.df
#'
#' # Plot the wind fetch vectors ------------------------------------
#'
#' # Plot the fetch vectors in the projected space...
#' plot(my_fetch, philippines.proj, axes = TRUE)
#'
#' # ... or in the original coordinate reference system
#' plot(my_fetch, philippines.sp, axes = TRUE)
#'
#' # Output to KML --------------------------------------------------
#' \dontrun{
#'
#' # Save a KML file in the current working directory.
#' kml(my_fetch, colour = "white")
#' }
#' @export
fetch = function(polygon_layer, site_layer, max_dist = 300, n_directions = 9,
site_names, quiet = FALSE){
if (!is(polygon_layer, "SpatialPolygons"))
stop(paste("polygon_layer must be a SpatialPolygons object.\nSee",
"'?SpatialPolygons' for details on how to create a",
"SpatialPolygons object."), call. = FALSE)
if (!is(site_layer, "SpatialPoints"))
stop(paste("site_layer must be a SpatialPoints object.\nSee",
"'?SpatialPoints' for details on how to create a SpatialPoints",
"object."), call. = FALSE)
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)
stop("n_directions must be a single integer.", call. = FALSE)
n_directions = round(n_directions)
if (n_directions < 1 || n_directions > 20)
stop("n_directions must be between 1 and 20.", call. = FALSE)
if (!missing(site_names)){
site_names = as.character(site_names)
if (length(site_names) != length(site_layer)){
warning(paste("lengths differ for the number of sites and site names;",
"using default names instead."), call. = FALSE)
site_names = paste("Site", seq_along(site_layer))
}
} else {
if (any(grepl("^[Nn]ames{0,1}$", names(site_layer)))){
name_col = grep("^[Nn]ames{0,1}$", names(site_layer))
site_names = as.character(data.frame(site_layer)[, name_col[[1]]])
} else {
site_names = paste("Site", seq_along(site_layer))
}
}
quiet = as.logical(quiet[1])
## Check if the polygon and points layers are projected, and ensure they have
## the same CRS.
which_proj = c(is.projected(polygon_layer), is.projected(site_layer))
if (all(!which_proj))
stop("polygon_layer and/or site_layer must be projected to calculate fetch",
call. = FALSE)
if (all(which_proj) && !identicalCRS(polygon_layer, site_layer)){
warning("the CRS for polygon_layer and site_layer differ; transforming
site_layer CRS to match")
site_layer = spTransform(site_layer, CRS(proj4string(polygon_layer)))
}
if (!which_proj[1]){
if (!quiet)
message("projecting polygon_layer onto the site_layer CRS")
polygon_layer = spTransform(polygon_layer, CRS(proj4string(site_layer)))
}
if (!which_proj[2]){
if (!quiet)
message("projecting site_layer onto the polygon_layer CRS")
site_layer = spTransform(site_layer, CRS(proj4string(polygon_layer)))
}
if (!quiet)
message("checking site locations are not on land")
# Should remove these sites with warning instead of returning error.
if (any(!is.na(over(polygon_layer, site_layer))))
stop("at least one site location is on land")
# Convert max_dist to appropriate units.
# First of all convert max_dist to metres (default)
max_dist = max_dist * 1000
# Double check if metres are the correct units
proj_unit = strsplit(gsub("*.+units=", "",
CRSargs(CRS(proj4string(polygon_layer)))), " ")[[1]][1]
# If not, warn the user that the supplied max_dist should be scaled
# appropriately
if (proj_unit != "m")
warning("the PROJ.4 unit is not metres; ensure max_dist has been scaled
appropriately")
directions = head(seq(0, 360, by = 360 / (n_directions * 4)), -1)
# Return the quadrant the directions belong to
dirs = as.numeric(directions)
dirs_bin = findInterval(dirs, seq(45, 315, by = 90))
quadrant = rep("North", length(dirs))
quadrant[dirs_bin == 1] = "East"
quadrant[dirs_bin == 2] = "South"
quadrant[dirs_bin == 3] = "West"
# Rearrange sequence order to start at 90 degrees to match up with the output
# from gBuffer
directions = unlist(split(directions, directions < 90), use.names = FALSE)
# Create an empty list to store the Fetch objects
fetch_list = vector("list", length(site_layer))
for (i in seq_along(site_layer)){
message("calculating fetch for ", site_names[i], " (", i, " out of ",
length(site_layer), ")")
## Create polygon (approximating a circle) with a given radius. These
## vertices are used for creating the end points for the fetch vectors.
d_bff = gBuffer(site_layer[i, ], width = max_dist, quadsegs = n_directions)
## Calculate end points at the maximum distances.
fetch_ends = head(coordinates(d_bff@polygons[[1]]@Polygons[[1]]), -1)
fetch_ends = fetch_ends[order(directions), ]
## Create a list of Line objects radiating from the site location to the
## maximum distance for each direction.
line_list = create_line_list(site_layer[i, ], fetch_ends)
## Create a SpatialLines object (i.e. add in CRS information)
fetch_sp_lines = create_sp_lines(line_list, sort(directions), polygon_layer)
## Subset polygon coastline layer to only incorporate polygons that
## interfere with any of the fetch vectors. This will speed up computation
## times.
poly_layer_subset =
polygon_layer[which(!is.na(over(polygon_layer, fetch_sp_lines))), ]
if (length(poly_layer_subset) > 0){
## Calculate the fetch bearings that hit land
hit_land = !sapply(gIntersects(fetch_sp_lines, poly_layer_subset,
byid = c(TRUE, FALSE), returnDense = FALSE),
is.null)
## Calculate intersections and identify closest shoreline for those
## vectors that hit land
ints = gIntersection(fetch_sp_lines[hit_land], poly_layer_subset,
byid = c(TRUE, FALSE))
fetch_ends[hit_land, ] = t(sapply(ints@lines, function(x){
coordinates(x)[[1]][1, ]
}))
## Update the line list and spatialLines objects
line_list = create_line_list(site_layer[i, ], fetch_ends)
fetch_sp_lines = create_sp_lines(line_list, sort(directions), polygon_layer)
}
fetch.df = data.frame(site = site_names[i],
fetch = SpatialLinesLengths(fetch_sp_lines) / 1000,
direction = sort(directions),
quadrant = factor(quadrant,
levels = c("North", "East",
"South", "West")))
## Create a SpatialLinesDataFrame object to include the fetch lengths, and
## add it to the Fetch list
fetch_list[[i]] = SpatialLinesDataFrame(fetch_sp_lines, fetch.df)
}
new("Fetch", fetch_list, names = site_names, max_dist = max_dist / 1000)
}