前処理に使用するパッケージを読み込みます．未インストールの場合は，インストールを済ませておいてください．インストールの過程で「Do you want to install from sources the packages which need compilation? (Yes/no/cancel)」というメッセージが表示された場合には，とりあえず「no」を選択することをお勧めします．
library(tidyverse)
library(sf)
library(openxlsx)
library(janitor)
library(data.table)
library(osmdata)
library(terra)
library(stringi)動作環境によっては，sfパッケージの読み込み時にエラーが出る場合があります．その原因は，sfパッケージが依存しているRcppのバージョンが古いことが多いようです．もしエラーが出る場合には，以下のコードを実行して見てください（筆者のPCではエラーが出なかったので，コメント化してあります）．
# #Rcppパッケージをインストールする
# install.packages("Rcpp")
# #Rcppパッケージを読み込み
# library(Rcpp)取得元：鉄道時系列データ
##鉄道時系列（全国・駅データ）
rail_jp <- sf::read_sf(dsn="N05-20_Station2.shp",
options="ENCODING=CP932") %>%
#WGS84/UTM54Nに投影変換
sf::st_transform(crs=sf::st_crs(32654)) %>%
#現存する駅のレコードだけに絞る
dplyr::filter(N05_005e=="9999") %>%
#ID変数を追加
dplyr::mutate(SID=row_number()) %>%
#必要な変数のみ残す
dplyr::select(SID,N05_002,N05_003,N05_011) %>%
#変数名を分かりやすく変更
dplyr::rename(lname=N05_002,
company=N05_003,
sname=N05_011)
#加工したデータを書き出し
sf::write_sf(obj=rail_jp,
dsn="output/rail_jp.geojson",
delete_dsn=TRUE)取得元：ニュータウン
##ニュータウン（福岡県）
newTown_fukuoka <- sf::read_sf(dsn="P26-13_40.shp",
options="ENCODING=CP932") %>%
#投影法を定義
sf::st_set_crs(value=sf::st_crs(4612)) %>%
#WGS84/UTM54Nに投影変換
sf::st_transform(crs=sf::st_crs(32654)) %>%
#ID変数を追加
dplyr::mutate(TID=row_number()) %>%
#必要な変数のみ残す
dplyr::select(P26_001,P26_002,P26_004,P26_010) %>%
#変数名を分かりやすく変更
dplyr::rename(pref=P26_001,
city=P26_002,
ntname=P26_004,
year_c=P26_010)
#加工したデータを書き出し
sf::write_sf(obj=newTown_fukuoka,
dsn="output/newTown_fukuoka.geojson",
delete_dsn=TRUE)
#加工したデータを書き出し（shapefile形式で）
sf::write_sf(obj=newTown_fukuoka,
dsn="output/newTown_fukuoka.shp",
layer_options="ENCODING=CP932",
delete_dsn=TRUE)取得元：統計地理情報システム
##小地域（福岡市）
chome_fukuoka_c <- sf::read_sf(dsn="h27ka40.shp",
options="ENCODING=CP932") %>%
#WGS84/UTM54Nに投影変換
sf::st_transform(crs=sf::st_crs(32654)) %>%
#変数GST_NAME（市区町村名）が「福岡市」に等しいレコードのみ残す
dplyr::filter(GST_NAME=="福岡市") %>%
#小地域名がNAでない陸地ポリゴンのみに絞る
dplyr::filter(!is.na(MOJI)&!(MOJI%in%c("水面","博多湾"))) %>%
#ID変数を追加
dplyr::mutate(BID=row_number()) %>%
#必要な変数のみ残す
dplyr::select(KEN_NAME,GST_NAME,CSS_NAME,MOJI)
#加工したデータを書き出し
sf::write_sf(obj=chome_fukuoka_c,
dsn="output/chome_fukuoka_c.geojson",
delete_dsn=TRUE)取得元：統計地理情報システム
##500mメッシュ（福岡市）
mesh_fukuoka_0500 <- sf::read_sf(dsn="MESH05030.shp",
options="ENCODING=CP932") %>%
#WGS84/UTM54Nに投影変換
sf::st_transform(crs=sf::st_crs(32654)) %>%
#メッシュIDのみ残す
dplyr::select(KEY_CODE)
#メッシュの重心を作成
mesh_fukuoka_0500_cent <- sf::st_centroid(x=mesh_fukuoka_0500) %>%
#重心と小地域ポリゴンの共通部分を取る
sf::st_intersection(x=.,y=chome_fukuoka_c) %>%
#小地域ポリゴンの属性がついた重心のみ残す
dplyr::filter(!is.na(MOJI))
mesh_fukuoka_0500_c <- mesh_fukuoka_0500 %>%
#福岡市内のメッシュ（＝小地域ポリゴンの属性がついた重心を持つポリゴン）のみ残す
dplyr::filter(KEY_CODE%in%mesh_fukuoka_0500_cent$KEY_CODE)
#加工したデータを書き出し
sf::write_sf(obj=mesh_fukuoka_0500_c,
dsn="output/mesh_fukuoka_0500_c.geojson",
delete_dsn=TRUE)取得元：統計地理情報システム
#csvファイルの読み込み
census <- read.csv(file="tblT000847H5030.txt",
fileEncoding="CP932",
#1行目は飛ばす
skip=1) %>%
#列名をクレンジング
janitor::clean_names(case="old_janitor")
#細かい列名の修正
colnames(census) <- c("KEY_CODE","HTKSYORI","HTKSAKI","GASSAN",colnames(census)[5:length(colnames(census))])
colnames(census) <- gsub("x_","",colnames(census))
census_eld <- census %>%
#高齢夫婦のみの一般世帯数が秘匿でないレコードのみ残す
dplyr::filter(高齢夫婦のみの一般世帯数!="*") %>%
#不要な列を削除
dplyr::select(-c(HTKSYORI,HTKSAKI,GASSAN)) %>%
#変数型を数値に変換
dplyr::mutate(across(where(is.character),as.numeric))
#データ書き出し
openxlsx::write.xlsx(x=census_eld,
file="output/census_2015.xlsx")#「Overpass query unavailable without internet」エラーが出たら下の行を実行
assign("has_internet_via_proxy",TRUE,environment(curl::has_internet))
#日本・福岡県・福岡市に対応するバウンディング・ボックスを取得する
shop_fuk <- osmdata::getbb(place_name="Fukuoka-shi Fukuoka-ken Japan") %>%
#OSMから空間データにアクセス
osmdata::opq() %>%
#キーがshopの空間データを取得
osmdata::add_osm_feature(key="shop") %>%
#sf形式に変換
osmdata::osmdata_sf()
#list形式で返された店舗の空間データから，ポイントデータに該当する要素を取り出す．
shop_fuk_point <- shop_fuk[["osm_points"]] %>%
#施設種別に関する変数shopがNAでないレコードのみ残す
dplyr::filter(!is.na(shop)) %>%
#ID・施設名・施設種別の変数のみ残す
dplyr::select(osm_id,name,shop)
#list形式で返された店舗の空間データから，ポリゴンデータに該当する要素を取り出す．
shop_fuk_polygon <- shop_fuk[["osm_polygons"]] %>%
#施設種別に関する変数shopがNAでないレコードのみ残す
dplyr::filter(!is.na(shop)) %>%
#ID・施設名・施設種別の変数のみ残す
dplyr::select(osm_id,name,shop) %>%
#ポリゴンの重心を取る
sf::st_centroid()
#取り出した要素を行方向に結合
shop_fuk_p <- rbind(shop_fuk_point,shop_fuk_polygon) %>%
#WGS84/UTM54Nに投影変換
sf::st_transform(crs=sf::st_crs(32654))
#加工したデータを書き出し
sf::write_sf(obj=(shop_fuk_p %>% dplyr::select(-name)),
dsn="output/shop_fuk_p.geojson",
delete_dsn=TRUE)
#データをテキストに変換
shop_fuk_p_lonlat <- shop_fuk_p %>%
#WGS84に投影変換
st_transform(crs=sf::st_crs(4326)) %>%
#ポイントデータのX座標（経度）・Y座標（緯度）を取得し，変数として追加する
dplyr::mutate(lon=sf::st_coordinates(.)[,"X"],
lat=sf::st_coordinates(.)[,"Y"]) %>%
#geometry属性を削除する
sf::st_drop_geometry()
#テキストを書き出し
write.csv(x=shop_fuk_p_lonlat,
#書き出す際のファイル名を指定
file="output/fuk_p_lonlat.csv",
#行名は書き出さない
row.names=FALSE,
#書き出す際の文字コードはUTF-8
fileEncoding="UTF-8")取得元：A harmonized global nighttime light dataset 1992–2018
#全世界の夜間光データを読み込み
ntl_09 <- terra::rast("ntl/Harmonized_DN_NTL_2009_calDMSP.tif")
##シリア
ntl_09_sy <- ntl_09 %>%
#シリア国内の緯度経度範囲でデータを切り出し
terra::crop(y=ext(c(35.0,42.3,32.3,37.4)))
#切り出したデータを書き出し
terra::writeRaster(x=ntl_09_sy,
filename="output/ntl_09_sy.tif",
overwrite=TRUE)取得元：A harmonized global nighttime light dataset 1992–2018
#全世界の夜間光データを読み込み
ntl_18 <- terra::rast("ntl/Harmonized_DN_NTL_2018_simVIIRS.tif")
##シリア
ntl_18_sy <- ntl_18 %>%
#シリア国内の緯度経度範囲でデータを切り出し
terra::crop(y=ext(c(35.0,42.3,32.3,37.4)))
#切り出したデータを書き出し
terra::writeRaster(x=ntl_18_sy,
filename="output/ntl_18_sy.tif",
overwrite=TRUE)取得元：地球地図全球版
lc_13 <- terra::rast("gm_lc_v3.tif")
##ベトナム
lc_13_vn <- lc_13 %>%
#ベトナム国内の緯度経度範囲でデータを切り出し
terra::crop(y=ext(c(102,110,8,24)))
#切り出したデータを書き出し
terra::writeRaster(x=lc_13_vn,
filename="output/lc_13_vn.tif",
overwrite=TRUE)取得元：地球地図全球版
el <- terra::rast("gm_el_v2_1_4.tif")
##ベトナム
el_vn <- el %>%
#ベトナム国内の緯度経度範囲でデータを切り出し
terra::crop(y=ext(c(102,110,8,24))) %>%
#標高が-100m以上8000m以下のピクセルのみ残す
terra::clamp(lower=-100,upper=8000,values=TRUE)
#切り出したデータを書き出し
terra::writeRaster(x=el_vn,
filename="output/el_vn.tif",
overwrite=TRUE)取得元：DIVA-GIS
#シリアの行政区域のshpを読み込み
ab_sy <- sf::read_sf(dsn="SYR_adm2.shp") %>%
#必要な変数のみ残す
dplyr::select(NAME_1,NAME_2) %>%
#行政区域IDを作成
dplyr::mutate(ID=row_number())
#データを書き出し
sf::write_sf(obj=ab_sy,
dsn="output/ab_sy.geojson",
delete_dsn=TRUE)取得元：DIVA-GIS
#ベトナムの行政区域のshpを読み込み
ab_vn <- sf::read_sf(dsn="VNM_adm2.shp") %>%
#必要な変数のみ残す
dplyr::select(NAME_1,NAME_2) %>%
#ベトナム語の声調記号を削除する
dplyr::mutate(across(starts_with("NAME"),~stri_trans_general(str=.,id="Latin-ASCII")))
#北部に属す県のリスト
north <- c("Bac Giang",
"Bac Kan",
"Cao Bang",
"Ha Giang",
"Lang Son",
"Phu Tho",
"Quang Ninh",
"Thai Nguyen",
"Tuyen Quang",
"Bac Ninh",
"Ha Nam",
"Hai Duong",
"Hung Yen",
"Nam Dinh",
"Ninh Binh",
"Thai Binh",
"Vinh Phuc",
"Ha Noi",
"Hai Phong")
ab_vn_n <- ab_vn %>%
#北部地方リストに含まれる県のみ取り出し
dplyr::filter(NAME_1%in%north) %>%
#県IDを作成
dplyr::mutate(ID=row_number())
#データを書き出し
sf::write_sf(obj=ab_vn_n,
dsn="output/ab_vn_n.geojson",
delete_dsn=TRUE)