前回までで都道府県の輪郭を取得することは出来たのですが、レコード数が多すぎて(500万件近く)白地図作成には適さないので、今回はそのレコード数を絞る作業を行います。
前回取得したデータを検証すると、領域数が7万以上ありました。これらのほとんど島で、しかも島というより岩礁に近いような領域もかなり含まれていてデータを肥大化させています。
そのため、こういった小さい島々を白地図作成対象から除外すればデータはスッキリしそうです。最初に考えた除外方法は、「レコード数（=登録点数）の少ない領域を対象外にする」というもの。
簡単な方法なので早速試しました。登録点数5000件以下の領域を除外すると、残りの領域数は100くらいになり、ちょうど良いと思ったのですが・・・。これだと、日本の島しょ（本州・北海道・九州・四国を除く）の中で一番大きい択捉島が漏れてしまうことが判明。次に大きい国後島も同様で、北方領土の島々は測量が出来ないためかデータ点数が少ないようです。これは北方領土のみを例外対応とすれば済むのですが、他の領域についても面積とデータ点数が比例しているかは分からないので、正攻法で領域の面積を求めることにしました。本来の目的からは脱線している気がしますが、せっかくの機会なので学習を兼ねて・・・。
「面積を求める」といっても、領域の形は円や四角形といった単純なものではありません。簡単に計算出来るのか？と思っていたのですが、Web検索したところ以下サイトに便利な計算式が紹介されているのを見つけました。これは使えます！
|
| 図-1 多角形の例 |
上図の例では、多角形の面積は、
= (30-20-10+5+33+26)÷2
= 32
となります。※反時計回りに点を結ぶと、計算値は正の数値となります。
ただし、この計算を行うには輪郭を形成する各点の緯度・経度を座標値（x,y）に変換しなければなりません。
この作業には結構手間がかかりました・・・。
緯度・経度を座標値（平面直角座標と言うようです）に変換する方法（計算式）は、国土地理院のサイトに掲載されています。
・国土地理院 測量計算サイト:http://surveycalc.gsi.go.jp/sokuchi/main.html
＜緯度・経度から平面直角座標への変換＞
計算式 http://surveycalc.gsi.go.jp/sokuchi/surveycalc/algorithm/bl2xy/bl2xy.htm
上記の国土地理院サイトでは、緯度・経度→平面直角座標の変換計算が出来ますが、計算プログラムは公開されていません。
それでも「きっと何処か他のサイトで紹介されているだろうから、それを参考に作ってみよう！」なんて軽く考えていたのですが、どうやら計算式が最近変わったようで・・・。古い式での計算プログラム掲載のサイトはいくつかあったのですが、新しい計算式でのプログラムは探しても見つかりません。
仕方がないので、国土地理院サイトのやたら複雑な計算式と睨めっこしながら１からプログラミングです。
そして・・・、
なかなか正しい結果が得られず苦労しましたが、何とか出来上がりました。特に嵌ったのが以下の２点
・tan-1(x)、tanh-1(x) といった三角関数の「-1」の部分の解釈
「tan-1とはタンジェントの逆数だろう」と思い、1/tan(x)としてプログラムを組んだ所、正しい計算結果を得られない。いろいろ調べて「tan-1はタンジェントの逆関数である」ということを知りました。昔、習っただろうか？ 三角関数なんて久しく使っていないので・・・。
PHPにはtan()の逆関数atan()、及びtanh()の逆関数atanh()が存在しているので、需要はそれなりにあるのでしょうね。
・国土地理院サイトの計算式に登場する、謎の変数「ρ″」
ρ″とはラジアンから秒への変換係数とのこと。当係数は計算式の分母に使われているので逆に秒からラジアンへの変換を意味します。
しかし、変換元の緯度（φ0）は三角関数で使用するため既にラジアン変換されているはず！ しかも何故に「度」ではなくて「秒」からの変換なのか？？？？
測量ではよく使う係数のようですが、門外漢の私には訳が分からず・・・。
各種資料をあさった末、最新版の理科年表に同様の計算式を発見。こちらの式ではρ″の表記はないので、それに従いこの係数を無視して（つまりρ″=1とみなして）計算した所、上手くいきました！
以下のソースが緯度・経度を平面直角座標に変換する関数です。検証のため、いくつかの地点について国土地理院の計算ページと比較してみましたが、結果は一致したのでおそらく正しいかと思われます。
<?php
//
// 緯度・経度を平面直角座標に変換
//
// 国土地理院サイト（http://surveycalc.gsi.go.jp/sokuchi/main.html）の計算式を適用
//
// $p0:座標原点の緯度・経度 array("Lat"=>latitude,"Lng"=>longitude)
// $p1:対象点の緯度・経度 array("Lat"=>latitude,"Lng"=>longitude)
//
function LatLng_to_XY($p0, $p1) {
$a = 6378137.0; //長半径
$F = 298.257222101; //逆扁平率
$m0 = 0.9999; //平面直角座標系のX軸上における縮尺係数
$n = 1/(2*$F-1);
//緯度経度をラジアン変換
$radLat0 = deg2rad($p0["Lat"]); //座標原点の緯度
$radLng0 = deg2rad($p0["Lng"]); //座標原点の経度
$radLat1 = deg2rad($p1["Lat"]); //変換対象点の緯度
$radLng1 = deg2rad($p1["Lng"]); //変換対象点の経度
$t = sinh(atanh(sin($radLat1))-(2*sqrt($n))/(1+$n)*(atanh((2*sqrt($n))/(1+$n)*sin($radLat1))));
$tbar = sqrt(1+pow($t,2));
$lmdc = cos($radLng1-$radLng0);
$lmds = sin($radLng1-$radLng0);
$xi = atan($t/$lmdc);
$eta = atanh($lmds/$tbar);
$alp1 = 1/2*$n - 2/3*pow($n,2) + 5/16*pow($n,3) + 41/180*pow($n,4) - 127/288*pow($n,5);
$alp2 = 13/48*pow($n,2) - 3/5*pow($n,3) + 557/1440*pow($n,4) + 281/630*pow($n,5);
$alp3 = 61/240*pow($n,3) - 103/140*pow($n,4) + 15061/26880*pow($n,5);
$alp4 = 49561/161280*pow($n,4) - 179/168*pow($n,5);
$alp5 = 34729/80640*pow($n,5);
$A0 = 1 + pow($n,2)/4 + pow($n,4)/64;
$A1 = -3/2 * ($n - pow($n,3)/8 - pow($n,5)/64);
$A2 = 15/16 * (pow($n,2) - pow($n,4)/4);
$A3 = -35/48 * (pow($n,3) - 5/16*pow($n,5));
$A4 = 315/512 * pow($n,4);
$A5 = -693/1280 * pow($n,5);
$rho2d = 1; // 国土地理院サイトの計算式にある「ρ''」。
// ラジアン→秒 変換係数なので本来の値は(180/π)*3600であるが、理科年表の式に当表記はない
// ρ''=１として計算した所、正しいと思われる変換結果が得られた
$Sbar0 = (($m0*$a)/(1+$n))*($A0*$radLat0/$rho2d + $A1*sin(2*$radLat0) + $A2*sin(4*$radLat0) + $A3*sin(6*$radLat0)
+ $A4*sin(8*$radLat0) + $A5*sin(10*$radLat0)); //赤道から座標系原点までの子午線弧長
$Abar = (($m0*$a)/(1+$n))*$A0;
$coordinate = array();
$x = $Abar * ($xi + $alp1*sin(2*$xi)*cosh(2*$eta) + $alp2*sin(4*$xi)*cosh(4*$eta) + $alp3*sin(6*$xi)*cosh(6*$eta)
+ $alp4*sin(8*$xi)*cosh(8*$eta) + $alp5*sin(10*$xi)*cosh(10*$eta)) - $Sbar0;
$y = $Abar * ($eta + $alp1*cos(2*$xi)*sinh(2*$eta) + $alp2*cos(4*$xi)*sinh(4*$eta) + $alp3*cos(6*$xi)*sinh(6*$eta)
+ $alp4*cos(8*$xi)*sinh(8*$eta) + $alp5*cos(10*$xi)*sinh(10*$eta));
//日本の測量ではX軸が北方向、Y軸が東方向となっていて数学座標とはX,Yが逆のため、
//上記のx,y（測量座標）を入れ替えて数学座標として返す
$coordinate["x"] = $y;
$coordinate["y"] = $x;
return $coordinate;
}
?>
当関数は、最後（62～63行め）にxとyの値を入れ替えたものを返り値としています。これは日本の測量では何故かX軸が北方向、Y軸が東方向となっていて数学座標とはx,yの関係が逆のため、求めた座標（測量座標）を数学座標に合わせるための処理です。
次は座標の原点となる位置の指定。日本では19箇所の基準点が「平面直角座標系」として定められていて、座標を求める地点によって適用する座標系原点が異なります。法令はこちら
以下は、対象点の都道府県コード及び緯度・経度を元に適用すべき座標系原点を選択し、その情報を返す関数のソースです。
<?php
// 座標原点とする日本の平面直角座標系を選択
// 平成十四年国土交通省告示第九号 より
// ※北海道などは簡易的な選択を行っているため、本来の座標系とは異なる座標が選択される場合がある
//
// $pref:都道府県コード $lat：緯度 $lng:経度
//
function selectCoordinateOrigin($pref,$lat,$lng) {
$jpc = array();
$jpc[1] = array("33:00:00","129:30:00");
$jpc[2] = array("33:00:00","131:00:00");
$jpc[3] = array("36:00:00","132:10:00");
$jpc[4] = array("33:00:00","133:30:00");
$jpc[5] = array("36:00:00","134:20:00");
$jpc[6] = array("36:00:00","136:00:00");
$jpc[7] = array("36:00:00","137:10:00");
$jpc[8] = array("36:00:00","138:30:00");
$jpc[9] = array("36:00:00","139:50:00");
$jpc[10] = array("40:00:00","140:50:00");
$jpc[11] = array("44:00:00","140:15:00");
$jpc[12] = array("44:00:00","142:15:00");
$jpc[13] = array("44:00:00","144:15:00");
$jpc[14] = array("26:00:00","142:00:00");
$jpc[15] = array("26:00:00","127:30:00");
$jpc[16] = array("26:00:00","124:00:00");
$jpc[17] = array("26:00:00","131:00:00");
$jpc[18] = array("20:00:00","136:00:00");
$jpc[19] = array("26:00:00","154:00:00");
if ($pref==42 || ($pref==46 && ($lat>=27 && $lat<=32 && $lng>=128.3 && $lng<=130))) {
//長崎県、鹿児島県の一部
$gNo = 1;
} else if ($pref>=40 && $pref<=46) { //福岡県、佐賀県、熊本県、大分県、宮崎県、鹿児島県（I系区域を除く)
$gNo = 2;
} else if ($pref==32 || $pref==34 || $pref==35) { //島根県、広島県、山口県
$gNo = 3;
} else if ($pref>=36 && $pref<=39) { //徳島県、香川県、愛媛県、高知県
$gNo = 4;
} else if ($pref==28 || $pref==31 || $pref==33) { //兵庫県、鳥取県、岡山県
$gNo = 5;
} else if ($pref==18 || ($pref>=24 && $pref<=27) || $pref==29 || $pref==30) { //福井県、三重県、京都府、大阪府、奈良県、和歌山県
$gNo = 6;
} else if ($pref==16 || $pref==17 || $pref==21 || $pref==23) { //富山県、石川県、岐阜県、愛知県
$gNo = 7;
} else if ($pref==15 || $pref==19 || $pref==20 || $pref==22) { //新潟県、山梨県、長野県、静岡県
$gNo = 8;
} else if ($pref==13) { //東京都
if ($lat<=28) {
if ($lng<=140.5) {
$gNo = 18;
} else if ($lng>=143) {
$gNo = 19;
} else {
$gNo = 14;
}
} else {
$gNo = 9;
}
} else if ($pref>=7 && $pref<=14) { //福島県、茨城県、栃木県、群馬県、埼玉県、千葉県、神奈川県
$gNo = 9;
} else if ($pref>=2 && $pref<=6) { //青森県、岩手県、宮城県、秋田県、山形県
$gNo = 10;
} else if ($pref==1) { //北海道（本来は市区町村により振り分けるが、簡易的に経度で分けている）
if ($lng<=141) {
$gNo = 11;
} else if ($lng>=143.5) {
$gNo = 13;
} else {
$gNo = 12;
}
} else if ($pref==47) { //沖縄県
if ($lng<=126) {
$gNo = 16;
} else if ($lng>=130) {
$gNo = 17;
} else {
$gNo = 15;
}
} else { //デフォルト値
$gNo = 8;
}
$coordinateorg = array();
$hmsLat = explode(":",$jpc[$gNo][0]);
$hmsLng = explode(":",$jpc[$gNo][1]);
return array("Lat"=> intval($hmsLat[0]) + intval($hmsLat[1])/60 + intval($hmsLat[2])/3600,
"Lng"=> intval($hmsLng[0]) + intval($hmsLng[1])/60 + intval($hmsLng[2])/3600,
"GNo"=> $gNo);
}
?>
基本的には対象点の都道府県によって座標系原点が決まりますが、北海道、東京都、鹿児島県、沖縄県は場所により適用する座標系原点が変わります。うち、東京都、鹿児島県、沖縄県は緯度・経度により原点を決めることになっているので、上のソースでほぼ正しく割り当て出来ますが、北海道に関しては法令では市区町村単位で原点を割り当てています。しかし、北海道本島全体の輪郭を座標化する際に市区町村により原点を分けてしまうのは都合が悪いので、簡易的に経度で適用原点を振り分けることにしました。そのため、一部本来とは異なる原点を選択してしまう場合もありますが（小さな島等）、今回の目的（=面積計算）上は特に影響ないと思われます。
次に、領域の面積の計算結果を収納するテーブルを用意します。
[SQL]
CREATE TABLE IF NOT EXISTS t_PrefectureBorderArea (
PrefectureCode smallint(2) unsigned zerofill NOT NULL, — 都道府県コード
BorderCode int(5) NOT NULL, — 領域コード
Area decimal(20,8) DEFAULT NULL, — 領域面積
ActiveFlag tinyint(2) NOT NULL DEFAULT ‘0’, — 白地図描画の対象とするか否かのフラグ
PRIMARY KEY (PrefectureCode,BorderCode)
) ENGINE=MyISAM DEFAULT CHARSET=utf8;
[/SQL]
今回求める面積で、各領域毎に白地図出力の対象とするか否かを判断するので、その処理用のフラグ（ActiveFlag）も設けました。
以上でようやく準備は整ったので、今回のメインパートである面積計算へ・・・。
<?php
ini_set( 'error_reporting', ~E_WARNING );
ini_set( 'memory_limit', '2048M' );
set_time_limit(0);
$errMsg;
$mysqli = connectDBi($errMsg);
if ($mysqli) {
echo "MySQL接続OK<br>";
} else {
echo "MySQL接続NG<br>".$errMsg;
exit();
}
echo str_pad(" ",4096)."<br>";
for ($PrefCode = 1; $PrefCode <= 47; $PrefCode++) {
echo "Pref=".$PrefCode."<br>";
ob_flush();
flush();
$border = array();
$strSQL = "SELECT PrefectureCode, BorderCode, Avg(Lat) AS AvgLat, Avg(Lng) AS AvgLng";
$strSQL .= " FROM t_PrefectureBorder";
$strSQL .= " WHERE PrefectureCode = ".$PrefCode;
$strSQL .= " GROUP BY BorderCode";
$strSQL .= " ORDER BY BorderCode";
$rst = $mysqli->query($strSQL);
while($col = $rst->fetch_array(MYSQLI_ASSOC)) {
array_push($border,$col);
}
$rst->close();
for ($b=0;$b<count($border);$b++) {
$s = calBorderArea($mysqli,$border[$b]);
echo " BorderCode=".$border[$b]["BorderCode"]." Area=".$s."<br>";
ob_flush();
flush();
$strSQL = "INSERT INTO t_PrefectureBorderArea (PrefectureCode,BorderCode,Area)";
$strSQL .= " VALUES (".$border[$b]["PrefectureCode"].",".$border[$b]["BorderCode"].",".$s.")";
if (!$mysqli->query($strSQL)) {
echo "失敗！". $strSQL ."<br>";
exit();
}
}
}
if ($mysqli) {
$mysqli->close();
}
echo "終了しました。<br>";
ini_restore('error_reporting');
ini_restore('memory_limit');
//領域の面積を計算
function calBorderArea($mysqli, $BorderInfo) {
$strSQL = "SELECT * FROM t_PrefectureBorder";
$strSQL .= " WHERE PrefectureCode = ".$BorderInfo["PrefectureCode"];
$strSQL .= " AND BorderCode = ".$BorderInfo["BorderCode"];
$strSQL .= " ORDER BY PointNo";
$rst = $mysqli->query($strSQL);
$point = array();
$currentP = array();
$formerP = array();
$S = 0;
if ($rst->num_rows > 0) {
$p0 = selectCoordinateOrigin($BorderInfo["PrefectureCode"],$BorderInfo["AvgLat"],$BorderInfo["AvgLng"]); //座標系原点を取得
while($col = $rst->fetch_array(MYSQLI_ASSOC)) {
$currentP = LatLng_to_XY($p0,$col);
if (isset($formerP["x"]) && isset($formerP["y"])) {
$S += ($formerP["x"]-$currentP["x"])*($formerP["y"]+$currentP["y"]) / 2;
}
$formerP = $currentP;
}
}
$rst->close();
return $S;
}
//MySQLへ接続
function connectDBi(&$err) {
$MySQL_SERVER = "サーバー名";
$MySQL_USER = "ユーザー名";
$MySQL_PASSWORD = "パスワード";
$MySQL_DBNAME = "データベース名";
$err = "";
$mysqli = new mysqli($MySQL_SERVER, $MySQL_USER, $MySQL_PASSWORD, $MySQL_DBNAME);
if ($mysqli->connect_errno) {
$err = "データベース接続に失敗しました。";
} else {
//文字化け対策
if (!$mysqli->set_charset("utf8")) {
$err = "文字コードセットに失敗しました。";
}
}
if (strlen($err) > 0) {
return false;
} else {
return $mysqli;
}
}
?>
まずは領域情報を配列$borderへ丸ごと読み込んでいます（20行目～31行目）。レコードセットを開いたまま次の作業に移るよりも、一旦変数に読み込んでしまう方が段違いに速度が速いためです。
元テーブル（t_PrefectureBorder）にはない緯度・経度の平均値も取得していますが（21行目）、これは前述した座標系原点の選択時に使用します。
面積の計算は、関数「calBorderArea」内で行っています（54～82行目)。
まずは、前記した関数「selectCoordinateOrigin」を呼び出して、対象領域の都道府県コード及び平均緯度・経度をもとに原点となる座標系を決定し（69行目）、各点の緯度・経度を座標数値に変換しています（71行目）。
面積の計算は73行目の箇所。本記事の冒頭部に記載した計算式では、一番最後（n番目）の点まで到達後に最初の点に戻り、
（Xn-X1) (Yn+Y1)
の結果も加算する必要があるのですが、対象のデータ（テーブル：t_PrefectureBorder）は各領域とも最終レコードに最初のレコードと同一点を収納しているので、今回はこの処理は不要です。
また、対象データは基本的※に反時計回りに収納されているので計算結果は正の数になります。
※領域内に隣の都道府県の飛び地があるケースでは時計回りで収納されていますが、この領域の面積は引き算の対象になりますので、結果が負の数となるのはむしろ好都合です。
そして、面積計算の結果を領域毎にテーブルt_PrefectureBorderAreaへ書き込んでいます。
計算結果を検証するため、都道府県単位を集計し、国土地理院による「平成25年全国都道府県市区町村別面積調」の数値と比較してみました（表-1）。
| 面積の単位：km2 | |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
|
|
||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
| ※面積は国土交通省国土地理院『平成25年全国都道府県市区町村別面積調』の数値 | |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
一番誤差が大きいのが大阪府で0.53％。この他、東京都、神奈川県など埋立地の多い都府県の面積が実際よりも大きめの数値となっています。一方、実際よりも小さい数値となったのは長崎県で-0.33%。これらの誤差は潮汐など海岸線のデータ取得方法の違いによるものなのかもしれません。
全体的には実際の面積よりも若干小さめな計算結果となっていますが、大きなズレはなく概ね正しく計算出来ていると判断して良さそうです。
また、主要４島（本州、北海道、九州、四国）を除く主な島についても、国土地理院の数値と比較してみました（表-2）。
| 面積の単位：km2 | |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
|
|
||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
| ※面積は国土交通省国土地理院『平成25年全国都道府県市区町村別面積調』の数値 ※面積20km2を超える島しょについて掲載 |
|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
上表でセルに色がついているのは、実際の面積と計算値との誤差が１％以上となった島です。
特に差異が大きかったのが対馬で36％以上。これは明らかに計算誤差の範囲を超えているので調べてみた所、対馬は元の領域データが下図のように３分割されていました。ストリートビューで確認すると、橋が架かっていて確かに領域としては分かれています。一方、国土地理院調査では３領域を合わせて対馬の面積として扱っているようなので、これが差異の原因。３領域の計算結果を足すと誤差は0.35％程度となるので、計算自体は問題ないと思われます。
 図-2 対馬 ※領域が３つに分割されています（分かりやすいよう色分けしています） |
