1 目的

The Earth Observation Group (EOG), Payne Institute for Public Policy, Colorado School of Minesが提供する夜間光データを地図上に可視化．夜間光データは経済活動の代理変数として経済学でも利用されている．特に地理（空間）的な経済活動を表すことに適している．
- 夜間光データはラスター型．
  - ラスターデータ：セル（格子状の区画）で地理空間情報を表現するデータ形式．

パッケージ:tidyverse，sf，terra，rnaruralearth，rmapshaper

library(tidyverse) 
library(sf)

#夜間光データの読込など
library(terra) 
#日本地図利用のため
library(rnaturalearth)
#sfオブジェクトの編集，簡略化
library(rmapshaper)

2 1992年の夜間光データの読込と日本周辺の可視化

2.1 夜間光データのダウンロードと読込

夜間光データは地球全体をカバーしているため，ファイルサイズが非常に大きく，ダウンロードに時間がかかる．¹
- そこで，最初にDMSP-F10によって観測された1992年の夜間光データを使用，日本地図に可視化．
  - 解像度は30秒グリッド（約1km）のため，ファイルサイズが相対的に小さい（それでもPCの性能によっては時間がかかることに注意）．

夜間光データのダウンロード

DMSP-OLS　>　Average Visible, Stable Lights, & Cloud Free Coverages　>　F101992
- Average Visible, Stable Lights, & Cloud Free Coverages：雲のない夜の夜間光の平均値を記録．
- F101992：DMSP-F10によって観測された1992年の夜間光データ．

ダウンロードしたZIPファイルは作業フォルダ（ディレクトリ）に納め，ZIPファイルを選択　>　右クリック　>　すべて展開　>　展開先の選択とファイルの選択　>　展開を実行する．展開後は最初にダウンロードしたZIPファイルを削除．
- 展開したファイルにはさらにたZIPファイルが格納されている．そのうち，F101992.v4b_web.stable_lights.avg_vis.tif.gzを「すべて展開」し作業フォルダに格納．
- 利用するのはF101992.v4b_web.stable_lights.avg_vis.tif（夜間光の平均強度を示すデータ）．
- .tif：画像ファイル形式の一つ．
terraパッケージに含まれるrast()でデータを読込．オブジェクト名をnight1992とする．

#夜間光データの読込
night1992 <-
  rast("F101992.v4b_web.stable_lights.avg_vis.tif")

2.2 日本周辺の可視化

日本周辺の地図はNatural Earthを利用．²
- インストールしていない場合は，install.packages("rnaturalearth")，install.packages("rnaturalearthdata")．
- 世界地図（World_map）から日本地図（Japan_map）を切り抜き．

World_map <-
  ne_countries(scale="medium")

日本周辺（極東地域）を示す．東経（x）110～155度，北緯（y）20度～50度までの範囲．
- st_bbox()：境界ボックスの作成．
- crsは世界測地系（4326）． 110, 155

#緯度経度（数値のリスト）
fareast_box <- st_bbox(c(xmin=110, 
                     xmax=155, 
                     ymin=20, 
                     ymax=50), crs=4326)

fareast_boxをst_as_sfc()でポリゴンにする．
st_intersects()：World_mapからfareast_box_sfcに収まる国（ポリゴン）を検索．st_intersects()を利用するには，オブジェクトが（マルチ）ポイント，（マルチ）ライン，（マルチ）ポリゴンである必要．
lengths()>0: 重なるポリゴンを抽出しTRUEに．
オブジェクト名をfareastとする．

#sfcオブジェクト（四角形のポリゴン）の作成
fareast_box_sfc <- fareast_box %>% 
  st_as_sfc()

#重なるポリゴン抽出
fareast <-
  lengths(st_intersects(World_map, fareast_box_sfc))>0

World_mapからfareast==TRUEを満たすメッシュを抽出（filter()）．
- 第1引数.ですべての変数を残す．
- オブジェクト名をfareast_mapとする．
- 異なるパッケージ間で関数名が競合しているため，dplyr::filter()と指示．

fareast_map <- World_map %>%
  dplyr::filter(., fareast==TRUE)

日本周辺の可視化³

