0

I have a shapefile consisting of hundreds of polygons depicting building footprints in a city. I want to fill these polygons with grid points. I am doing this by making a grid of points and then clipping the points grid to the polygons so that I can get a new clipped layer that is the grid points but with the invisible boundaries of the polygons. Here is my code:

library(sp)
library(rgdal)
library(raster)
library(rgeos)

polygons <- readOGR("C:/PATH_TO_MY_FILE/POLYGONS.shp")
polygons <- gBuffer(polygons, width=0, byid=TRUE)
grid <- makegrid(polygons)
grid <- SpatialPoints(grid, proj4string = CRS(proj4string(polygons)))
grid <- spTransform(grid, CRSobj = crs(polygons))
grid.clip <- grid.t[polygons, ]
grid.clip

This results in the following:

class       : SpatialPoints 
features    : 0 
crs         : +proj=lcc +lat_0=40.1666666666667 +lon_0=-74 +lat_1=41.0333333333333 +lat_2=40.6666666666667 +x_0=300000 +y_0=0 +datum=NAD83 +units=us-ft +no_defs

Running plot(grid.clip) yields nothing but blankness.

To identify the CRS of the polygons:

crs(polygons)

Which yields:

CRS arguments:
 +proj=lcc +lat_0=40.1666666666667 +lon_0=-74 +lat_1=41.0333333333333 +lat_2=40.6666666666667 +x_0=300000 +y_0=0 +datum=NAD83
+units=us-ft +no_defs

I am not sure why I am getting no features from the clipped grid points, but my guess is because there were no points identified within any of the polygons. According to the following picture of my data, it looks like this is because the dimensions of the grid are significantly far off from the scale of the polygons.

I run the following to see what I was working with before the clip:

library(mapview)
mapview(polygons) + mapview(grid)

And I see this: enter image description here

As you can see, the points are spread out so far beyond the scale of the polygons. How can I shrink down the spacing between the grid points so that I can actually see the points within the polygons? I recognize that there is no direct answer for what the scale should be, but I am looking for something reasonable. Since these are building footprints, and since this polygon layer has units in feet, I was thinking that maybe the grid cell size should be within feet. How do I set this? I think this involves the cell.size argument within makegrid(), but I am not sure how to use this. Also, what might be a good way to decide between 1-foot spacing vs 5-foot spacing vs 10-foot spacing for the grid points? I want the best resolution/highest reasonable point density without overloading RStudio and crashing.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.