[R-sig-Geo] Spatial join/intersection in R
Edzer Pebesma
edzer@pebe@m@ @end|ng |rom un|-muen@ter@de
Mon May 6 12:23:48 CEST 2019
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: rvalavi using student.unimelb.edu.au
> <mailto:rvalavi using student.unimelb.edu.au>
>
>
> On Wed, May 1, 2019 at 3:16 PM Tom Philippi <tephilippi using gmail.com
> <mailto:tephilippi using gmail.com>> 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
> <rvalavi using student.unimelb.edu.au
> <mailto:rvalavi using student.unimelb.edu.au>> 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
> R-sig-Geo using r-project.org <mailto:R-sig-Geo using r-project.org>
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>
>
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo using r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>
--
Edzer Pebesma
Institute for Geoinformatics
Heisenbergstrasse 2, 48151 Muenster, Germany
Phone: +49 251 8333081
-------------- next part --------------
A non-text attachment was scrubbed...
Name: Rplots.pdf
Type: application/pdf
Size: 5374 bytes
Desc: not available
URL: <https://stat.ethz.ch/pipermail/r-sig-geo/attachments/20190506/29fbcd05/attachment.pdf>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: pEpkey.asc
Type: application/pgp-keys
Size: 2472 bytes
Desc: not available
URL: <https://stat.ethz.ch/pipermail/r-sig-geo/attachments/20190506/29fbcd05/attachment.bin>
More information about the R-sig-Geo
mailing list