coord_sf()で範囲を指定しないと，東西の広い範囲を示してしまう．

ggplot()+
  geom_sf(data=fareast_map, fill="white")+
  coord_sf(xlim = c(110, 155), ylim = c(20, 50))+
  labs(caption="出典：Natural Earth")+
  ggtitle("日本周辺地図")

2.3 日本周辺の夜間光

読み込んだ夜間光データnight1992の観測数は膨大なため，これを可視するのも時間がかかる．
- まず，範囲（先ほど作成ししたfareast_box）内でnight1992を切り抜き，データフレームに変換して可視する．
- crop()：terraパッケージに含まれるcrop()で切り抜くデータを読込．こうすることで観測数を大きく減らせる．
- ext()：切り抜く範囲（数値のリスト）の指示．ここでは，すでに作成したfareast_boxを利用．

night1992 <- night1992 %>% 
  crop(ext(fareast_box))

ggplot()で可視するため，night1992をデータフレーム（as.data.frame()）に変更．
- xy=TRUEを指定することで，座標情報をデータフレームの列として取り込む．
- 夜間光強度の相対値（0～63）を示す3列目の列名が長いため，lightに変更．
- 観測数は900万．

night1992_df <- night1992 %>% 
  as.data.frame(xy=TRUE) %>%
  rename(light=3)

日本周辺の夜間光データの可視化

geom_raster()：夜間光強度（night1992_dfのlight列）の可視化．
tran=sqrtで夜間光強度の平方根を求める．平方根変換により暗いエリアの微妙な差違も可視．

ggplot() +
  geom_raster(data=night1992_df, 
              aes(x=x, y=y, fill=light))+
  scale_fill_viridis_c(option="G", trans="sqrt")+
  geom_sf(data=fareast_map, fill=NA, color="white")+  
  theme_void()+
  coord_sf(xlim=c(110, 155), ylim=c(20, 50)) +
  labs(fill="夜間光強度",
       caption="出典：Natural Earth, Earth Observation Group")+
  ggtitle("日本周辺夜間光（1992年）")

3 2020年の夜間光データの読込と1都3県の可視化

3.1 夜間光データのダウンロードと読込

次にVIIRS（2012年以降）が観測した夜間光データを使用，1都3県（2020年3月）に可視化．
- 解像度は15秒グリッド（約500m）のため，ファイルサイズは大きくなった．
- ただし，年平均だけではなく，月平均のデータ，地球全体（ファイルサイズが大きい）だけではなく，一部の地域だけのデータがダウンロード可能に⇒ファイルサイズが小さい（それでも時間は相当かかる）．
- 利用にはログインが必要に．

夜間光データのダウンロード

EOG　>　Monthly Cloud-free DNB Composite　>　Download Tiled　>　2020　>　202003　>　vcmslcfg　>　SVDNB_npp_20200301-20200331_75N060E_vcmslcfg_v10_c202007042300.tgz （612M）　
- Monthly Cloud-free DNB Composite：雲のない夜の夜間光の月平均値を記録．
- 75N060E：東京を含む日本の夜間光データ．東経（x）60～180度，北緯（y）0度～60度までの範囲のため日本全体を含んでいる．
- 612M：612メガバイト．

ダウンロードしたZIPファイルは作業フォルダ（ディレクトリ）に納め，ZIPファイルを選択　>　右クリック　>　すべて展開　>　展開先の選択とファイルの選択　>　展開を実行する．展開後は最初にダウンロードしたZIPファイルを削除．
- 利用するのはSVDNB_npp_20200301-20200331_75N060E_vcmslcfg_v10_c202007042300.avg_rade9h.tif（夜間光の平均強度を示すデータ）．
- オブジェクト名をnight203とする．

# 夜間光データの読み込み
night203 <- 
  rast("SVDNB_npp_20200301-20200331_75N060E_vcmslcfg_v10_c202007042300.avg_rade9h.tif")

3.2 1都3県（東京都離島除く）の可視化

1都3県の行政区域のシェープファイルを国土数値情報からダウンロード．⁴

#行政区域の読込
Tokyo_map <- read_sf("N03-20220101_13_GML/N03-22_13_220101.shp")