他の小豆島等も同様で公表の面積と計算結果のズレは領域が分かれていることによるものでした。
この程度の精度であれば、今回の目的（領域を白地図作成時に対象とするか否かの判断）には使用出来そうです。出力対象領域の面積のボーダーラインは上表掲載の20km2以上でちょうど良い感じなので、以下SQLコマンドでフラグ（ActiveFlag）を更新。
[SQL]
UPDATE t_PrefectureBorderArea
SET ActiveFlag = 1
WHERE Area > 20000000 — 単位：平米
[/SQL]
ただし、これだと前記した対馬の真ん中の領域は20km2未満で除外されてしまいます。
この他、数件は手作業でフラグ調整しました。
・対馬、小豆島、倉橋島、伯方島の面積20km2未満の領域→フラグON
・出力地図のスペースの都合上、以下の領域はフラグOFFに
・東京都：八丈島・小笠原諸島
・鹿児島県：トカラ列島・奄美諸島（屋久島よりも南の諸島）
・沖縄県：沖縄本島以外
以上の作業で、出力対象の領域数はかなり絞れましたが、レコード数は約260万件と期待ほど減りませんでした。まだまだ多すぎるので、データ間引き作業が必要です。それは次回で・・・。
前回、国土交通省のページより市区町村境のデータを取得し、データベース（MySQL）への収納、GoogleMap 市区町村領域のポリゴン表示（GoogleMapsAPI）を行いました。今回はこの市区町村境の点情報のうち、都道府県境に該当する点のみをピックアップしてみます。結果からいうと、予想に反して結構手間がかかったのですが・・・。
まずは基本考察。前回取得した「市区町村境界」の点には、次の３タイプがあります。
(1) 都道府県境界
(2) 同一都道府県内の市区町村境界
(3) 他の市区町村と接していない点＝海岸線
都道府県の輪郭を取得するためには、上記のうち(2)の点を除外すれば良いことになります。
※厳密には「(2)の他に(1)(3)の条件も兼ね備えている」点は除外対象とはなりませんが・・・
前回作成した市区町村境の点情報テーブル（t_Point）はレコード数が1,000万件近く、このままでは重すぎるので、今回の作業に当たってテーブルを以下のコード（PHP）で都道府県毎に分割しました。
<?php
ini_set( 'error_reporting', ~E_WARNING );
ini_set( 'memory_limit', '1024M' );
set_time_limit(0);
$errMsg;
$mysqli = connectDBi($errMsg);
if ($mysqli) {
echo "MySQL接続OK<br>";
} else {
echo "MySQL接続NG<br>".$errMsg;
exit();
}
echo str_pad(" ",4096)."<br>";
for ($PrefCode = 1; $PrefCode <= 47; $PrefCode++) {
$strSQL = "CREATE TABLE IF NOT EXISTS t_Point_Pref" . sprintf("%02d",$PrefCode) . " (";
$strSQL .= " PrefectureCode smallint(2) unsigned zerofill NOT NULL,";
$strSQL .= " AdministrativeAreaCode int(5) unsigned zerofill NOT NULL DEFAULT '00000',";
$strSQL .= " SurfaceCode int(5) NOT NULL,";
$strSQL .= " CurveCode int(5) NOT NULL,";
$strSQL .= " PointNo int(5) NOT NULL DEFAULT '0',";
$strSQL .= " Lat decimal(15,8) DEFAULT NULL,";
$strSQL .= " Lng decimal(15,8) DEFAULT NULL,";
$strSQL .= " CityBorderFlag tinyint(2) NOT NULL DEFAULT '0',";
$strSQL .= " PRIMARY KEY (PrefectureCode,AdministrativeAreaCode,SurfaceCode,CurveCode,PointNo),";
$strSQL .= " UNIQUE (SurfaceCode,CurveCode,PointNo),";
$strSQL .= " KEY LatLng (Lat,Lng)";
$strSQL .= ") ENGINE=MyISAM DEFAULT CHARSET=utf8";
if (!$mysqli->query($strSQL)){
echo "失敗！". $strSQL ."<br>";
exit();
}
echo "Create Table : t_Point_Pref".sprintf("%02d",$PrefCode)."<br>";
ob_flush();
flush();
$strSQL = "INSERT INTO t_Point_Pref".sprintf("%02d",$PrefCode)." (PrefectureCode,AdministrativeAreaCode,SurfaceCode,CurveCode,PointNo,Lat,Lng)";
$strSQL .= " SELECT PrefectureCode,AdministrativeAreaCode,SurfaceCode,CurveCode,PointNo,Lat,Lng";
$strSQL .= " FROM t_Point WHERE PrefectureCode = ". $PrefCode;
if (!$mysqli->query($strSQL)){
echo "失敗！". $strSQL ."<br>";
exit();
}
echo "INSERT DATA TO : t_Point_Pref".sprintf("%02d",$PrefCode)."<br>";
ob_flush();
flush();
}
if ($mysqli) {
$mysqli->close();
}
echo "終了しました。<br>";
ini_restore('error_reporting');
ini_restore('memory_limit');
//MySQLへ接続
function connectDBi(&$err) {
$MySQL_SERVER = "サーバー名";
$MySQL_USER = "ユーザー名";
$MySQL_PASSWORD = "パスワード";
$MySQL_DBNAME = "データベース名";
$err = "";
$mysqli = new mysqli($MySQL_SERVER, $MySQL_USER, $MySQL_PASSWORD, $MySQL_DBNAME);
if ($mysqli->connect_errno) {
$err = "データベース接続に失敗しました。";
}
if (strlen($err) > 0) {
return false;
} else {
return $mysqli;
}
}
?>
時間のかかる処理なので、途中でブラウザがタイムアウトになってしまわないよう、途中経過を出力しob_flush()関数で出力バッファをフラッシュしています。
新規テーブル「t_point_pref_XX」（XXには都道府県コードが入る）には元テーブルt_pointにはないフィールド「CityBorderFlag」を追加しました（上記ソース24行目）。同一都道府県内の市区町村境界の点の場合に値=1となるようにして、都道府県境界データピックアップ時の除外対象とします。また、この市区町村境界点のフラグ更新は同じ緯度・経度のレコードが他にあるかどうかで判断することになるので、緯度・経度のインデックスも追加しましたが（27行目）、効果の程は比較していないので不明。
次はこのCityBorderFlagの更新処理。同一都道府県内の市区町村境界であるか否かのチェックです。
for ($PrefCode = 1; $PrefCode <= 47; $PrefCode++) {
$strSQL = "UPDATE t_Point_Pref".sprintf("%02d",$PrefCode);
$strSQL .= " INNER JOIN t_Point_Pref".sprintf("%02d",$PrefCode)." AS t_Point_LatLng";
$strSQL .= " ON t_Point_Pref".sprintf("%02d",$PrefCode).".Lat = t_Point_LatLng.Lat";
$strSQL .= " AND t_Point_Pref".sprintf("%02d",$PrefCode).".Lng = t_Point_LatLng.Lng";
$strSQL .= " AND (t_Point_Pref".sprintf("%02d",$PrefCode).".SurfaceCode <> t_Point_LatLng.SurfaceCode";
$strSQL .= " OR t_Point_Pref".sprintf("%02d",$PrefCode).".CurveCode <> t_Point_LatLng.CurveCode)";
$strSQL .= " SET t_Point_Pref".sprintf("%02d",$PrefCode).".CityBorderFlag = 1";
if (!$mysqli->query($strSQL)){
echo "失敗！". $strSQL ."<br>";
exit();
}
echo "UPDATE CityBorderFlag : t_Point_Pref".sprintf("%02d",$PrefCode)."<br>";
ob_flush();
flush();
}
そして、都道府県の輪郭形成の点を収納するテーブルを用意
・t_PrefectureBorder
[SQL]
CREATE TABLE IF NOT EXISTS `t_PrefectureBorder` (
`PrefectureCode` smallint(2) unsigned zerofill NOT NULL,
`BorderCode` int(5) NOT NULL,
`PointNo` int(5) NOT NULL DEFAULT ‘0’,
`Lat` decimal(15,8) DEFAULT NULL,
`Lng` decimal(15,8) DEFAULT NULL,
PRIMARY KEY (`PrefectureCode`,`BorderCode`,`PointNo`),
KEY `LatLng` (`Lat`,`Lng`)
) ENGINE=MyISAM DEFAULT CHARSET=utf8;
[/SQL]
さて、ここからが本番。
(1) 都道府県境界
(2) 同一都道府県内の市区町村境界
(3) 他の市区町村と接していない点＝海岸線
前述のとおり、全データから(2)のデータを除外して、残りの点を繋げば都道府県の輪郭が出来上がるのですが、この繋ぎ方が問題です。正しい順序で繋いでいかないと、滅茶苦茶な形になってしまうので・・・。
そこで、以下のようロジックにて点を正しく繋がないかと考えました。
1. 上記の(2)に該当しない（CityBorderFlag=0）任意の点を始点として設定
2. 始点と同じ領域に属する点を順に繋いでいく。
3. CityBorderFlag=1となる点まで来たら、その点と同じ緯度・経度を持つ点を探索する。
その点が隣接市区町村との「接続点」になるので、今度はそこから2と同様に順に繋いでいく。
4. ある時点で始点に戻ってくるので、そこまでを全て繋げば都道府県の輪郭となる（島しょを除く）。
※島しょの場合は、大抵ひとつの市区町村のみに属しているので、上記1～4のうち3を経ずに1周して終了となる。
このロジックを実現すべく試行錯誤の結果、以下のコードで上手くいきました（多分）。
<?php
ini_set( 'error_reporting', ~E_WARNING );
ini_set( 'memory_limit', '2048M' );
set_time_limit(0);
$errMsg;
$mysqli = connectDBi($errMsg);
if ($mysqli) {
echo "MySQL接続OK<br>";
} else {
echo "MySQL接続NG<br>".$errMsg;
exit();
}
echo str_pad(" ",4096)."<br>";
for ($PrefCode = 1; $PrefCode <= 47; $PrefCode++) {
$strSQL = "SELECT t_Point_Pref".sprintf("%02d",$PrefCode).".* FROM t_Point_Pref".sprintf("%02d",$PrefCode);
$strSQL .= " INNER JOIN";
$strSQL .= " (SELECT PrefectureCode, AdministrativeAreaCode, SurfaceCode, CurveCode, Min(CityBorderFlag), Max(CityBorderFlag)";
$strSQL .= " FROM t_Point_Pref".sprintf("%02d",$PrefCode);
$strSQL .= " GROUP BY PrefectureCode, AdministrativeAreaCode, SurfaceCode, CurveCode";
$strSQL .= " HAVING Min(CityBorderFlag)=0 AND Max(CityBorderFlag)=1) AS q_BorderCurve";
$strSQL .= " ON t_Point_Pref".sprintf("%02d",$PrefCode).".PrefectureCode = q_BorderCurve.PrefectureCode";
$strSQL .= " AND t_Point_Pref".sprintf("%02d",$PrefCode).".AdministrativeAreaCode = q_BorderCurve.AdministrativeAreaCode";
$strSQL .= " AND t_Point_Pref".sprintf("%02d",$PrefCode).".SurfaceCode = q_BorderCurve.SurfaceCode";
$strSQL .= " AND t_Point_Pref".sprintf("%02d",$PrefCode).".CurveCode = q_BorderCurve.CurveCode";
$strSQL .= " ORDER BY SurfaceCode, CurveCode, PointNo";
echo "Pref=".$PrefCode."<br>";
echo " ----- Reading Point Data -----<br>";
ob_flush();
flush();
$rst = $mysqli->query($strSQL);
$rcount = $rst->num_rows;
$pointInfo = array();
$connectpoint = array();
$i = 0;
$j = 0;
$formerpoint = array();
if ($rcount > 0) {
while($col = $rst->fetch_array(MYSQLI_ASSOC)) {
if (isset($formerpoint["PrefectureCode"]) && $formerpoint["PrefectureCode"] == $PrefCode
&& ($formerpoint["SurfaceCode"] <> $col["SurfaceCode"] || $formerpoint["CurveCode"] <> $col["CurveCode"])) {
$i++;
$j = 0;
}
$col["iNo"] = $i;
$col["jNo"] = $j;
$col["CheckedFlagP"] = 0;
$col["CheckedFlagS"] = 0;
$pointInfo[$i][$j] = $col;
if ($j>0) {
if (!$pointInfo[$i][$j-1]["CityBorderFlag"] && $pointInfo[$i][$j]["CityBorderFlag"]) {
array_push($connectpoint,$pointInfo[$i][$j]);
} else if ($pointInfo[$i][$j-1]["CityBorderFlag"] && !$pointInfo[$i][$j]["CityBorderFlag"]) {
array_push($connectpoint,$pointInfo[$i][$j-1]);
}
}
$formerpoint = $col;
$j++;
}
}
$border = array();
$formerp = Null;
$p = get_BorderStartPoint($pointInfo);
$bNo = 0;
while (is_array($p)) {
if (!isset($border[$bNo])) {
$border[$bNo] = array();
}
if (!isset($formerp) || $p["Lat"] <> $formerp["Lat"] || $p["Lng"] <> $formerp["Lng"]) {
array_push($border[$bNo],$p); //前の点と同じ緯度経度の場合は登録しない
}
$pointInfo[$p["iNo"]][0]["CheckedFlagS"] = 1; //チェック（領域単位）
$pointInfo[$p["iNo"]][$p["jNo"]]["CheckedFlagP"] = 1; //チェック（ポイント単位）
$formerp = $p;
$p = get_NextPoint($pointInfo,$connectpoint,$p);
if (is_null($p)) {
$p = get_BorderStartPoint($pointInfo);
$bNo++;
}
}
for ($i=0;$i<count($border);$i++) {
//開始点と終了点の座標が同じになるよう調整
if (count($border[$i])>1) {
if ($border[$i][0]["Lat"] <> $border[$i][count($border[$i])-1]["Lat"]
|| $border[$i][0]["Lng"] <> $border[$i][count($border[$i])-1]["Lng"]) {
array_push($border[$i],$border[$i][0]);
}
}
echo "INSERT to t_PrefectureBorder: BorderCode=".$i." record count:".count($border[$i])."<br>";
ob_flush();
flush();
for ($j=0;$j<count($border[$i]);$j++) {
$strSQL = "INSERT INTO t_PrefectureBorder (PrefectureCode,BorderCode,PointNo,Lat,Lng)";
$strSQL .= " VALUES (".$PrefCode.",".$i.",".$j.",";
$strSQL .= $border[$i][$j]["Lat"].",".$border[$i][$j]["Lng"].")";
if (!$mysqli->query($strSQL)) {
echo "失敗！". $strSQL ."<br>";
exit();
}
}
}
//境界線が全て都道府県境となる領域（島など）
$strSQL = "SELECT PrefectureCode, AdministrativeAreaCode, SurfaceCode, CurveCode, Min(CityBorderFlag), Max(CityBorderFlag)";
$strSQL .= " FROM t_Point_Pref".sprintf("%02d",$PrefCode);
$strSQL .= " GROUP BY PrefectureCode, AdministrativeAreaCode, SurfaceCode, CurveCode";
$strSQL .= " HAVING Max(CityBorderFlag)=0";
$rst = $mysqli->query($strSQL);
$rcount = $rst->num_rows;
if ($rcount > 0) {
$bNo = count($border);
while($col = $rst->fetch_array(MYSQLI_ASSOC)) {
$strSQL = "INSERT INTO t_PrefectureBorder (PrefectureCode,BorderCode,PointNo,Lat,Lng)";
$strSQL .= " SELECT PrefectureCode,".$bNo.",PointNo,Lat,Lng";
$strSQL .= " FROM t_Point_Pref".sprintf("%02d",$PrefCode);
$strSQL .= " WHERE PrefectureCode = ". $col["PrefectureCode"];
$strSQL .= " AND AdministrativeAreaCode = ". $col["AdministrativeAreaCode"];
$strSQL .= " AND SurfaceCode = ". $col["SurfaceCode"];
$strSQL .= " AND CurveCode = ". $col["CurveCode"];
echo "INSERT to t_PrefectureBorder: BorderCode=".$bNo." (islands and others)<br>";
ob_flush();
flush();
if (!$mysqli->query($strSQL)){
echo "失敗！". $strSQL ."<br>";
exit();
}
$bNo++;
}
}
}
if ($mysqli) {
$mysqli->close();
}
echo "終了しました。<br>";
ini_restore('error_reporting');
ini_restore('memory_limit');
//輪郭の始点となる点を設定
function get_BorderStartPoint($pointInfo) {
echo " ----- Get Border Start Point -----<br>";
ob_flush();
flush();
for ($i=0; $i<count($pointInfo); $i++) {
if (!$pointInfo[$i][0]["CheckedFlagS"]) { //チェック済みの点のある領域がスタート地点になることはない
for ($j=0; $j<count($pointInfo[$i]); $j++) {
if (!$pointInfo[$i][$j]["CityBorderFlag"]) {
//出力
foreach ($pointInfo[$i][$j] as $key => $value) {
echo $value ." ";
}
echo "<br>";
ob_flush();
flush();
return $pointInfo[$i][$j];
}
}
}
}
echo " not found<br>";
ob_flush();
flush();
return NULL; //見つからない（＝全ての点をチェックした）場合
}
//都道府県輪郭線の次の点を探す
function get_NextPoint($pointInfo, $connectpoint, $p) {
$i = $p["iNo"];
$j = $p["jNo"] + 1;
if ($j >= count($pointInfo[$i])) {
$j = 0;
}
if (!$pointInfo[$i][$j]["CheckedFlagP"]) {
if (!$pointInfo[$i][$j]["CityBorderFlag"]) {
return $pointInfo[$i][$j];
} else {
return get_ConnectCityBorder($connectpoint, $pointInfo[$i][$j]);
}
} else {
return NULL; //チェック済みの点に戻ってきた（=１周した）
}
}
//市区町村の変わり目（接続点）を取得
function get_ConnectCityBorder($connectpoint, $p) {
echo "----- Connect City Border -----<br>";
$cp = array();
for ($c=0; $c<count($connectpoint); $c++) {
if ($connectpoint[$c]["Lat"] == $p["Lat"] && $connectpoint[$c]["Lng"] == $p["Lng"]) {
array_push($cp,$connectpoint[$c]);
}
}
//３つ以上の市町村境である点を考慮し、接続元の市区町村の次に登録されている点を返す
if (count($cp)>0) {
$chk = 0;
$cc = 0;
for ($c=0; $c<count($cp); $c++) {
if ($cp[$c]["iNo"] == $p["iNo"]) {
$chk = 1;
} else {
if ($chk == 1) {
$cc = $c;
break;
}
}
}
//出力
foreach ($cp[$cc] as $key => $value) {
echo $value ." ";
}
echo "<br>";
ob_flush();
flush();
return $cp[$cc];
}
return NULL; //見つからなかった場合（通常はあり得ないはず）
}
//MySQLへ接続
function connectDBi(&$err) {
$MySQL_SERVER = "サーバー名";
$MySQL_USER = "ユーザー名";
$MySQL_PASSWORD = "パスワード";
$MySQL_DBNAME = "データベース名";
$err = "";
$mysqli = new mysqli($MySQL_SERVER, $MySQL_USER, $MySQL_PASSWORD, $MySQL_DBNAME);
if ($mysqli->connect_errno) {
$err = "データベース接続に失敗しました。";
} else {
//文字化け対策
if (!$mysqli->set_charset("utf8")) {
$err = "文字コードセットに失敗しました。";
}
}
if (strlen($err) > 0) {
return false;
} else {
return $mysqli;
}
}
?>
データベース上で操作するのは難しそうなので、まずは変数（配列）に読み込んでから作業することにしました。そのSQL文が17～27行目の部分。
※以下は東京都（都道府県コード=13）の場合
SELECT t_Point_Pref13.* FROM t_Point_Pref13
INNER JOIN
(SELECT PrefectureCode, AdministrativeAreaCode, SurfaceCode, CurveCode, Min(CityBorderFlag), Max(CityBorderFlag) FROM t_Point_Pref13 GROUP BY PrefectureCode, AdministrativeAreaCode, SurfaceCode, CurveCode HAVING Min(CityBorderFlag)=0 AND Max(CityBorderFlag)=1) AS q_BorderCurve
ON t_Point_Pref13.PrefectureCode = q_BorderCurve.PrefectureCode
AND t_Point_Pref13.AdministrativeAreaCode = q_BorderCurve.AdministrativeAreaCode
AND t_Point_Pref13.SurfaceCode = q_BorderCurve.SurfaceCode
AND t_Point_Pref13.CurveCode = q_BorderCurve.CurveCode
ORDER BY SurfaceCode, CurveCode, PointNo
パフォーマンス向上のため読み込むレコード数はなるべく少なくしたいので、サブクエリ（緑フォント部分）を使用して対象となる領域（市区町村）を絞りました。
それが、Min(CityBorderFlag)=0 AND Max(CityBorderFlag)=1 の部分。CityBorderFlagは同一都道府県内で他の市区町村と接している場合に1となるので、Min(CityBorderFlag)=1となる領域は同一都道府県の市区町村に完全に囲まれていることを意味し、都道府県の輪郭にはなり得ません。
一方、Max(CityBorderFlag)=0となるのは同一都道府県の市区町村と全く接していない領域で（主に島しょ。一部、都道府県またぎの飛び地のケースもあり）、都道府県の輪郭になりますが全点が輪郭形成対象（=選定作業不要）なので後回し。結果、上記の式を満たす残りの領域の点を配列$pointInfo[$i][$j]に読み込んでいます（53行目）。作業対象の点は領域毎に管理したいので2次元配列とし、同一領域の点を$pointInfo[$i]に収納しています。
この他「CityBorderFlag=1かつ前後どちらかの点のCityBorderFlag=0である」点は隣接市区町村との接続点なので、$connectpointに別途収納しています（55～59行目）。
※この変数がなくても接続点探査は出来るのですが、$pointInfoで探査するとデータの絶対数が多い北海道などではパフォーマンスが落ちることが判明し、後で追加しました。
対象点を配列に収納した後の流れは、
(1) 輪郭の始点となるポイントを探査（関数:get_BorderStartPoint 148～173行目）
(2) 同一領域の点を順に繋ぐ（関数:get_NextPoint 175～191行目）
(3) 隣接市区町村境界点（CityBorderFlag=1）に達したら、接続すべき次の領域の点を緯度・経度を元に探査する。（関数:get_ConnectCityBorder 193～226行目）
これらの処理で繋いだ点を順番に配列$borderに収納していき（75行目）、対象の点が見つからなくなった時点では輪郭を1周して始点に戻ってきているはずです。 これを2つ目の輪郭、3つ目の輪郭・・・と繰り返し全部の対象点をチェックしたら、それらをテーブルt_PrefectureBorderに収納しています（99～105行目）。
最後に、除外していた島しょ等のデータを丸ごとt_PrefectureBorderに収納して（109～136行目）、１都道府県についての作業は終了です。
これを47都道府県分繰り返した結果、t_PrefectureBorderのレコード数は4,763,506・・・予想より随分多くなってしまいました。一番データ数の多い都道府県は長崎県で、レコード数は633,523と全体の13%を占めています。面積はむしろ狭い方ですが島の数がとにかく多く、6,284もの領域に分かれているのがデータ数が多くなっている要因です。
前回同様、今回も取得データの検証のためポリゴン表示のテスト地図をGoogleMapsAPIにて作成しました（ＭＡＰ表示）。埼玉県や奈良県などデータ数が少ない所はさほど待たずに表示されますが（図-1）、島の多い東京都や広大な北海道の表示には時間がかかります（図-2）。さらにデータ数の多い長崎県や鹿児島県に至ってはメモリオーバーやタイムアウトで表示出来ません。
このままでは、都道府県別の白地図表示など到底出来ません。日本全体の地図表示では今より輪郭の精度が落ちても全然問題ないですし、島しょについても淡路島や佐渡島などはともかく小さい島々は表示する必要もないので、次回はデータ数削減に取り組みます。
 図-1 選択した都道府県をポリゴン表示（埼玉県） |
 図-2 選択した都道府県をポリゴン表示（北海道） |
