将空间数据分解为更小的网格单元

问题描述 投票：0回答：1

我正在寻求帮助，将一些空间数据分解为更小的网格单元。

structure(list(GEOID = c("17031010100", "17031010200"), NAME.x = c("Census Tract 101, Cook County, Illinois", 
"Census Tract 102, Cook County, Illinois"), median_inc = c(30708, 
35932), pop = c(4835, 8830), share25_edu = c(17.1124828532236, 
21.5430163706026), pop_poverty = c(1719, 2139), share_poverty = c(35.5532574974147, 
24.2242355605889), median_age = c(30.3, 32.6), geometry = structure(list(
    structure(list(list(structure(c(-87.677765, -87.678074, -87.675257, 
    -87.673257, -87.672251, -87.671657, -87.668657, -87.668657, 
    -87.666357, -87.6650920979757, -87.6662895604484, -87.6665966344484, 
    -87.67245, -87.677357, -87.677765, 42.02303, 42.023011, 42.020531, 
    42.019231, 42.019294, 42.019331, 42.019431, 42.019231, 42.019231, 
    42.0193213501446, 42.0223445264575, 42.0231197815191, 42.023031, 
    42.02303, 42.02303), dim = c(15L, 2L)))), class = c("XY", 
    "MULTIPOLYGON", "sfg")), structure(list(list(structure(c(-87.676457, 
    -87.679457, -87.680757, -87.683157, -87.684958, -87.684557, 
    -87.683957, -87.683357, -87.680457, -87.678457, -87.674857, 
    -87.672657, -87.670557, -87.670657, -87.670657, -87.670757, 
    -87.671557, -87.671657, -87.672251, -87.673257, -87.676457, 
    42.019131, 42.019531, 42.019531, 42.019431, 42.019431, 42.016431, 
    42.014131, 42.012331, 42.012531, 42.012631, 42.012731, 42.012731, 
    42.012731, 42.013931, 42.016031, 42.018331, 42.018031, 42.019331, 
    42.019294, 42.019231, 42.019131), dim = c(21L, 2L)))), class = c("XY", 
    "MULTIPOLYGON", "sfg"))), class = c("sfc_MULTIPOLYGON", "sfc"
), precision = 0, bbox = structure(c(xmin = -87.684958, ymin = 42.012331, 
xmax = -87.6650920979757, ymax = 42.0231197815191), 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        MEMBER[\"World Geodetic System 1984 (G2139)\"],\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)), row.names = c(NA, 
-2L), tigris = "tract", sf_column = "geometry", agr = structure(c(GEOID = NA_integer_, 
NAME.x = NA_integer_, median_inc = NA_integer_, pop = NA_integer_, 
share25_edu = NA_integer_, pop_poverty = NA_integer_, share_poverty = NA_integer_, 
median_age = NA_integer_), class = "factor", levels = c("constant", 
"aggregate", "identity")), class = c("sf", "tbl_df", "tbl", "data.frame"
))

variable_tracts 包含伊利诺伊州库克县所有人口普查区的某些变量的数据。我对位于芝加哥的那些感兴趣，我是通过

找到的 
chi_tracts <- variables_tract |> 
  st_filter(chi_map)

哪里

chi_map_comm <- read_sf("https://raw.githubusercontent.com/thisisdaryn/data/master/geo/chicago/Comm_Areas.geojson") 
chi_map <- st_union(chi_map_comm)

这剩下大约 800 个人口普查区。

**我想做什么**

我想将 chi_tracts 中的所有人口普查区域分解为更小的网格单元。我想获取每个网格单元的几何形状。

我担心这个过程可能需要多长时间 - 800 个人口普查区太多了。

我对如何解决这个问题感到非常困惑，所以任何帮助将不胜感激！如果格式不正确，我深表歉意 - 我在这个论坛上不是很有经验。

我尝试了以下从chatgpt无耻地获取的代码，该代码运行了很长时间没有效果：


split_into_grid <- function(tract, cellsize) {
  grid <- st_make_grid(tract, cellsize = 0.01, what = "polygons", square = TRUE) |> 
    st_as_sf()  
 grid <- st_intersection(grid, tract)  
 grid <- grid |> 
    mutate(
      median_inc = tract$median_inc,         
      share25_edu = tract$share25_edu,       
      pop = tract$pop / nrow(grid), 
      share_poverty = tract$share_poverty,             
      median_age = tract$median_age     
    )
  
  return(grid) 
}

cell_size <- 0.01  


chi_grid <- chi_tracts |> 
  group_split() |> 
  map_dfr(~split_into_grid(.x, cell_size))
r rstudio
1个回答
0
投票

您可以在您的区域 {sf} 对象和网格之间仅使用 

st_intersection()
。

gr <- sf::st_make_grid(variables_tract, cellsize = c(0.01, 0.01)) |>
  sf::st_as_sf()
gr$id <- seq(nrow(gr))

tmap::qtm(variables_tract) +
  tmap::qtm(gr, fill = "id", fill_alpha = 0.4)

当您与它们相交时，两个 sf 数据帧中的所有变量都将被复制到结果对象中：


p <- variables_tract |>
  sf::st_intersection(gr)

p
#> Simple feature collection with 6 features and 9 fields
#> Geometry type: POLYGON
#> Dimension:     XY
#> Bounding box:  xmin: -87.68496 ymin: 42.01233 xmax: -87.66509 ymax: 42.02312
#> Geodetic CRS:  WGS 84
#> # A tibble: 6 × 10
#>   GEOID NAME.x median_inc   pop share25_edu pop_poverty share_poverty median_age
#> * <chr> <chr>       <dbl> <dbl>       <dbl>       <dbl>         <dbl>      <dbl>
#> 1 1703… Censu…      30708  4835        17.1        1719          35.6       30.3
#> 2 1703… Censu…      35932  8830        21.5        2139          24.2       32.6
#> 3 1703… Censu…      30708  4835        17.1        1719          35.6       30.3
#> 4 1703… Censu…      35932  8830        21.5        2139          24.2       32.6
#> 5 1703… Censu…      30708  4835        17.1        1719          35.6       30.3
#> 6 1703… Censu…      30708  4835        17.1        1719          35.6       30.3
#> # ℹ 2 more variables: id <int>, geometry <POLYGON [°]>

tmap::qtm(p, fill = "id")

现在您可以按 

id
（来自网格）或来自 Variables_tract 的 
GEOID
进行分组并执行分析。

创建于 2024 年 12 月 13 日，使用 reprex v2.1.1.9000

最新问题
© www.soinside.com 2019 - 2024. All rights reserved.