Saitama_map <-
  read_sf("N03-20220101_11_GML/N03-22_11_220101.shp")

Chiba_map <-
  read_sf("N03-20220101_12_GML/N03-22_12_220101.shp")

Kanagagawa_map <-
  read_sf("N03-20220101_14_GML/N03-22_14_220101.shp")

東京都の離島を除き，1都3県を合併．
- bind_rows()：列名の並び順が異なっても，揃えて結合．データは下に加えられてゆく．

#離島の市町村ダミーの作成
Tokyo_map <- Tokyo_map %>%
  mutate(island=
           ifelse(N03_007>=13361 & N03_007<=13421, 
                  1, 0)) 

#離島の市町村を除いたシェープファイルの作成
Tokyo_map <- Tokyo_map %>%
  subset(island!=1)

#1都3県合併
Greater_Tokyo <-
  bind_rows(Saitama_map, Chiba_map, 
        Tokyo_map, Kanagagawa_map)

このままだと，可視化した際の境界線は市区町村単位．境界線を県単位にする（見やすい？）．
- rmapshaperパッケージのms_dissolve()を使用．
- ms_dissolve()：N03_001列に県名が記載されているので，県単位で境界を統合．

#県単位で境界を統合
Greater_Tokyo <- Greater_Tokyo %>% 
  ms_dissolve(field="N03_001")

crsは夜間光データ（nightGTY203）に揃えて世界測地系（4326）．

#CRSをJSD2011→WGS84（4326）
Greater_Tokyo <- Greater_Tokyo %>% 
  st_transform(crs=4326)

1都3県の可視化

ggplot()+
  geom_sf(data=Greater_Tokyo, fill="white")+
  labs(caption="出典：国土数値情報")+
  ggtitle("1都3県")

3.3 1都3県の夜間光

1都3県を含む範囲（GreaterTY_box）でnightGTY203を切り抜き，データフレームに変換して可視する．
st_bbox()：境界ボックスの作成．
mask()：ポリゴンに沿って不要な部分をNAに変更．
vect()：Greater_Tokyo（sfオブジェクト）をterraの SpatVectorに変換．

#緯度経度（数値のリスト） 
GreaterTY_box <- Greater_Tokyo %>% 
  st_bbox()  

#夜間光データを東京都にあわせて切り抜き
nightGTY203  <- night203 %>% 
  crop(ext(GreaterTY_box)) 

#東京都の境界線に合わせて不要な部分を削除
nightGTY203  <- nightGTY203 %>% 
  mask(vect(Greater_Tokyo))

#ggpolt()で可視するためデータフレームに変更
nightGTY203_df <- nightGTY203 %>% 
  as.data.frame(xy=TRUE)%>%
  rename(light=3)

1都3県の夜間光データの可視化

夜間光強度はnW/㎠/sr（ナノワット／平方センチメートル／ステラジアン）．範囲は0以上の正の値であり，絶対的な値．

ggplot()+
  geom_tile(data=nightGTY203_df, 
            aes(x=x, y=y, fill=light))+
  scale_fill_viridis_c(option="G", trans="sqrt")+ 
  geom_sf(data=Greater_Tokyo, fill=NA, color="white")+  
  theme_void()+ 
  labs(fill="夜間光強度（nW/㎠/sr）",
       caption="出典：国土数値情報, Earth Observation Group")+
  ggtitle("1都3県夜間光（2020年3月）")

4 応用：変化率の可視化

2020年3月と4月の夜間光強度の変化率を計算し，可視化．
2020年4月の夜間光データをダウンロードし，1都3県で切り抜き，データフレームに変換し，オブジェクト名をnightGTY204_dfとする．
- nightGTY203_dfと合併し（オブジェクト名：night_data），変化率（change_rate）を求める．

nightGTY203_df <- nightGTY203_df %>%
  rename(light_203=light) 

nightGTY204_df <- nightGTY204_df %>%
  rename(light_204=light) 

night_data <-
  left_join(nightGTY203_df, nightGTY204_df, 
            by=c("x", "y"))

night_data <- night_data %>%
  mutate(change_rate=((light_204-light_203)/light_203)*100)

変化率の記述統計を示し，おおよその分布を把握．

