Rで地図を描く方法
研究用に（あるいは他の用途でも）地図を描画する際、著作権の問題や、正確性の問題を考慮する必要がある。以下では、国土交通省・国土数値情報を用いて、自力で描画できるようにする方法を紹介する。なお、R言語を想定しているが、Pythonでも（ほぼ）同様にできる。
Required packages
install.packages("tidyverse", repos="http://cran.rstudio.com/")install.packages("sf", repos="http://cran.rstudio.com/")library(tidyverse)library(sf)日本全体の地図
日本地図（県境あり）を描く最もシンプルで効率的な方法はnaturalearthのデータを使ったやり方。国土地理院のウェブサイトからAdmin 1 states and provincesのセクションで必要な.zipデータがダウンロードできる。解凍して、R fileと同じフォルダにおいておく。
https://www.naturalearthdata.com/downloads/10m-cultural-vectors/
states01 <- sf::read_sf(".//ne_10m_admin_1_states_provinces") states01#データを確認## Simple feature collection with 4596 features and 121 fields## Geometry type: MULTIPOLYGON## Dimension: XY## Bounding box: xmin: -180 ymin: -90 xmax: 180 ymax: 83.6341## Geodetic CRS: WGS 84## # A tibble: 4,596 × 122## featurecla scalerank adm1_code diss_me iso_3166_2 wikipedia iso_a2 adm0_sr## <chr> <int> <chr> <int> <chr> <chr> <chr> <int>## 1 Admin-1 stat… 3 ARG-1309 1309 AR-E <NA> AR 1## 2 Admin-1 stat… 6 URY-8 8 UY-PA <NA> UY 1## 3 Admin-1 stat… 2 IDN-1185 1185 ID-KI <NA> ID 5## 4 Admin-1 stat… 6 MYS-1186 1186 MY-12 <NA> MY 5## 5 Admin-1 stat… 3 CHL-2694 2694 CL-AP <NA> CL 1## 6 Admin-1 stat… 4 BOL-1936 1936 BO-L <NA> BO 1## 7 Admin-1 stat… 4 BOL-1937 1937 BO-O <NA> BO 1## 8 Admin-1 stat… 3 CHL-2693 2693 CL-TA <NA> CL 1## 9 Admin-1 stat… 4 BOL-1939 1939 BO-P <NA> BO 1## 10 Admin-1 stat… 3 CHL-2695 2695 CL-AN <NA> CL 1## # ℹ 4,586 more rows## # ℹ 114 more variables: name <chr>, name_alt <chr>, name_local <chr>,## # type <chr>, type_en <chr>, code_local <chr>, code_hasc <chr>, note <chr>,## # hasc_maybe <chr>, region <chr>, region_cod <chr>, provnum_ne <int>,## # gadm_level <int>, check_me <int>, datarank <int>, abbrev <chr>,## # postal <chr>, area_sqkm <int>, sameascity <int>, labelrank <int>,## # name_len <int>, mapcolor9 <int>, mapcolor13 <int>, fips <chr>, …japan01 <- states01 %>% dplyr::filter(adm0_a3 =='JPN')#日本地図に必要なデータをjapan01に格納japan01#データを確認## Simple feature collection with 47 features and 121 fields## Geometry type: MULTIPOLYGON## Dimension: XY## Bounding box: xmin: 122.9382 ymin: 24.2121 xmax: 153.9856 ymax: 45.52041## Geodetic CRS: WGS 84## # A tibble: 47 × 122## featurecla scalerank adm1_code diss_me iso_3166_2 wikipedia iso_a2 adm0_sr## * <chr> <int> <chr> <int> <chr> <chr> <chr> <int>## 1 Admin-1 stat… 2 JPN-3501 3501 JP-46 <NA> JP 5## 2 Admin-1 stat… 6 JPN-1835 1835 JP-44 <NA> JP 1## 3 Admin-1 stat… 6 JPN-1829 1829 JP-40 <NA> JP 1## 4 Admin-1 stat… 6 JPN-1827 1827 JP-41 <NA> JP 1## 5 Admin-1 stat… 2 JPN-3500 3500 JP-42 <NA> JP 3## 6 Admin-1 stat… 6 JPN-1830 1830 JP-43 <NA> JP 1## 7 Admin-1 stat… 6 JPN-1831 1831 JP-45 <NA> JP 1## 8 Admin-1 stat… 6 JPN-1836 1836 JP-36 <NA> JP 1## 9 Admin-1 stat… 6 JPN-1833 1833 JP-37 <NA> JP 1## 10 Admin-1 stat… 6 JPN-1832 1832 JP-38 <NA> JP 4## # ℹ 37 more rows## # ℹ 114 more variables: name <chr>, name_alt <chr>, name_local <chr>,## # type <chr>, type_en <chr>, code_local <chr>, code_hasc <chr>, note <chr>,## # hasc_maybe <chr>, region <chr>, region_cod <chr>, provnum_ne <int>,## # gadm_level <int>, check_me <int>, datarank <int>, abbrev <chr>,## # postal <chr>, area_sqkm <int>, sameascity <int>, labelrank <int>,## # name_len <int>, mapcolor9 <int>, mapcolor13 <int>, fips <chr>, …#ggplotとgeom_sfで描画japanmap <- ggplot(data=japan01)+geom_sf()+theme_void()+ theme_void()japanmap
ケーススタディ1 九州全体の地図
九州の地図は、国土交通省・国土数値情報ダウンロードサイト（ここ、新規タブで開く）のデータを使って描画するのが良い。「九州」のファイル（GML）をダウンロード。解凍したフォルダ（GML）を、Rファイルと同じフォルダに入れておく。
#ggplotとgeom_sfで九州地図を描画kyushumap <- ggplot(data =kyushu)+geom_sf()+theme_void()kyushumap
特定の地域（ポイント）を追加
ここでは、九州地図に椎葉村と鹿児島県いちき串木野市の位置をポイントで表示してみる。位置情報（緯度・経度）をここで取得。
#椎葉と串木野ポイント指定した地図shiibakushikinomap <- ggplot(data =kyushu)+geom_sf()+ geom_sf(data = shiiba, color ="black", size =3)+ geom_sf(data = kushikino, color ="black", size =3)+ theme_void()#ポイントの色と大きさはcolorとsizeで調整。shiibakushikinomap
ケーススタディ2 琉球列島の地図
琉球列島の場合、鹿児島県のうち大島郡のみ（喜界島から与論まで）の地図情報と沖縄県の地図情報を合算させる必要がある。一番シンプルなやり方を以下で紹介する。
九州の場合と同様、国土交通省・国土数値情報ダウンロードサイト（ここ、新規タブで開く）のデータを使って描画する。「鹿児島県」「沖縄県」のファイル（GML）をダウンロード。解凍したフォルダ（GML）を、Rファイルと同じフォルダに入れておく。
奄美データの取得と描画
まず鹿児島県のデータを加工して奄美データを抜き取る。
#まず鹿児島の地図情報を読み込むkagoshima <- st_read("N03-20210101_46_GML")## Reading layer `N03-21_46_210101' from data source## `/Users/michinorishimoji2023/Documents/4_stats/Geo/N03-20210101_46_GML'## using driver `ESRI Shapefile'## Simple feature collection with 8136 features and 5 fields## Geometry type: POLYGON## Dimension: XY## Bounding box: xmin: 128.3954 ymin: 27.01868 xmax: 131.2054 ymax: 32.3106## Geodetic CRS: JGD2011str(kagoshima)## Classes'sf' and'data.frame':8136 obs. of6 variables:## $ N03_001 : chr"鹿児島県""鹿児島県""鹿児島県""鹿児島県" ...## $ N03_002 : chr NA NA NA NA ...## $ N03_003 : chr NA NA NA NA ...## $ N03_004 : chr"鹿児島市""鹿児島市""鹿児島市""鹿児島市" ...## $ N03_007 : chr"46201""46201""46201""46201" ...## $ geometry:sfc_POLYGON of length8136; first list element: List of1## ..$ : num [1:7,1:2]131131131131131 ...## ..- attr(*,"class")= chr [1:3]"XY""POLYGON""sfg"## - attr(*,"sf_column")= chr"geometry"## - attr(*,"agr")= Factor w/3 levels"constant","aggregate",..: NA NA NA NA NA## ..- attr(*,"names")= chr [1:5]"N03_001""N03_002""N03_003""N03_004" ...#そのうち奄美市、大島郡のデータだけを抜き取るamami <- kagoshima[kagoshima$N03_007%in%c("46523","46524","46525","46222","46527","46222","46529","46530","46531","46532","46533","46534","46535"), ]#このデータで奄美列島のみの地図を描画できるamamimap <- ggplot(data =amami)+geom_sf()+ theme_void()amamimap
沖縄県のデータと合わせて琉球全体を描画
#奄美データの地図と沖縄県データの地図を重ねてryukyumapとして描画ryukyumap <- ggplot(data =amami)+geom_sf()+ geom_sf(data=okinawa)+theme_void()ryukyumap
宮古の地図
#宮古諸島の地図を描画miyakomap <- ggplot(data =miyako)+geom_sf()+theme_void()miyakomap
