3

I'm trying to get the coordinates of a set of points defining a grid within a polygon (which I have a shapefile for). It seemed like the simplest thing to do would be to create a grid of points, and then filter those points down to only the ones within the polygon. I looked at https://gis.stackexchange.com/questions/133625/checking-if-points-fall-within-polygon-shapefile and Convert a shapefile from polygons to points?, and based on the answers there I tried this:

library(rgdal)
city_bdry <- readOGR("Boundaries - City",
                     "geo_export_32ded882-2eab-4eaa-b9da-a18889600a40")
res <- 0.01
bb <- bbox(city_bdry)
gt <- GridTopology(cellcentre.offset = bb[,1], cellsize = c(res, res), 
                   cells.dim = c(diff(bb[,1]), diff(bb[2,])) / res + 1)
pts <- SpatialPoints(gt, proj4string = CRS(proj4string(city_bdry)))
ov <- over(pts, city_bdry)

The result, however, doesn't include the actual coordinates of the points that overlap the polygon, so it's useless to me. How can I get that information to be included in the dataframe? Or, is there a simpler way to do what I'm trying to do?

The shapefile I'm using can be downloaded from https://data.cityofchicago.org/Facilities-Geographic-Boundaries/Boundaries-City/ewy2-6yfk

Hack-R
  • 22,422
  • 14
  • 75
  • 131
Empiromancer
  • 3,778
  • 1
  • 22
  • 53
  • 1
    Use `splancs::inout()`. See the last answer on this thread: https://stackoverflow.com/questions/43436466/create-grid-in-r-for-kriging-in-gstat/45948442#45948442 I've dealt with this exact problem before, and `inout()` is the most simple solution I've found. – Rich Pauloo Sep 05 '17 at 19:23
  • 1
    @RichPauloo This solved my problem, thank you. If you make it an answer, I'll mark it as accepted. – Empiromancer Sep 05 '17 at 22:11
  • I'm glad it worked! This is such a common task with geospatial data, and it's one example where an answer wasn't as easy to find online as I would have thought. Extracting the outline and then using inout() is a two step process, and I wonder if anyone out there (maybe on GIS SE) has a more simple one-line solution within one of the spatial packages like sp. Anyone? – Rich Pauloo Sep 05 '17 at 22:28

3 Answers3

3

If I got you right, you could try

library(rgdal)
download.file("https://data.cityofchicago.org/api/geospatial/ewy2-6yfk?method=export&format=Shapefile", tf<-tempfile(fileext = ".zip"), mode="wb")
unzip(tf, exdir=dir<-file.path(tempdir(), "Boundaries - City"))
city_bdry <- readOGR(dir, tools::file_path_sans_ext((list.files(dir)[1])))
res <- 0.01
bb <- bbox(city_bdry)
gt <- GridTopology(cellcentre.offset = bb[,1], cellsize = c(res, res), 
                   cells.dim = c(diff(bb[,1]), diff(bb[2,])) / res + 1)
pts <- SpatialPoints(gt, proj4string = CRS(proj4string(city_bdry)))

ov <- sp::over(pts, as(city_bdry, "SpatialPolygons"))
pts_over <- pts[!is.na(ov)]

plot(city_bdry)
points(pts_over)
coordinates(pts_over)
lukeA
  • 53,097
  • 5
  • 97
  • 100
  • @Empiromancer Hm which point is reordered? I guess if the order was changed, `plot(city_bdry);points(pts[!is.na(ov)])` and `plot(city_bdry);points(pts[is.na(ov)])` would not work. – lukeA Sep 05 '17 at 19:45
  • I don't know for sure it does reorder anything, it just seems to rely on order being preserved and I'm wondering if that behavior is guaranteed by the functions involved. – Empiromancer Sep 05 '17 at 22:13
  • @Empiromancer There's no reordering. It just "Find[s] coordinates of grid points within polygon", as the title asked. – lukeA Sep 05 '17 at 22:25
2

Use splancs::inout().

1. Get the outline of your polygon

outline <- mySpatialPolygonsDataFrame@polygons[[2]]@Polygons[[1]]@coords

2. Use inout() to find what points are within the outline 

library(splancs)
pts_in_polygon <- points[inout(points,outline), ]

Note: very similar to the answer I provide to create an irregularly shaped grid (especially for kriging.)

Rich Pauloo
  • 7,734
  • 4
  • 37
  • 69
1

You can also achieve this by simply subsetting the data with [ like yourPoints[yourPolygons, ]: 

library(raster)
bra <- getData(country = "BRA", level = 1)

pts <- makegrid(bra, 100)
pts <- SpatialPoints(pts, proj4string = CRS(proj4string(bra)))

coordinates(pts[bra, ])
loki
  • 9,816
  • 7
  • 56
  • 82