如何使用 r 的 sf 包在 sf 多边形内创建 500m * 500m 的网格？

Question

我正在尝试在 class 为“sf”的纽约萨福克郡边界内创建网格。我将图层命名为“SUFF”。通过st_area(SUFF)得知县城面积为6136105813平方米

所以，我决定创建大小为 500 米 * 500 米的矩形网格。我写了代码： fishnet <- st_make_grid(st_transform(SUFF, crs=st_crs(4326)),cellsize = 500, square = TRUE) %>% st_sf().

然而，我只有一个格子。 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

按照我的理解，像元大小的单位应该和SUFF图层一样，是米，对吗？你介意给我一些指导，告诉我如何使用 st_make_grid() 创建 500m * 500m 的网格吗？

Answer 1

欢迎使用 Whosebug。我在下面的数据中有一个萨福克县。正如@d-j 指出的那样，在这种情况下，边界框是包含所有 Suffolk 的最小正方形。根据县的形状，相对于正方形，您的许多网格可能不会超过 Suffolk，如果您将值插值到 Suffolk 上（例如，在陆地上），则应考虑这一点。使用更大的网格单元尺寸：

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)

以及您可能感兴趣的与萨福克交叉的网格：

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

当你想做一个 st_interpolate_aw loosely like this.

时，考虑这个差异很有用

数据：

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"))

如何使用 r 的 sf 包在 sf 多边形内创建 500m * 500m 的网格？

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

r

sf