最終更新：2025-12-30

1 VIIRS夜間光データ

  • VIIRS夜間光データ（Annual VNL V2）を用いる。

  • このデータは、2012年から2020年にかけての年間グローバル夜間光の時系列データである。一貫した処理方法で生成されており、経年比較が可能である。

  • 日光・月光・雲を除去した後、野火などの外れ値を除外している。

  • V2では12ヶ月の中央値を基準に外れ値を除去し、複数年のデータを使用して統一的な検出閾値を設定することで、年度間の一貫性を確保している。

  • 本ページで用いるデータ：https://ayumu-tanaka.github.io/teaching/NightLightData.zip

2 RでのVIIRS夜間光データの読み込み

  • terraパッケージのrast関数を用いて、GeoTIFF形式のVIIRS夜間光データを読み込む。
library(terra)
## terra 1.8.86
NightLight <- terra::rast("Data_output/NightLightJapan2021.tif")

3 国勢調査のシェープファイル

  • 東京都の市区町村データを用いる。
  • 国勢調査のシェープファイルを用いる。シェープファイルは、e-Statの統計地理情報システムの境界データダウンロードのページからダウンロードできる。
  • 「境界データダウンロード」→「小地域」「国勢調査」→「2015年」→「小地域（町丁・字等）（JGD2000）」→「世界測地系平面直角座標系・Shapefile」→「東京都」→「東京都全域」（世界測地系平面直角座標系・Shapefile）を選択してダウンロードする。

3.1 座標参照系（CRS: Coordinate Reference System）について

  • 世界測地系（WGS 84 / JGD 2000など）には2つの主な座標系がある
項目 緯度経度 平面直角座標系
形状 球面座標（3次元） 平面座標（2次元）
単位 度・分・秒 または 十進度 メートル (m)
原点 地球の中心 各系ごとに異なる
標準名 EPSG:4326 (WGS84) EPSG:6668～6677（日本）
用途 グローバル、汎用 地形図、測量、計画
  • 今回は、「平面直角座標系」のshpファイルをダウンロードして用いて、夜間光の衛星データ（緯度経度）と接続する際に、座標変換を行う流れにする。
  • 「平面直角座標系」の利点は、距離の計算が簡単で正確になることである。
  • 「緯度経度」のshpファイルをダウンロードして、用いた場合は、座標変換は不要になる。

3.2 東京都市区町村のシェープファイルの読み込みと表示

3.2.1 全体

library(sf)
## Linking to GEOS 3.13.0, GDAL 3.8.5, PROJ 9.5.1; sf_use_s2() is TRUE
tokyo <- sf::st_read("Data_raw/Kokusei2015/h27ka13.shp")
## Reading layer `h27ka13' from data source 
##   `/Users/ayumu/Documents/GitHub/research/4_GIS/NightLight/Data_raw/Kokusei2015/h27ka13.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 6010 features and 35 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: -392946.9 ymin: -1721444 xmax: 1446129 ymax: -10968.95
## Projected CRS: JGD2000 / Japan Plane Rectangular CS IX
# plotする
plot(sf::st_geometry(tokyo))

3.2.2 島嶼部を除く東京都の地図を表示

tokyo_mainland <- subset(tokyo, !CITY_NAME %in% c("利島村", "大島町", "小笠原村", "御蔵島村", "新島村", "神津島村", "青ヶ島村", "三宅村", "八丈町"))

plot(sf::st_geometry(tokyo_mainland))

4 座標参照システムの確認と統一

  • crs関数を用いて、夜間光データと市区町村データのCRSを確認する。
  • WGS 84（EPSG:4326）であることを確認する。
