GIS Lecture 3: R で GIS を使う。
| OS | MacOSX 10.12 |
|---|---|
| Fink | |
| R | 3.5 |
3.1. Rの概要
略
3.2. R のインストール
$ fink install r-base35 $ fink install cran-rgdal-r35
3.3. R の使い方
すると、以下のようなメッセージが表示されます。
統計局から、04101 仙台市青葉区の世界測地系緯度経度・Shape形式をダウンロードし、解凍しておきます。
$ R-3.5
WARNING: ignoring environment value of R_HOME R version 3.5.2 (2018-12-20) -- "Eggshell Igloo" Copyright (C) 2018 The R Foundation for Statistical Computing Platform: x86_64-apple-darwin16.7.0 (64-bit) R は、自由なソフトウェアであり、「完全に無保証」です。 一定の条件に従えば、自由にこれを再配布することができます。 配布条件の詳細に関しては、'license()' あるいは 'licence()' と入力してください。 R は多くの貢献者による共同プロジェクトです。 詳しくは 'contributors()' と入力してください。 また、R や R のパッケージを出版物で引用する際の形式については 'citation()' と入力してください。 'demo()' と入力すればデモをみることができます。 'help()' とすればオンラインヘルプが出ます。 'help.start()' で HTML ブラウザによるヘルプがみられます。 'q()' と入力すれば R を終了します。
> library(rgdal)
要求されたパッケージ sp をロード中です rgdal: version: 1.3-6, (SVN revision 773) Geospatial Data Abstraction Library extensions to R successfully loaded Loaded GDAL runtime: GDAL 2.3.1, released 2018/06/22 Path to GDAL shared files: /sw/share/gdal GDAL binary built with GEOS: TRUE Loaded PROJ.4 runtime: Rel. 4.9.2, 08 September 2015, [PJ_VERSION: 492] Path to PROJ.4 shared files: (autodetected) Linking to sp version: 1.3-1
では、ダウンロードした仙台市宮城野区の国勢調査データを読み込んで見ましょう。 (注: dsn がフォルダ名、layer がファイル名 (拡張子を除いたもの)です。)
> lnd<-readOGR(dsn=".",layer="h27ka04101")
OGR data source with driver: ESRI Shapefile Source: "/Users/baba/gisdata", layer: "h27ka04101" with 286 features It has 35 fields Integer64 fields read as strings: JINKO SETAI >
ここに書かれている情報は重要です。JINKOとSETAIを文字列と認識しているようです。
> lnd$JINKO
[1] 697 798 2535 2550 1793 2080 1078 0 220 127 47 2148 874 434 597 [16] 698 2585 1338 40 31 1986 2529 271 508 581 653 708 605 3763 1257 [31] 716 1076 139 1912 1504 1383 1282 1362 933 2144 554 337 504 532 932 [46] 429 934 221 961 683 115 679 1148 780 2863 2991 1862 1601 1938 0 [61] 1063 6017 1820 2016 2162 1267 1523 1952 33 490 1160 2419 1393 476 514 [76] 1140 1698 109 1411 174 1427 1212 1259 1378 2294 2769 2486 1020 1123 940 [91] 714 1282 620 97 1573 1311 1201 542 932 2968 667 364 119 1544 579 [106] 693 1819 1191 2048 1048 1189 665 1882 1382 1103 884 368 1830 1003 1254 [121] 965 1324 1257 1490 0 29 2275 2097 867 1400 557 1187 1592 1119 729 [136] 200 976 1248 1536 445 2536 3050 723 778 175 1974 1441 1253 1431 2902 [151] 325 1493 1116 2504 1315 1869 1236 848 18 45 73 600 1145 1498 2375 [166] 1569 1249 1104 1385 1476 1681 315 1406 1449 1087 1501 891 760 392 1114 [181] 1366 698 2537 1241 1102 1069 1602 841 998 2269 1813 905 355 1970 1566 [196] 1837 2483 1633 777 817 571 1807 623 1099 1144 2252 1266 615 1497 441 [211] 2582 1300 444 239 951 1129 1698 31 895 1002 1023 732 569 1712 1345 [226] 1308 1548 1256 451 92 2791 955 766 727 269 622 1448 381 978 843 [241] 522 1452 487 417 924 777 1072 655 1201 886 771 566 337 460 1085 [256] 615 533 18 648 756 465 571 577 199 517 0 4 0 174 1011 [271] 324 1388 542 42 147 348 484 582 1 224 652 1275 614 718 297 [286] 506 268 Levels: 0 1 1002 1003 1011 1020 1023 1048 1063 1069 1072 1076 1078 ... 998
JINKO を整数値に変換します。
> lnd$JINKO <- as.numeric(as.character(lnd$JINKO)) > lnd$JINKO
[1] 697 798 2535 2550 1793 2080 1078 0 220 127 47 2148 874 434 597 [16] 698 2585 1338 40 31 1986 2529 271 508 581 653 708 605 3763 1257 [31] 716 1076 139 1912 1504 1383 1282 1362 933 2144 554 337 504 532 932 [46] 429 934 221 961 683 115 679 1148 780 2863 2991 1862 1601 1938 0 [61] 1063 6017 1820 2016 2162 1267 1523 1952 33 490 1160 2419 1393 476 514 [76] 1140 1698 109 1411 174 1427 1212 1259 1378 2294 2769 2486 1020 1123 940 [91] 714 1282 620 97 1573 1311 1201 542 932 2968 667 364 119 1544 579 [106] 693 1819 1191 2048 1048 1189 665 1882 1382 1103 884 368 1830 1003 1254 [121] 965 1324 1257 1490 0 29 2275 2097 867 1400 557 1187 1592 1119 729 [136] 200 976 1248 1536 445 2536 3050 723 778 175 1974 1441 1253 1431 2902 [151] 325 1493 1116 2504 1315 1869 1236 848 18 45 73 600 1145 1498 2375 [166] 1569 1249 1104 1385 1476 1681 315 1406 1449 1087 1501 891 760 392 1114 [181] 1366 698 2537 1241 1102 1069 1602 841 998 2269 1813 905 355 1970 1566 [196] 1837 2483 1633 777 817 571 1807 623 1099 1144 2252 1266 615 1497 441 [211] 2582 1300 444 239 951 1129 1698 31 895 1002 1023 732 569 1712 1345 [226] 1308 1548 1256 451 92 2791 955 766 727 269 622 1448 381 978 843 [241] 522 1452 487 417 924 777 1072 655 1201 886 771 566 337 460 1085 [256] 615 533 18 648 756 465 571 577 199 517 0 4 0 174 1011 [271] 324 1388 542 42 147 348 484 582 1 224 652 1275 614 718 297 [286] 506
上二つをよく見比べてください。左寄せは文字列、右寄せは数値です。
> head(lnd@data, n=2)
KEY_CODE PREF CITY S_AREA PREF_NAME CITY_NAME S_NAME KIGO_E HCODE 0 041010010 04 101 001000 宮城県 青葉区 青葉町8101 1 041010020 04 101 002000 宮城県 青葉区 あけぼの町 8101 AREA PERIMETER H27KAxx_ H27KAxx_ID KEN KEN_NAME SITYO_NAME GST_NAME 0 132477.0 1660.062 5063 5062 04 宮城県 仙台市 1 98653.1 1768.070 2634 2633 04 宮城県 仙台市 CSS_NAME KIHON1 DUMMY1 KIHON2 KEYCODE1 KEYCODE2 AREA_MAX_F KIGO_D N_KEN 0 青葉区 0010 - 00 101001000 1010010 M 1 青葉区 0020 - 00 101002000 1010020 M N_CITY KIGO_I MOJI KBSUM JINKO SETAI X_CODE Y_CODE KCODE1 0 青葉町 18 697 402 140.8648 38.28277 0010-00 1 あけぼの町 14 798 438 140.8643 38.28817 0020-00
では、地図としてプロットして見ます。
> plot(lnd) > plot(lnd, col = "lightgrey") # 仙台市青葉区の人口をプロット > sel <- lnd$JINKO > 200 > plot(lnd[ sel, ], col = "turquoise", add = TRUE) # 人口200より上をターコイズで上書き
一見、成功したように見えます。次の QGIS で作った地図と比較すると、、
めでたし、めでたし。