将有关四键的信息加入平铺 shapefile

Question

Facebook为全球超过 1900 万个微区域（2.4 公里网格单元）构建了所谓的相对财富指数。他们在 csv 文件中共享数据( zip )，其中列出了四键 ID、纬度/经度（我认为是图块单元格的左上角）以及图块的索引值。它看起来像这样：

在他们的技术论文中，他们指出这些 2.4 公里的网格单元对应于 Bing 磁贴级别 14。

我以前没有使用过Bing 磁贴。a) 创建或访问覆盖多边形（例如肯尼亚）的 2.4 平铺网格和 b) 将财富指数值从 csv 连接到此网格形状文件的最佳方法是什么？我想要一个具有此财富指数属性的网格多边形，我可以在将来的分析中使用它，通过网格单元从栅格中提取信息。

到目前为止我所知道/认为我所知道的：

• sf::st_make_grid()会创建一个网格，但我不认为它会是Bing网格。
• 像 { } 这样的包rosm将绘制 bing 磁贴，但这并不是我想要的。
• 人们已经创建了接受四键输入并返回左上角坐标的函数，例如https://gis.stackexchange.com/a/359636/22560。我不确定我能用这个做什么（如果有的话）。

[从 gis.stackexchange.com 移动问题]

编辑 1：RWI csv 文件不再包含四键，但您可以使用上面链接的 python 包来计算它。这里有一个有用的教程。

Answer 1

这是墨西哥的示例，但这是调整 csv 读取的问题。看起来网格对齐得很好（参见最后一张图），但是滑移数学是错误的，数据指的是细胞中心而不是左上角。当然，结果可能会更快，但似乎足够快（墨西哥是较大的国家之一）。在第一部分中，我探索在第二部分中实际读取数据时创建网格（本例放大 4）。请注意，网格的尺寸需要固定，因为一个是规则的，而另一个是不规则的。这会导致st_as_sf.

require(stars)
#> Loading required package: stars
#> Loading required package: abind
#> Loading required package: sf
#> Linking to GEOS 3.8.0, GDAL 3.0.4, PROJ 6.3.1
require(ggplot2)
#> Loading required package: ggplot2
require(tidyverse)
#> Loading required package: tidyverse
z<-4
d<-data.frame(quadkey=apply(as.matrix(expand.grid(0:3,0:3,0:3,0:3)),1, paste0, collapse=''), val=runif(n=(2^z)^2))
lons<-slippymath::tilenum_to_lonlat(0:(2^z),1,z)$lon
lats<-slippymath::tilenum_to_lonlat(1,0:(2^z),z)$lat
l<-lapply(strsplit(d$quadkey,''), as.numeric)
d$x<-unlist(lapply(lapply(lapply(l,'%%',2), '*',2^((z-1):0) ), sum))
d$y<-unlist(lapply(lapply(lapply(l,'%/%',2), '*',2^((z-1):0) ), sum))
m<-matrix(d %>% arrange(x,y) %>% pull(val), ncol=2^z)
grd<-st_as_stars(list(m=m),dimensions=st_dimensions(x=lons, y=lats))
st_crs(grd)=4326
worldMap = rnaturalearth::ne_countries(scale = "medium", returnclass = "sf")
ggplot()+geom_stars(data=grd)+geom_sf(data=worldMap)

# now with the data
f<-tempfile(fileext = '.zip')
download.file('https://data.humdata.org/dataset/76f2a2ea-ba50-40f5-b79c-db95d668b843/resource/bff723a4-6b55-4c51-8790-6176a774e13c/download/relative-wealth-index-april-2021.zip',f)
#unzip(f, list=T)
d<-read.csv(unz(f,'relative-wealth-index-april-2021/MEX_relative_wealth_index.csv'))
z<-14
l<-lapply(strsplit(formatC(d$quadkey, width=z, flag='0', digits = z),''), as.numeric)
d$x<-unlist(lapply(lapply(lapply(l,'%%',2), '*',2^((z-1):0) ), sum))
d$y<-unlist(lapply(lapply(lapply(l,'%/%',2), '*',2^((z-1):0) ), sum))
head(d)
#>        quadkey latitude  longitude    rwi error    x    y
#> 1 2.331030e+12 18.63583  -97.98706 -0.395 0.495 3732 7328
#> 2 2.313220e+12 24.35711 -100.42603 -0.045 0.396 3621 7048
#> 3 2.330113e+12 19.46659 -101.37085 -0.239 0.391 3578 7288
#> 4 2.331032e+12 17.08829  -97.83325 -0.893 0.513 3739 7402
#> 5 2.331030e+12 18.17717  -98.03101  0.053 0.493 3730 7350
#> 6 2.331000e+12 21.21770 -100.29419 -0.065 0.450 3627 7203
lons<-slippymath::tilenum_to_lonlat(min(d$x):(max(d$x)+1),1,z)$lon
lats<-slippymath::tilenum_to_lonlat(1,min(d$y):(max(d$y)+1),z)$lat
require(raster)
#> Loading required package: raster
#> Loading required package: sp
#> 
#> Attaching package: 'raster'
#> The following object is masked from 'package:dplyr':
#> 
#>     select
#> The following object is masked from 'package:tidyr':
#> 
#>     extract
m<-(as.matrix(rasterFromXYZ(d[,c('x','y','rwi')])))
m<-t(m[nrow(m):1,])
grd<-st_as_stars(list(rwi=m),dimensions=st_dimensions(x=lons, y=lats))
st_crs(grd)=4326
# manipulate the dimensions to fix them
dm<-st_dimensions(grd)
dm$x$delta<-NA
dm$x$offset<-NA
ll<-list(start=head(lons,-1), end=tail(lons,-1))
class(ll)<-'intervals'
dm$x$values<-ll
st_dimensions(grd)<-dm
worldMap = rnaturalearth::ne_countries(scale = "medium", returnclass = "sf")
ggplot()+geom_stars(data=grd)+geom_sf(data=worldMap, fill=NA)+coord_sf(xlim=range(lons), ylim=range(lats))

s<-st_as_sf(d,coords = c('longitude','latitude'), crs=4326)
ggplot()+geom_stars(data=grd)+geom_sf(data=worldMap, fill=NA)+geom_sf(data=s, aes(fill=rwi), shape=21)+coord_sf(xlim=-100+0:1, ylim=20+0:1)

st_as_sf(grd)
#> Simple feature collection with 77083 features and 1 field
#> Geometry type: POLYGON
#> Dimension:     XY
#> Bounding box:  xmin: -117.1362 ymin: 14.51978 xmax: -86.7041 ymax: 32.73184
#> Geodetic CRS:  WGS 84
#> First 10 features:
#>       rwi                       geometry
#> 1   0.458 POLYGON ((-114.7632 32.7318...
#> 2   0.642 POLYGON ((-114.7412 32.7318...
#> 3   0.048 POLYGON ((-114.7632 32.7133...
#> 4   0.437 POLYGON ((-114.7412 32.7133...
#> 5  -0.158 POLYGON ((-114.8071 32.6948...
#> 6  -0.031 POLYGON ((-114.7852 32.6948...
#> 7   0.425 POLYGON ((-114.7632 32.6948...
#> 8  -0.070 POLYGON ((-115.5762 32.6763...
#> 9  -0.250 POLYGON ((-115.5542 32.6763...
#> 10  0.371 POLYGON ((-115.5322 32.6763...

^{由reprex 包(v2.0.1)于 2021-09-30 创建}