#r #sf
Вопрос:
Facebook построил так называемый индекс относительного благосостояния для >19 млн микрорегионов (2,4 км ячеек сетки) по всему миру. Они поделились данными (zip) в файле csv, в котором указан идентификатор ключа quad, lat/long (который, как я полагаю, находится в верхнем левом углу ячейки плитки) и значение индекса для плитки. Это выглядит так:
В своем техническом документе они отмечают, что эти ячейки сетки длиной 2,4 км соответствуют 14-му уровню плитки Bing.
Я раньше не работал с плитками Bing. Каков лучший способ а) создать или получить доступ к сетке листов 2.4, которая охватывает полигон (например, Кения), и б) объединить значения индекса благосостояния из csv в этот шейп-файл сетки? Я хотел бы иметь многоугольник сетки с этим атрибутом индекса богатства, который я мог бы использовать в будущем анализе, извлекающем информацию из растра по ячейкам сетки.
Что я знаю/думаю, что знаю до сих пор:
sf::st_make_grid()создал бы сетку, но я не думаю, что это будет сетка Bing.- Такие пакеты, как {
rosm}, будут отображать плитки bing, но это не совсем то, что я ищу. - Люди создали функции, которые принимают ввод квадроключем и возвращают координату верхнего левого угла, например, https://gis.stackexchange.com/a/359636/22560. Я не уверен, что смогу с этим что-нибудь сделать, если вообще что-нибудь смогу.
[перенесенный вопрос из gis.stackexchange.com]
Правка 1: csv-файлы RWI больше не содержат квадрокоптер, но вы можете использовать пакет python, связанный выше, для его вычисления. Здесь есть полезный урок.
Ответ №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 (v2.0.1)
Комментарии:
1. Это исключительное событие, @Bart. Я многому научился из вашего ответа. Как вы можете видеть, существует много отсутствующих значений ячеек, поэтому мне нужно изучить, как интерполировать
rwiзначения. Как только я это сделаю, я хотел бы преобразовать объект starsgrdв шейп-файл. Это последняя часть моего первоначального вопроса. Я думал, что это может бытьst_as_sf(grd, as_points = FALSE, merge = FALSE), но я получаю неверную тусклую ошибку.2. Я обновил ответ, похоже, это проблема,
starsно можно исправить размеры.3. PS интерполяция сетки кажется довольно сложной, по крайней мере, есть регионы, где данные очень скудны, а пространственные корреляции на первый взгляд не кажутся такими высокими. Я думаю, что состоятельные соседи-капюшоны могут граничить совсем с другими. Я бы, по крайней мере, проверил, на каких расстояниях интерполяция все еще дает разумные результаты.
4. Это хороший довод. Похоже, что проблема в Мексике хуже, чем в других наборах данных, но хорошо бы это учесть!