6月14日（土）
六本木ヒルズの森アーツセンターギャラリーにて「こども展」を観覧してきました。
こども展 ちらし
本展覧会の核となっているのはベルト・モリゾの娘ジュリー。ジュリーを描いた作品の数々の他、後に画家となった本人が甥を描いた作品も出品されていました。ポスターやちらしの表紙にも使用されているルノアールの「ジュリー・マネの肖像、あるいは猫を抱く子ども」も良いのですが、やはりモリゾの作品が娘への愛情が感じられて一番でした。
ジュリー関連以外では、何といってもモーリス・ドニ！ 「子どもを描かせたらドニの右に出る者はいないのではないか」と個人的には思います。
この他、ピカソの作品も・・・。キュビズムで描かれた子ども達（クロード、パロマ）は、その絵を見て一体どのような感想を持ったのでしょうか？
教室にて、子どもたちの学習
＜アンリ・ジュール・ジャン・ジョフロワ＞
庭のウジェーヌ・マネとその娘
＜ベルト・モリゾ＞
トランペットを吹くアコ
＜モーリス・ドニ＞
本展は、単に「子ども」が描かれた作品を寄せ集めただけ？という気がしないでもありませんが、個々の作品は良いものが多く楽しめました。見たことのある作品が多く、真新しさは感じませんでしたが・・・。
6月29日（日）まで。