print(crs(NightLight))
## [1] "GEOGCRS[\"WGS 84\",\n    DATUM[\"World Geodetic System 1984\",\n        ELLIPSOID[\"WGS 84\",6378137,298.257223563,\n            LENGTHUNIT[\"metre\",1,\n                ID[\"EPSG\",9001]]]],\n    PRIMEM[\"Greenwich\",0,\n        ANGLEUNIT[\"degree\",0.0174532925199433,\n            ID[\"EPSG\",9122]]],\n    CS[ellipsoidal,2],\n        AXIS[\"latitude\",north,\n            ORDER[1],\n            ANGLEUNIT[\"degree\",0.0174532925199433,\n                ID[\"EPSG\",9122]]],\n        AXIS[\"longitude\",east,\n            ORDER[2],\n            ANGLEUNIT[\"degree\",0.0174532925199433,\n                ID[\"EPSG\",9122]]]]"
  • 市区町村データのCRSを確認する。
  • JGD2000（EPSG:4612）であることを確認する。
print(st_crs(tokyo_mainland))
## Coordinate Reference System:
##   User input: JGD2000 / Japan Plane Rectangular CS IX 
##   wkt:
## PROJCRS["JGD2000 / Japan Plane Rectangular CS IX",
##     BASEGEOGCRS["JGD2000",
##         DATUM["Japanese Geodetic Datum 2000",
##             ELLIPSOID["GRS 1980",6378137,298.257222101,
##                 LENGTHUNIT["metre",1]]],
##         PRIMEM["Greenwich",0,
##             ANGLEUNIT["degree",0.0174532925199433]],
##         ID["EPSG",4612]],
##     CONVERSION["Japan Plane Rectangular CS zone IX",
##         METHOD["Transverse Mercator",
##             ID["EPSG",9807]],
##         PARAMETER["Latitude of natural origin",36,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8801]],
##         PARAMETER["Longitude of natural origin",139.833333333333,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8802]],
##         PARAMETER["Scale factor at natural origin",0.9999,
##             SCALEUNIT["unity",1],
##             ID["EPSG",8805]],
##         PARAMETER["False easting",0,
##             LENGTHUNIT["metre",1],
##             ID["EPSG",8806]],
##         PARAMETER["False northing",0,
##             LENGTHUNIT["metre",1],
##             ID["EPSG",8807]]],
##     CS[Cartesian,2],
##         AXIS["northing (X)",north,
##             ORDER[1],
##             LENGTHUNIT["metre",1]],
##         AXIS["easting (Y)",east,
##             ORDER[2],
##             LENGTHUNIT["metre",1]],
##     USAGE[
##         SCOPE["Cadastre, engineering survey, topographic mapping (large and medium scale)."],
##         AREA["Japan - onshore - Honshu - Tokyo-to. (Excludes offshore island areas of Tokyo-to covered by Japan Plane Rectangular Coordinate System zones XIV, XVIII and XIX)."],
##         BBOX[29.31,138.4,37.98,141.11]],
##     ID["EPSG",2451]]

4.1 CRSを統一（EPSG:4612）

  • st_transform関数を用いて、市区町村データtokyo_mainlandのCRSを夜間光データNightLightのCRSに変換する。
tokyo_mainland <- sf::st_transform(tokyo_mainland, crs = crs(NightLight))
print(st_crs(tokyo_mainland))
## Coordinate Reference System:
##   User input: GEOGCRS["WGS 84",
##     DATUM["World Geodetic System 1984",
##         ELLIPSOID["WGS 84",6378137,298.257223563,
##             LENGTHUNIT["metre",1,
##                 ID["EPSG",9001]]]],
##     PRIMEM["Greenwich",0,
##         ANGLEUNIT["degree",0.0174532925199433,
##             ID["EPSG",9122]]],
##     CS[ellipsoidal,2],
##         AXIS["latitude",north,
##             ORDER[1],
##             ANGLEUNIT["degree",0.0174532925199433,
##                 ID["EPSG",9122]]],
##         AXIS["longitude",east,
##             ORDER[2],
##             ANGLEUNIT["degree",0.0174532925199433,
##                 ID["EPSG",9122]]]] 
##   wkt:
## GEOGCRS["WGS 84",
##     DATUM["World Geodetic System 1984",
##         ELLIPSOID["WGS 84",6378137,298.257223563,
##             LENGTHUNIT["metre",1,
##                 ID["EPSG",9001]]]],
##     PRIMEM["Greenwich",0,
##         ANGLEUNIT["degree",0.0174532925199433,
##             ID["EPSG",9122]]],
##     CS[ellipsoidal,2],
##         AXIS["latitude",north,
##             ORDER[1],
##             ANGLEUNIT["degree",0.0174532925199433,
##                 ID["EPSG",9122]]],
##         AXIS["longitude",east,
##             ORDER[2],
##             ANGLEUNIT["degree",0.0174532925199433,
##                 ID["EPSG",9122]]]]
  • 夜間光データNightLightのCRSを市区町村データtokyo_mainlandのCRSに変換する場合は、terraパッケージのproject関数を用いた以下のコードを用いることで変換できるはずだが、エラーがでて、変換に失敗した。そのため、上記のように、市区町村データtokyo_mainlandのCRSを夜間光データNightLightのCRSに変換した。
