How to create 500m * 500m grids inside a sf polygon by using sf package of r?

Question

I'm trying to create grids inside the boundary of Suffolk County, NY whose class is "sf". I named the layer "SUFF". Through using st_area(SUFF), I got to know that the area of the county is 6136105813 square meter.

So, I decided to create the rectangular grid with the size of 500 meter * 500 meter. I wrote the code: fishnet <- st_make_grid(st_transform(SUFF, crs=st_crs(4326)),cellsize = 500, square = TRUE) %>% st_sf().

However, I only got one grid. Fishnet for cellsize = 500 And then I tried many different cell size values and I found that I would got 1 grid if cellsize >= 1, 4 grids if cellsize = 0.5, 32 grids if cellsize = 0.25... Fishnet for cellsize = 0.25

In my understanding, the unit of the cell size should be the same as the SUFF layer, which is meter, is that correct? Would you mind giving me some guidance that how I can create 500m * 500m grids by using st_make_grid()?

Answer 1

I have a Suffolk county in data below. AS @d-j points out, the bounding box is the smallest square that contains all of Suffolk, in this case. Depending on the shape of the county, relative to a square, many of your grids likely will not be over Suffolk, and this should be considered if you're interpolating values onto Suffolk that make sense to be, say, on land. Using a larger grid cellsize:

suff_grid <- st_transform(suffolk_cty, crs='EPSG:5070') |>
+ st_make_grid(cellsize = 2000, square = TRUE)
plot(suff_grid)
suffolk_5070 <- st_transform(suffolk_cty, crs='EPSG:5070') #data
plot(suffolk_5070, add = TRUE)

And the grids you're probably interested in, ones that intersect with Suffolk:

plot(suff_grid[suffolk_5070, col='#ff000088', add = TRUE)

This difference is useful to consider when you're thinking to do a st_interpolate_aw loosely like this.

data:

structure(list(statefp = "36", countyfp = "103", countyns = "00974149", 
    affgeoid = "0500000US36103", geoid = "36103", name = "Suffolk", 
    namelsad = "Suffolk County", stusps = "NY", state_name = "New York", 
    lsad = "06", aland = 2359930478, awater = 3786764809, state_name.1 = "New York", 
    state_abbr = "NY", jurisdiction_type = "state", geometry = structure(list(
        structure(list(list(structure(c(-72.018926, -71.926802, 
        -71.917281, -72.034754, -72.018926, 41.274114, 41.290122, 
        41.251333, 41.234818, 41.274114), .Dim = c(5L, 2L))), 
            list(structure(c(-73.485365, -73.436664, -73.392862, 
            -73.33136, -73.2358274066785, -73.229285, -73.148994, 
            -73.144673, -73.110368, -73.040445, -72.859831, -72.708069, 
            -72.585327, -72.504305, -72.445242, -72.389809, -72.354123, 
            -72.291109, -72.189163, -72.182033, -72.254704, -72.283093, 
            -72.217476, -72.162898, -72.126704, -72.084207, -72.095711, 
            -72.051928, -71.959595, -71.919385, -71.856214, -71.936977, 
            -72.097369, -72.298727, -72.39585, -72.757176, -72.923214, 
            -73.012545, -73.1460808692108, -73.20844, -73.306396, 
            -73.351465, -73.4241151206597, -73.423275, -73.435582, 
            -73.438617, -73.462259, -73.4973510386263, -73.485365, 
            40.946397, 40.934897, 40.955297, 40.929597, 40.9066897675323, 
            40.905121, 40.928898, 40.955842, 40.971938, 40.964498, 
            40.966088, 40.977851, 40.997587, 41.043329, 41.086116, 
            41.108304, 41.139952, 41.155874, 41.193549, 41.178345, 
            41.110852, 41.067874, 41.040611, 41.053187, 41.115139, 
            41.101524, 41.05402, 41.020506, 41.071237, 41.080517, 
            41.070598, 41.006137, 40.95888, 40.903151, 40.86666, 
            40.764371, 40.713282, 40.679651, 40.6464079681013, 
            40.630884, 40.620756, 40.6305, 40.6132119188686, 
            40.670483, 40.734635, 40.751287, 40.86671, 40.9231822732944, 
            40.946397), .Dim = c(49L, 2L)))), class = c("XY", 
        "MULTIPOLYGON", "sfg"))), class = c("sfc_MULTIPOLYGON", 
    "sfc"), precision = 0, bbox = structure(c(xmin = -73.4973510386263, 
    ymin = 40.6132119188686, xmax = -71.856214, ymax = 41.290122
    ), class = "bbox"), crs = structure(list(input = "EPSG:4326", 
        wkt = "GEOGCRS[\"WGS 84\",\n    ENSEMBLE[\"World Geodetic System 1984 ensemble\",\n        MEMBER[\"World Geodetic System 1984 (Transit)\"],\n        MEMBER[\"World Geodetic System 1984 (G730)\"],\n        MEMBER[\"World Geodetic System 1984 (G873)\"],\n        MEMBER[\"World Geodetic System 1984 (G1150)\"],\n        MEMBER[\"World Geodetic System 1984 (G1674)\"],\n        MEMBER[\"World Geodetic System 1984 (G1762)\"],\n        ELLIPSOID[\"WGS 84\",6378137,298.257223563,\n            LENGTHUNIT[\"metre\",1]],\n        ENSEMBLEACCURACY[2.0]],\n    PRIMEM[\"Greenwich\",0,\n        ANGLEUNIT[\"degree\",0.0174532925199433]],\n    CS[ellipsoidal,2],\n        AXIS[\"geodetic latitude (Lat)\",north,\n            ORDER[1],\n            ANGLEUNIT[\"degree\",0.0174532925199433]],\n        AXIS[\"geodetic longitude (Lon)\",east,\n            ORDER[2],\n            ANGLEUNIT[\"degree\",0.0174532925199433]],\n    USAGE[\n        SCOPE[\"Horizontal component of 3D system.\"],\n        AREA[\"World.\"],\n        BBOX[-90,-180,90,180]],\n    ID[\"EPSG\",4326]]"), class = "crs"), n_empty = 0L), 
    ling_zo = 4L), sf_column = "geometry", agr = structure(c(NA_integer_, 
NA_integer_, NA_integer_, NA_integer_, NA_integer_, NA_integer_, 
NA_integer_, NA_integer_, NA_integer_, NA_integer_, NA_integer_, 
NA_integer_, NA_integer_, NA_integer_, NA_integer_, NA_integer_
), .Names = c(NA, NA, NA, NA, NA, NA, NA, NA, "state_name", NA, 
NA, NA, NA, "state_abbr", "jurisdiction_type", NA), .Label = c("constant", 
"aggregate", "identity"), class = "factor"), row.names = 2932L, class = c("sf", 
"data.frame"))