I am trying to create coordinate points based on a polygon that I've created in Google Maps and imported into R as an sf object. Specifically I have sample points of physical samples I've taken, but want the coordinates of every vine in the vineyard I sampled. The vines are spaced regularly in an 8ft * 8ft square over a rectangle with 26 rows and 59 columns. How can I create a dataframe of points?

Here's what I've done thus far to create regularly spaced points within the polygon extent.

i <- grep("California", epsg$note, ignore.case=TRUE)
   # first three
epsg[i[1:15], ]
# 717                +proj=lcc +lat_1=41.66666666666666 +lat_2=40       +lat_0=39.33333333333334 +lon_0=-122 +x_0=2000000.0001016 +y_0=500000.0001016001 +datum=NAD83 +units=us-ft +no_defs
# 718 +proj=lcc +lat_1=39.83333333333334 +lat_2=38.33333333333334 +lat_0=37.66666666666666 +lon_0=-122 +x_0=2000000.0001016 +y_0=500000.0001016001 +datum=NAD83 +units=us-ft +no_defs
# 719            +proj=lcc +lat_1=38.43333333333333 +lat_2=37.06666666666667 +lat_0=36.5 +lon_0=-120.5 +x_0=2000000.0001016 +y_0=500000.0001016001 +datum=NAD83 +units=us-ft +no_defs
f <- readOGR("dem_shape.kml")
#-122.8037, -122.8019, 38.6641, 38.66549  (xmin, xmax, ymin, ymax)
newcrs <- CRS("+proj=lcc +lat_1=39.83333333333334 +lat_2=38.33333333333334 +lat_0=37.66666666666666 +lon_0=-122 +x_0=2000000 +y_0=500000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs")
rob <- spTransform(f, newcrs)
#=9900.83 m^2
plot(rob) ; points(spsample(rob, n=1664.649, type='regular'), col='red', pch=3, cex=0.5)```

The issue is that, if each cross in the plotted image were a vine, they are not spaced in a square pattern in realtion to the polygon extent. 
  • 1
    Some example data and information on what you have tried would be appreciated. Please read our FAQ for information on how to ask a good question. There is a basic expectation that you have already made a good faith effort to tackle your problem and that we are not just acting as a coding service. May 1 '19 at 16:28
  • Thank you very much for the feedback, I hadn't read the FAQ.
    – getthere
    May 4 '19 at 17:08

You really need the coordinate of at least one sample point to extrapolate the sample grid. You can try sp::spsample with type = "regular" and offset = 8 or sp::makegrid but that will not necessarily align your results with the field samples.

You can also create an empty raster with row=26, col=59 and then turn it into a SpatialPolygonsDataFrame where each polygon represents a sample location. Please note that if your polygon is irregular you will likely have to rasterize it to correctly represent the sample frame. You could also try the method below and use something like raster::intersect to clip the grid to an irregular polygon.

Here is a quick example using a square polygon. First, add libraries and create a polygon to act as your sample frame

r <- raster(nrows=10, ncols=10, xmn=571823.6, xmx=616763.6, ymn=4423540, 
            ymx=4453690, crs = CRS("+proj=utm +zone=12 +datum=NAD83 
            +units=m +no_defs +ellps=GRS80 +towgs84=0,0,0"))
 r[] <- round(runif(ncell(r)),0) 
    p <- rasterToPolygons(r, fun=function(x){x == 1})[1,]

Then we can create a raster, using your sample polygon, and coerce it to a polygon grid. We can then take the centroid to create points. Getting the fixed sample distance would be dependent on your polygon dimensions with the nrow=26 and ncol=59 working out to the desired sampling distance.

r <- raster(extent(p), nrows=26, ncols=59, crs = CRS("+proj=utm +zone=12 +datum=NAD83 
            +units=m +no_defs +ellps=GRS80 +towgs84=0,0,0"))
  r[] <- rep(1,ncell(r))
g <- as(r, "SpatialPolygonsDataFrame")
  g.pts <- SpatialPoints(coordinates(rp))

  points(g.pts, pch=19, cex=0.5)

Same thing but, using sf (coercing sp polygon object to sf on the fly). Adding st_intersection makes it so that the resulting sample grid is clipped to the sample polygon thus, dealing with an irregular polygon. However, this could potentially truncate the number of desired rows and columns.

p.grid <- st_make_grid(st_as_sf(p), square = TRUE, n = c(59,26) ) %>% 
          st_intersection(st_as_sf(p)) %>% 
  • Thank you very much for this solution. I've added some code that I tried previously above for clarification.
    – getthere
    May 4 '19 at 17:09

Your Answer

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

Not the answer you're looking for? Browse other questions tagged or ask your own question.