How do I generate a grid of coordinates within specified polygons in R?

Question

I have a list of two-dimensional polygons defined as two-column matrices of x and y coordinates in R. They completely fill a square area and are mutually exclusive. I want to use these polygon definitions to generate a fine grid of x,y coordinate values in which each value is identified by which polygon it falls into.

I have explored the sp package and can get my polygons into an object of class SpatialPolygons, but I don't know if that gets me closer to my goal. With my polygons in a dataframe, I can use ggplot with geom_polygon(aes(fill=ID)) to generate a plot of the polygons with coloring based on polygon ID.

I can see several paths forward, but don't know how to accomplish any of them:

1. A function that takes a polygon and generates a uniform grid of coordinates within the polygon boundaries. (My polygons are quite irregular, with many sides, so creating a custom function for them would be painful and error prone.)

2. A function that takes a pair of x, y coordinates and my list of polygons and outputs which polygon the coordinates fall into.

3. A function that takes my ggplot-generated plot and converts the colors into a grid of numeric coordinate values that I could read back into R.

There may be also be other approaches that I'm not imagining. I have to believe that other people have had this same need before, but extensive searching has not led me to any existing functions that do what I need.

Answer 1

[mumble about spsample deleted]

Hmm in the cold light of day it seems you want something else:

IF all your polygons make a rectangle AND you want a regular grid of points over that rectangle THEN create a SpatialPoints object of grid coordinates (see 'expand.grid' for part of the solution to that sub-problem) AND THEN use 'overlay' from package:sp to test what polygon your grid points are in.

You might also want to use bbox to get the extent of your polygons.

Answer 2

It sounded like you were doing this positioning on a square grid so it may be simpler than the more general polygon approach would require. Let's say your coordinates for this grid-on-the-square are two vectors, 'xx' and 'yy', and you have list of points in a data.frame or matrix named 'mypoints'. This will create a matrix of row-col-indices to look up the proper sub-square:

 xx <- seq(0,1,by=.1)
 yy <- seq(0,1,by=.1)
 mypoints <- matrix(runif(10), ncol=2)
 head(mypoints)
#---------------
          [,1]      [,2]
[1,] 0.7731868 0.2707768
[2,] 0.7005779 0.7881789
[3,] 0.9520941 0.6661852
[4,] 0.4625906 0.9176813
[5,] 0.4550811 0.5017386
#---------------
 findInterval(mypoints[1:5,1], xx)
#[1]  8  8 10  5  5
 findInterval(mypoints[1:5,2], yy)
#[1]  3  8  7 10  6
 pointidxs <- matrix( c( findInterval(mypoints[,1], xx), 
                         findInterval(mypoints[,2], yy) ), ncol=2)
 head(pointidxs)
#--------------
     [,1] [,2]
[1,]    8    3
[2,]    8    8
[3,]   10    7
[4,]    5   10
[5,]    5    6

Answer 3

I haven't thought very much, but here is a before-coffee idea: I understand that your polygons form a Voronoi tesselation. Now, it is supposed to be easy to obtain the corresponding Delaunay triangulation, which should give you a straight-forward way to decide whether a particular point belongs to the corresponding polygon.

Hope that makes sense?