[Microsoft Bing Maps] 世界の個別建物フットプリントを R で取得する

2022年後半にかけて、Microsoft (Bing Maps) から世界版の個別建物フットプリントデータのリリースがあり、手元で利用できるようになった。DL して R で使うまでを試してみよう。

#-----------------
Bing Maps is releasing open building footprints around the world. We have detected 999M buildings from Bing Maps imagery between 2014 and 2022 including Maxar, Airbus, and IGN France imagery. The data is freely available for download and use under ODbL.

#-----------------

図1：現時点で公開されている個別建物フットプリントの対象範囲と，Bangkokをカバーする範囲を選択する様子．
残念ながら日本は（まだ？）ない。 公開されました（2023年6月に発見）。
PLATEAU（というより国土地理院の数値地図（国土基本情報）、というより基盤地図情報）と比べてみるとどうなるのだろうか...。

図2：OSMとの比較．
Bangkokの市域でclipして、OSMの建物と比較すると中心部でも結構な違いがみえる。このデータはフットプリントだけで属性が付いていないが空間分布は網羅的で、OSMは建物用途などの属性が付いた建物も多いが疎な地域も多々あるといった具合で、一長一短の様相だ。

ここからコードです。途中、データや図の読み込みでやや時間がかかる（止まったようにみえる）かもしれない。

library(mapview)
library(sf)
library(data.table)
library(dplyr)
library(mapedit)
## library(geojsonR)
## library(geojsonsf)

#---------------------------

# check the current working directory (wd)
getwd()
# make a directory for managing data
dir.create("MSBuildingFootprint")
# change wd to the directory
setwd("./MSBuildingFootprint/")
getwd()
# make directories for saving data
dir.create("Input")
dir.create("Output")

#---------------------------

# read building coverage data
building_coverage <- "https://minedbuildings.blob.core.windows.net/global-buildings/buildings-coverage.geojson" %>%
  sf::st_read()
# check the map of the building coverage 
mapview::mapview(building_coverage)
# select target grids from map viewer (in this example, we select grids covering Bangkok city, Thailand)
## see Fig. 1 （図1）
bc_Bangkok <- mapedit::selectFeatures(building_coverage)
## Check the selected data
## > bc_Bangkok$QuadKey
## [1] "132203132" "132203133" "132203310" "132203311"

#---------------------------

# read url-links data to building footprint
links <- "https://minedbuildings.blob.core.windows.net/global-buildings/dataset-links.csv" %>%
  data.table::fread()
# select url-links of Bangkok
links_Bangkok <- links %>%
  dplyr::filter(Location=="Thailand") %>%
  dplyr::filter(QuadKey %in% bc_Bangkok$QuadKey)

#---------------------------

# loop for downloading, decompressing, and reading data
building_footprint <- NULL
for(i in 1:nrow(links_Bangkok)){
  # download csv.gz file
  download.file(links_Bangkok$Url[i], paste0("Input/d",i,".csv.gz"))
  # unzip csv.gz file
  R.utils::gunzip(paste0("Input/",list.files("./Input")[i]))
  # change file-extension from csv to geojson <see the Note below>
  file.rename(paste0("Input/d",i,".csv"), paste0("Input/d",i,".geojson"))
  # read geojson file
  temp <- sf::st_read(paste0("Input/d",i,".geojson"))
  # bind geojson files
  building_footprint <- dplyr::bind_rows(building_footprint, temp)
  rm(temp)
}

#---------------------------

# write binded building footprint data as geojson 
sf::st_write(building_footprint,"Output/BF_Bangkok.geojson")
# write binded building footprint data as shapfile
sf::st_write(building_footprint,"Output/BF_Bangkok.shp")

Note
geojson のファイルが csv.gz で圧縮されているので、解凍すると csv になる。
中身をみると csv に geojson が入っているあべこべな状態なので、拡張子を変える。これで解決する。