NightLight2 <- terra::project(NightLight, "EPSG:4612")

5 市区町村ごとに夜間光の統計量を抽出

  • terraパッケージのextract関数を用いて、市区町村ごとに夜間光データの統計量を抽出する。
  • 市区町村ポリゴンを使用してマスク処理し、統計量を計算する。
# 統計量の抽出
NightLightStat <- terra::extract(
  NightLight,
  vect(tokyo_mainland),
  fun = c("mean", "median", "sum", "sd"),
  na.rm = TRUE
)

# 変数名の変更
colnames(NightLightStat)[-1] <- c("ntl_mean", "ntl_median", "ntl_sum", "ntl_sd")

# 結果をデータフレームに変換（不要かも）
NightLightStat <- as.data.frame(NightLightStat)

# 結果の表示
head(NightLightStat)
##   ID  ntl_mean ntl_median   ntl_sum   ntl_sd
## 1  1 114.42669  130.31349 343.28007 31.76616
## 2  2  94.03635   94.03635  94.03635      NaN
## 3  3 122.22236  135.11537 611.11181 25.35337
## 4  4  74.27013   74.27013  74.27013      NaN
## 5  5 109.82556  109.82556 109.82556      NaN
## 6  6 145.42760  145.42760 145.42760      NaN

6 データの結合

  • cbind関数を用いて、市区町村データと夜間光の統計量を結合する。
# 市区町村IDと結合
tokyo_with_ntl <- cbind(tokyo_mainland, NightLightStat)

7 地図表示

7.1 plot関数による表示

7.1.1 シンプルなコード

  • plot関数を用いて、東京都の市区町村別夜間光強度の地図を表示する。
plot(tokyo_with_ntl["ntl_mean"],
     main = "Night Light Intensity by Municipality in Tokyo (2021, Mean)")

7.1.2 色分け数を指定して表示するコード

library(RColorBrewer)
# 色分け数
nc <- 7 
# 配色
pal <- brewer.pal(nc, "YlOrRd")
# データの最小値と最大値から breaks を設定
min_val <- min(tokyo_with_ntl$ntl_mean, na.rm = TRUE)
max_val <- max(tokyo_with_ntl$ntl_mean, na.rm = TRUE)
breaks <- seq(min_val, max_val, length.out = nc + 1)

plot(tokyo_with_ntl["ntl_mean"],
     main = "Night Light Intensity by Municipality in Tokyo (2021, Mean)",
     pal = pal,
     breaks = breaks)

7.1.3 等分位数ベースの breaks を設定

  • 等分位数（quantile-based）でbreaksを設定することで、各カテゴリーに均等にサンプルを配分できる。
  • ここでは、quantile関数を用いて、等分位数ベースの breaks を設定し、地図を表示する。
# 色分け数
nc <- 7 
# 配色
pal <- brewer.pal(nc, "YlOrRd")
# 等分位数ベースの  breaks を設定
breaks <- quantile(
  tokyo_with_ntl$ntl_mean,
  probs = seq(0, 1, length.out = nc + 1),
  na.rm = TRUE
)

plot(tokyo_with_ntl["ntl_mean"],
     main = "Night Light Intensity by Municipality in Tokyo (2021, Mean)",
     pal = pal,
     breaks = breaks)

