Skip to content

Commit

Permalink
let st_make_grid() cover geometry, not bounding box
Browse files Browse the repository at this point in the history
addresses #1260
  • Loading branch information
edzer committed Feb 9, 2020
1 parent 6351350 commit 89a4788
Show file tree
Hide file tree
Showing 7 changed files with 27 additions and 15 deletions.
2 changes: 1 addition & 1 deletion R/join.R
Expand Up @@ -104,7 +104,7 @@ st_join = function(x, y, join, ...) UseMethod("st_join")
#' nc <- st_transform(st_read(system.file("shape/nc.shp", package="sf")), 2264)
#' gr = st_sf(
#' label = apply(expand.grid(1:10, LETTERS[10:1])[,2:1], 1, paste0, collapse = " "),
#' geom = st_make_grid(nc))
#' geom = st_make_grid(st_as_sfc(st_bbox(nc))))
#' gr$col = sf.colors(10, categorical = TRUE, alpha = .3)
#' # cut, to check, NA's work out:
#' gr = gr[-(1:30),]
Expand Down
17 changes: 12 additions & 5 deletions R/make_grid.R
@@ -1,6 +1,6 @@
#' Create a regular tesselation over the bounding box of an sf or sfc object
#'
#' Create a square or hexagonal grid over the bounding box of an sf or sfc object
#' Create a square or hexagonal grid covering the geometry of an sf or sfc object
#' @param x object of class \link{sf} or \link{sfc}
#' @param cellsize target cellsize
#' @param offset numeric of lengt 2; lower left corner coordinates (x, y) of the grid
Expand All @@ -11,6 +11,8 @@
#' @param flat_topped logical; if \code{TRUE} generate flat topped hexagons, else generate pointy topped
#' @return Object of class \code{sfc} (simple feature geometry list column) with, depending on \code{what} and \code{square},
#' square or hexagonal polygons, corner points of these polygons, or center points of these polygons.
#' @details to obtain a grid covering the bounding box of a set of geometries,
#' pass \code{st_as_sfc(st_bbox(x))} for argument \code{x}
#' @examples
#' plot(st_make_grid(what = "centers"), axes = TRUE)
#' plot(st_make_grid(what = "corners"), add = TRUE, col = 'green', pch=3)
Expand Down Expand Up @@ -93,10 +95,15 @@ st_make_grid = function(x,
} else
stop("unknown value of `what'")

ret = st_sfc(ret, crs = crs)

if (missing(x))
st_sfc(ret, crs = crs)
ret
else if (what != "polygons" || min(st_dimension(x)) < 2)
ret[x]
else
st_sfc(ret, crs = st_crs(x))
ret[lengths(st_relate(ret, x, "2********")) > 0]
# overlap dim at least equal to the mininmum of that of x geoms
}


Expand Down Expand Up @@ -159,8 +166,8 @@ make_hex_grid = function(obj, pt, dx, what, flat_topped = TRUE) {
else # points:
st_sfc(lapply(seq_len(nrow(centers)), function(i) mk_pol(centers[i,])), crs = st_crs(bb))

if (what == "points")
if (what == "points" || min(st_dimension(obj)) < 2)
ret[obj]
else
else
ret[lengths(st_relate(ret, obj, "2********")) > 0] # part or total overlap
}
2 changes: 1 addition & 1 deletion R/sfc.R
Expand Up @@ -135,7 +135,7 @@ sfg_is_empty = function(x) {
"[.sfc" = function(x, i, j, ..., op = st_intersects) {
if (!missing(i) && (inherits(i, "sf") || inherits(i, "sfc") || inherits(i, "sfg")))
i = lengths(op(x, i, ...)) != 0
st_sfc(NextMethod(), crs = st_crs(x), precision = st_precision(x))
st_sfc(unclass(x)[i], crs = st_crs(x), precision = st_precision(x))
}


Expand Down
2 changes: 1 addition & 1 deletion man/st_join.Rd

Some generated files are not rendered by default. Learn more about how customized files appear on GitHub.

6 changes: 5 additions & 1 deletion man/st_make_grid.Rd

Some generated files are not rendered by default. Learn more about how customized files appear on GitHub.

11 changes: 6 additions & 5 deletions tests/geos.R
Expand Up @@ -36,7 +36,8 @@ st_combine(nc)
st_dimension(st_sfc(st_point(0:1), st_linestring(rbind(c(0,0),c(1,1))),
st_polygon(list(rbind(c(0,0), c(1,0), c(1,1), c(0,1), c(0,0))))))

g = st_make_grid(nc)
ncbb = st_as_sfc(st_bbox(nc))
g = st_make_grid(ncbb)
x = st_intersection(nc, g)
x = st_intersection(g, nc)

Expand All @@ -48,7 +49,7 @@ set.seed(13531)

st_line_sample(ls, density = 1, type = "random")

g = st_make_grid(nc, n = c(20,10))
g = st_make_grid(ncbb, n = c(20,10))

a1 = st_interpolate_aw(nc["BIR74"], g, FALSE)
sum(a1$BIR74) / sum(nc$BIR74) # not close to one: property is assumed spatially intensive
Expand All @@ -62,8 +63,8 @@ length(g)
g = st_make_grid(what = "corners")
length(g)

g1 = st_make_grid(nc, 0.1, what = "polygons", square = FALSE)
g2 = st_make_grid(nc, 0.1, what = "points", square = FALSE)
g1 = st_make_grid(ncbb, 0.1, what = "polygons", square = FALSE)
g2 = st_make_grid(ncbb, 0.1, what = "points", square = FALSE)

# st_line_merge:
mls = st_multilinestring(list(rbind(c(0,0), c(1,1)), rbind(c(2,0), c(1,1))))
Expand Down Expand Up @@ -93,7 +94,7 @@ i = st_intersects(ncm, ncm[1:88,])
all.equal(i, t(t(i)))

# check use of pattern in st_relate:
sfc = st_sfc(st_point(c(0,0)), st_point(c(3,3)))
sfc = st_as_sfc(st_bbox(st_sfc(st_point(c(0,0)), st_point(c(3,3)))))
grd = st_make_grid(sfc, n = c(3,3))
st_intersects(grd)
st_relate(grd, pattern = "****1****")
Expand Down
2 changes: 1 addition & 1 deletion tests/sfc.R
Expand Up @@ -174,7 +174,7 @@ st_join(a, b, left = FALSE)
nc <- st_transform(st_read(system.file("shape/nc.shp", package="sf")), 2264)
gr = st_sf(
label = apply(expand.grid(1:10, LETTERS[10:1])[,2:1], 1, paste0, collapse = " "),
geom = st_make_grid(nc))
geom = st_make_grid(st_as_sfc(st_bbox(nc))))
gr$col = sf.colors(10, categorical = TRUE, alpha = .3)
# cut, to check, NA's work out:
gr = gr[-(1:30),]
Expand Down

1 comment on commit 89a4788

@michaelvalcic
Copy link

Choose a reason for hiding this comment

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

Thank you for making this change.

Please sign in to comment.