夜間光強度変化率の記述統計量
	平均	標準偏差	最小値	最大値	観察数
変化率（％）	0.90	19.09	-91.85	1244.68	77646

記述統計量を参考にカテゴリー化．⁵

night_data <- night_data %>%
  mutate(change_rate=case_when(
      change_rate>100 ~"(100, 1250)",
      change_rate>50 & change_rate<=100 ~ "(50, 100]",
      change_rate>10 & change_rate<=50 ~ "(10, 50]",
      change_rate>0 & change_rate<=10 ~ "(0, 10]",
      change_rate>-10 & change_rate<=0 ~ "(-10, 0]",
      change_rate>-50 & change_rate<=-10 ~ "(-50, -10]",
      change_rate>-100 & change_rate<=-50 ~ "(-100, -50]"))

factor()を用いて凡例の順番を指定．

night_data <- night_data %>%
  mutate(change_rate=factor(change_rate,
                            levels=c("(100, 1250)",
                                     "(50, 100]",
                                     "(10, 50]", 
                                     "(0, 10]",
                                     "(-10, 0]",
                                     "(-50, -10]",
                                     "(-100, -50]")))

geom_sf()を用いて可視するため，データフレームのnight_dataをsfオブジェクトに変換．

night_data <- night_data %>% 
  st_as_sf(coords=c("x", "y"), crs=4326)

夜間光強度変化率の可視化

scale_color_brewer()：カテゴリーデータ（離散値）の色分け．

ggplot() +
  geom_sf(data=night_data, aes(color=change_rate)) +
  scale_color_brewer(palette="RdBu")+
  geom_sf(data=Greater_Tokyo, fill=NA)+
  theme_void()+ 
  labs(color="夜間光強度変化率（％）",
       caption="出典：国土数値情報, Earth Observation Group") +
  ggtitle("1都3県夜間光（2020年3月～4月）")

Rによる地理空間データの可視化

チュートリアルホーム

アイコン

Font Awesome 6.7.2

生成AIの使用について

　作成過程において，ChatGPTを使用してコードの最適化に関する助言を得た．サービスを利用した後，内容を確認し，必要に応じて編集を行った．

Google Earth Engineを用いて，必要な範囲だけを取得可能な模様．Qitta←検索欄に検索ワード「Google Earth Engine　夜間光」を入力して検索．神奈川大学経済学部浦沢聡士研究室もこのサービスを利用して夜間光データを可視化．↩︎
Natural Earthの利用方法の詳細はNatural Earthの地理空間データの可視化参照．↩︎
地図の国境線については，学術誌（Elsevier，Wiley）では次のような注意書きが求められる場合がある．「Map lines delineate study areas and do not necessarily depict accepted national boundaries」．↩︎
市区町村データによる大学・大学院卒者割合の可視化参照↩︎
カテゴリー化についてはコミュニティバス路線上の人口密度の可視化を参照．↩︎

夜間光ラスターデータの可視化

日本周辺，1都3県の場合

岩田真一郎（　神奈川大学経済学部）

2025-02-25作成，2025-03-02改定

1 目的

2 1992年の夜間光データの読込と日本周辺の可視化

2.1 夜間光データのダウンロードと読込

2.2 日本周辺の可視化

2.3 日本周辺の夜間光

3 2020年の夜間光データの読込と1都3県の可視化

3.1 夜間光データのダウンロードと読込

3.2 1都3県（東京都離島除く）の可視化

3.3 1都3県の夜間光

4 応用：変化率の可視化

夜間光ラスターデータの可視化

日本周辺，1都3県の場合

岩田真一郎（ 神奈川大学経済学部）

2025-02-25作成，2025-03-02改定

1 目的

2 1992年の夜間光データの読込と日本周辺の可視化

2.1 夜間光データのダウンロードと読込

2.2 日本周辺の可視化

2.3 日本周辺の夜間光

3 2020年の夜間光データの読込と1都3県の可視化

3.1 夜間光データのダウンロードと読込

3.2 1都3県（東京都離島除く）の可視化

3.3 1都3県の夜間光

4 応用：変化率の可視化

　岩田真一郎（　神奈川大学経済学部）