7.2 ggplot2パッケージによる表示

  • ggplot2パッケージを用いて、東京都の市区町村別夜間光強度の地図を表示する。
library(ggplot2)

ggplot(data = tokyo_with_ntl) +
  geom_sf(aes(fill = ntl_mean), alpha = 0.5) +
  scale_fill_viridis_c(option = "plasma", guide = guide_colourbar(alpha = 0.5)) +
  labs(title = "Night Light Intensity by Municipality in Tokyo (2021, Mean)",
       fill = "Mean NTL") +
  theme_minimal()

8 データの保存

8.1 シェープファイルとして保存

  • st_write関数を用いて、東京都の市区町村別夜間光強度データをシェープファイルとして保存する。
sf::st_write(tokyo_with_ntl, "Data_output/tokyo_with_nightlight2021.shp", delete_layer = TRUE)
## Deleting layer `tokyo_with_nightlight2021' using driver `ESRI Shapefile'
## Writing layer `tokyo_with_nightlight2021' to data source 
##   `Data_output/tokyo_with_nightlight2021.shp' using driver `ESRI Shapefile'
## Writing 5609 features with 40 fields and geometry type Polygon.
## Warning in CPL_write_ogr(obj, dsn, layer, driver,
## as.character(dataset_options), : GDAL Message 1: One or several characters
## couldn't be converted correctly from UTF-8 to ISO-8859-1.  This warning will
## not be emitted anymore.

