加入有关四键的信息以平铺 shapefile
Join info about quadkeys to tile shapefile
Facebook 在一个 csv 文件中构建了它所谓的 relative wealth index for >19M micro regions (2.4km grid cells) around the world. They've shared the data (zip)，其中列出了四元键 ID lat/long（我认为它是图块单元格的左上角），以及图块的索引值。它看起来像这样：
在他们的 technical paper 中，他们注意到这些 2.4 公里的网格单元对应于 Bing 瓦片级别 14。
我之前没有与 Bing tiles 合作过。 a) 创建或访问覆盖多边形（例如肯尼亚）的 2.4 瓦格网格和 b) 将来自 csv 的财富指数值连接到此网格 shapefile 的最佳方法是什么？我想要一个具有此财富指数属性的网格多边形，我可以在未来的分析中使用它，通过网格单元从栅格中提取信息。
目前我know/think所知道的：
sf::st_make_grid() 会创建 a 网格，但我不认为它会是 the Bing网格.- 像 {
rosm} 这样的软件包将绘制 bing 个图块，但这并不是我想要的。
- 人们已经创建了接受四键输入和 return 左上角坐标的函数，例如 https://gis.stackexchange.com/a/359636/22560。我不确定我能用它做什么。
[从 gis.stackexchange.com 移动问题]
编辑 1：RWI csv 文件不再包含 quadkey，但您可以使用上面链接的 python 包来计算它。有一个有用的教程 here。
这是墨西哥的示例，但这是调整 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...
由 reprex package (v2.0.1)
于 2021-09-30 创建
Facebook 在一个 csv 文件中构建了它所谓的 relative wealth index for >19M micro regions (2.4km grid cells) around the world. They've shared the data (zip)，其中列出了四元键 ID lat/long（我认为它是图块单元格的左上角），以及图块的索引值。它看起来像这样：
在他们的 technical paper 中，他们注意到这些 2.4 公里的网格单元对应于 Bing 瓦片级别 14。
我之前没有与 Bing tiles 合作过。 a) 创建或访问覆盖多边形（例如肯尼亚）的 2.4 瓦格网格和 b) 将来自 csv 的财富指数值连接到此网格 shapefile 的最佳方法是什么？我想要一个具有此财富指数属性的网格多边形，我可以在未来的分析中使用它，通过网格单元从栅格中提取信息。
目前我know/think所知道的：
sf::st_make_grid()会创建 a 网格，但我不认为它会是 the Bing网格.- 像 {
rosm} 这样的软件包将绘制 bing 个图块，但这并不是我想要的。 - 人们已经创建了接受四键输入和 return 左上角坐标的函数，例如 https://gis.stackexchange.com/a/359636/22560。我不确定我能用它做什么。
[从 gis.stackexchange.com 移动问题]
编辑 1：RWI csv 文件不再包含 quadkey，但您可以使用上面链接的 python 包来计算它。有一个有用的教程 here。
这是墨西哥的示例，但这是调整 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...
由 reprex package (v2.0.1)
于 2021-09-30 创建