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

Question

Facebook 构建了它所谓的 relative wealth index适用于全局 > 1900 万个微型区域(2.4 公里网格单元)。他们已经 shared the data ( zip ) 在一个 csv 文件中，该文件列出了四元键 ID、纬度/经度(我认为是磁贴单元格的左上角)和磁贴的索引值。它看起来像这样:

在他们的 technical paper ，他们注意到这些 2.4 公里的网格单元对应于 Bing 瓦片级别 14。
我没有与 Bing tiles 一起工作过前。 a) 创建或访问覆盖多边形(例如，肯尼亚)的 2.4 平铺网格和 b) 将财富指数值从 csv 连接到此网格 shapefile 的最佳方法是什么？我想要一个具有此财富指数属性的网格多边形，我可以在以后的分析中使用它，该分析通过网格单元从栅格中提取信息。
到目前为止我所知道/认为我所知道的:

sf::st_make_grid()会创建一个网格，但我认为它不会是 Bing 网格。

像{ rosm这样的包将绘制 bing 瓷砖，但这不是我正在寻找的。

人们已经创建了接受四键输入并返回左上角坐标的函数，例如，https://gis.stackexchange.com/a/359636/22560 .我不确定我能用这个做什么，如果有的话。

[从 gis.stackexchange.com 移动问题]
编辑 1:RWI csv 文件不再包含 quadkey，但您可以使用上面链接的 python 包来计算它。有一个有用的教程 here .

Answer 1

这是墨西哥的一个例子，但这是调整 csv 读取的问题。似乎网格对齐得很好(见最后一个图)但是，slippymath 是错误的数据是指单元格中心而不是左上角。结果肯定会更快，但似乎足够快(墨西哥是较大的国家之一)。在第一个位中，我探索在第二个实际读取数据中创建一个网格(这种情况下缩放 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...

创建于 2021-09-30 由 reprex package (v2.0.1)