8.2 CSVとして保存（非地理情報）

  • rioパッケージの`export関数を用いて、東京都の市区町村別夜間光強度データをCSVとして保存する。
library(rio)
# 地理情報を除いたデータフレームを作成
tokyo_with_ntl_df <- st_set_geometry(tokyo_with_ntl, NULL)
# CSVとして保存
rio::export(tokyo_with_ntl_df, "Data_output/tokyo_nightlight2021.csv")

9 渋谷区の夜間光強度

  • subset関数を用いて、渋谷区のデータを抽出し、夜間光強度を表示する。
shibuya <- subset(tokyo_with_ntl, CITY_NAME == "渋谷区")

plot(shibuya["ntl_mean"],
     main = "Night Light Intensity in Shibuya Ward (2021, Mean)")

10 市区レベル集計

10.1 集計

  • dplyrパッケージを用いて、市区ごとに夜間光強度の平均値と合計値を集計する。
library(dplyr)
library(sf)

# 市区ごとに集計（実際の市区列名に合わせて変更）
tokyo_city <- tokyo_with_ntl %>%
  group_by(CITY_NAME) %>%  # CITY_NAME = 市区を識別する列
  summarise(
    geometry = st_union(geometry),
    ntl_mean = mean(ntl_mean, na.rm = TRUE),
    ntl_sum = sum(ntl_sum, na.rm = TRUE),
    n_chome = n(),
    .groups = 'drop'
  ) %>%
  st_as_sf()

10.2 英語の市区名の追加

  • dplyrパッケージのmutate関数を用いて、英語の市区名の列を追加する。
tokyo_city <- tokyo_city %>%
  mutate(
    CITY_NAME_EN = case_when(
      CITY_NAME == "千代田区" ~ "Chiyoda",
      CITY_NAME == "中央区" ~ "Chuo",
      CITY_NAME == "港区" ~ "Minato",
      CITY_NAME == "新宿区" ~ "Shinjuku",
      CITY_NAME == "文京区" ~ "Bunkyo",
      CITY_NAME == "台東区" ~ "Taito",
      CITY_NAME == "墨田区" ~ "Sumida",
      CITY_NAME == "江東区" ~ "Koto",
      CITY_NAME == "品川区" ~ "Shinagawa",
      CITY_NAME == "目黒区" ~ "Meguro",
      CITY_NAME == "大田区" ~ "Ota",
      CITY_NAME == "世田谷区" ~ "Setagaya",
      CITY_NAME == "渋谷区" ~ "Shibuya",
      CITY_NAME == "中野区" ~ "Nakano",
      CITY_NAME == "杉並区" ~ "Suginami",
      CITY_NAME == "豊島区" ~ "Toshima",
      CITY_NAME == "北区" ~ "Kita",
      CITY_NAME == "荒川区" ~ "Arakawa",
      CITY_NAME == "板橋区" ~ "Itabashi",
      CITY_NAME == "練馬区" ~ "Nerima",
      CITY_NAME == "足立区" ~ "Adachi",
      CITY_NAME == "葛飾区" ~ "Katsushika",
      CITY_NAME == "江戸川区" ~ "Edogawa",
      CITY_NAME == "三鷹市" ~ "Mitaka",
      CITY_NAME == "武蔵野市" ~ "Musashino",
      CITY_NAME == "府中市" ~ "Fuchu",
      CITY_NAME == "調布市" ~ "Chofu",
      CITY_NAME == "町田市" ~ "Machida",
      CITY_NAME == "小金井市" ~ "Koganei",
      CITY_NAME == "小平市" ~ "Kodaira",
      CITY_NAME == "東村山市" ~ "Higashimurayama",
      CITY_NAME == "国分寺市" ~ "Kokubunji",
      CITY_NAME == "国立市" ~ "Kunitachi",
      CITY_NAME == "昭島市" ~ "Akishima",
      CITY_NAME == "福生市" ~ "Fussa",
      CITY_NAME == "羽村市" ~ "Hamura",
      CITY_NAME == "あきる野市" ~ "Akiruno",
      CITY_NAME == "西東京市" ~ "Nishitokyo",
      CITY_NAME == "多摩市" ~ "Tama",
      CITY_NAME == "奥多摩町" ~ "Okutama",
      CITY_NAME == "日の出町" ~ "Hinode",
      CITY_NAME == "東久留米市" ~ "Higashikurume",
      CITY_NAME == "清瀬市" ~ "Kiyose",
      CITY_NAME == "稲城市" ~ "Inagi",
      CITY_NAME == "立川市" ~ "Tachikawa",
      CITY_NAME == "八王子市" ~ "Hachioji",
      CITY_NAME == "日野市" ~ "Hino",
      CITY_NAME == "東大和市" ~ "Higashiyamato",
      CITY_NAME == "檜原村" ~ "Hinohara",
      CITY_NAME == "武蔵村山市" ~ "Musahimurayama",
      CITY_NAME == "狛江市" ~ "Komae",
    CITY_NAME == "瑞穂町" ~ "Mizuho",
      CITY_NAME == "青梅市" ~ "Ome",
      TRUE ~ "Other"
    )
  )

10.3 地図

library(ggplot2)

tokyo_city <- tokyo_city %>%
  mutate(
    lon = st_coordinates(st_centroid(geometry))[, 1],
    lat = st_coordinates(st_centroid(geometry))[, 2]
  )

ggplot(data = tokyo_city) +
  geom_sf(aes(fill = ntl_mean), alpha = 0.2) +
  geom_text(aes(x = lon, y = lat, label = CITY_NAME_EN),
            size = 2, fontface = "bold", color = "black") +
  scale_fill_viridis_c(option = "plasma", 
                       guide = guide_colourbar(alpha = 0.2)) +
  theme_minimal()

11 東京23区の夜間光強度

  • subset関数を用いて、東京23区のデータを抽出し、夜間光強度を表示する。
# 東京23区のデータを抽出
tokyo_23 <- subset(tokyo_city , CITY_NAME %in% c("
千代田区", "中央区", "港区", "新宿区", "文京区", "台東区", "墨田区", "江東区", "品川区", "目黒区", "大田区", "世田谷区", "渋谷区", "中野区", "杉並区", "豊島区", "北区", "荒川区", "板橋区", "練馬区", "足立区", "葛飾区", "江戸川区"))


library(ggplot2)
ggplot(data = tokyo_23) +
  geom_sf(aes(fill = ntl_mean), alpha = 0.2) +
  geom_text(aes(x = lon, y = lat, label = CITY_NAME_EN),
            size = 2, fontface = "bold", color = "black") +
  scale_fill_viridis_c(option = "plasma", 
                       guide = guide_colourbar(alpha = 0.2)) +
  theme_minimal()

参考文献