I guess the problem is that if you represent a regular raster with polygons, and the polygons won't tell you for points falling on polygon edges or nodes that they should intersect with only one, rather than with two (edge) or four (nodes) polygons.
The reprex below shows the difference of doing an intersection of three
points with a raster represented as polygons (at the end), or as a stars
object (begin). The stars object considers only the left and upper edge
and upper-left corner as part of the grid cell; the resulting plot is
attached.
library(stars)
# Loading required package: abind
# Loading required package: sf
# Linking to GEOS 3.7.0, GDAL 2.4.0, PROJ 5.2.0
# create 10 x 10 raster
st_as_stars(matrix(1:100, 10, 10)) %>%
st_set_dimensions(names = c("x", "y")) %>%
st_set_dimensions("x", values = 0:9) %>%
st_set_dimensions("y", values = 9:0) -> r
# three points: at node, edge, and inside cell:
p3 = st_sfc(st_point(c(1,1)), st_point(c(1, 1.5)), st_point(c(1.1, 1.8)))
image(r, axes = TRUE, text_values = TRUE)
plot(p3, col = 'green', add = TRUE)
# query three points:
st_intersects(r, p3, transpose = TRUE)
# Sparse geometry binary predicate list of length 3, where the predicate
was `intersects'
# 1: 82
# 2: 72
# 3: 72
# query three points with same raster as polygons:
st_intersects(p3, st_as_sf(r))
# Sparse geometry binary predicate list of length 3, where the predicate
was `intersects'
# 1: 71, 72, 81, 82
# 2: 71, 72
# 3: 72
# raster does the same as stars:
library(raster)
# Loading required package: sp
extract(as(r, "Raster"), as(p3, "Spatial"))
# [1] 82 72 72
On 5/5/19 4:49 AM, Roozbeh Valavi wrote:
> Thank you for your explanation, Tom.
>
> Let me explain the case clearly.
> The point layer is a very dense grid field collection across Great
> Britain. I wanted to do spatial cross-validation by specifying each
> point in one of the polygons (and assign a group of the polygons to a
> fold). Since the points have been sampled on a regular grid they
> sometimes fall exactly on the border of the polygons. If I use any sort
> of intersection or spatial join, the points on the borders cause a problem.
>
> Currently, I use st_jitter() in sf package give a bit of randomness to
> the points. I wanted to see if there is another systematic way to avoid
> this.
>
> Thanks,
> Roozbeh
>
> PhD Candidate
> The Quantitative & Applied Ecology Group <http://qaeco.com/>
> School of BioSciences | Faculty of Science
> The University of Melbourne, VIC 3010, Australia
> M: +61 423 283 238 E: [email protected]
> <mailto:[email protected]>
>
>
> On Wed, May 1, 2019 at 3:16 PM Tom Philippi <[email protected]
> <mailto:[email protected]>> wrote:
>
> Roozbeh--
>
> What do you want to have happen with those points? Given the units
> in your figure, you could unambiguously assign them to the
> upper-right by adding a small value (0.000001) to both the X and Y
> coordinates of your points. Whether that is a sensible thing to do
> depends on what your data are.
>
> How or why do the points fall exactly on your polygon boundaries?
> What process puts them exactly on the polygon (grid cell)
> boundaries? It appears your polygon boundaries are a grid at
> multiples of 0.05, although not all potential cells in a rectangle
> are shown. Are there many other points not on polygon boundaries
> not shown in your figure, or are points only along grid lines? If
> the latter, does it make sense to assign them to _polygons_ rather
> than line segments?
>
> If you're doing a sample draw or measuring locations with only a few
> decimal points of accuracy, and generating your polygon boundaries
> at exact multiples of 0.05, then it might make sense to shift all of
> your points + 1/10th of your point location resolution so that their
> coordinates never fall on polygon boundaries. There _might_ be
> situations where instead of always adding (shifting up and right),
> you should randomly add or subtract in each direction for each
> point, but I'm stuck viewing this as points in sTomample units so I
> can't think of such a case.
>
> Tom
> .
>
> On Tue, Apr 30, 2019 at 7:18 PM Roozbeh Valavi
> <[email protected]
> <mailto:[email protected]>> wrote:
>
> Dear members,
>
> I want to do a spatial join/intersection in R. The problem is
> that some of my points are lying exactly at the border of the
> adjacent polygons (see the figure attached). So these points are
> either falling in both or none of the polygons! Is there any way
> to avoid this?
>
> Thanks in advance.
> Roozbeh
>
>
> image.png
>
>
>
> --
> *Roozbeh Valavi*
> PhD Candidate
> The Quantitative & Applied Ecology Group <http://qaeco.com/>
> School of BioSciences | Faculty of Science
> The University of Melbourne, VIC 3010, Australia
> Mobile: +61 423 283 238
> _______________________________________________
> R-sig-Geo mailing list
> [email protected] <mailto:[email protected]>
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>
>
> _______________________________________________
> R-sig-Geo mailing list
> [email protected]
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>
--
Edzer Pebesma
Institute for Geoinformatics
Heisenbergstrasse 2, 48151 Muenster, Germany
Phone: +49 251 8333081
Rplots.pdf
Description: Adobe PDF document
pEpkey.asc
Description: application/pgp-keys
_______________________________________________ R-sig-Geo mailing list [email protected] https://stat.ethz.ch/mailman/listinfo/r-sig-geo
