地図について#3

ひつじかいさんの「ひつじかいの雑記帳」というブログを拝見し、その中に、地図の描写において色々と興味深い記事を拝読しました。これは自分でもやってみたい、そう思い、その記録を記載してみたいと思います。

Web上で操作可能な日本の白地図（都道府県別）を作る（3）　― 領域の面積を計算 ―

今回は、都道府県データの数が9.6百万件ほどあるので、こちらのデータを間引く作業です。方法としては、面積の小さい島々を省くということで、多角形の面積を求める方法の記載があります。お恥ずかしながら、多角形の面積を求める方法については知らず、勉強になります。

まずはSQL文にて、データ領域を作成します。

[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]

解説頂いたプログラムは次の通りです。

<?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 接続情報
    $MySQL_SERVER = "localhost";
    $MySQL_USER = "[user_name]";
    $MySQL_PASSWORD = "[password]";
    $MySQL_DBNAME = "map";
    $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;
    }
}

//
//  緯度・経度を平面直角座標に変換
//
//  国土地理院サイト（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;
}


//  座標原点とする日本の平面直角座標系を選択
//     平成十四年国土交通省告示第九号 より
//     ※北海道などは簡易的な選択を行っているため、本来の座標系とは異なる座標が選択される場合がある
//
//  $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
    );
}

ここでひつじかいさんが計算されたときと今回の平成31年データで計算した結果は次の通りでした。誤差変化として、向上しているものを〇、そうでないものをー、ありませんでしたが同値のものを＝としています。〇は24、－は23でした。

次に主な島についても比較されています。ここは結果のみ記載されていますので、今まで教示いただいた内容から、自分で計算及び確認する必要がありそうです。まずは20km²以上の面積である一覧を表示させます。なお東京都・小笠原諸島の母島はこの条件では漏れてしまうようなので、19km²より大きい一覧を次のSQLで抜き出してみます。

[SQL]
SELECT * FROM t_prefectureborderarea
WHERE 19 * 1000 * 1000 < Area AND Area < 4000 * 1000 * 1000
ORDER BY Area DESC;
[/SQL]

ここで得られた結果を眺めると、諸島の名前とBorderCodeとの紐づけを行う必要がありそうです。少し考えましたが、都道府県とそれに紐づけられたBorderCodeを表示するソフトを作成して、GoogleMapで確認した方が早そうですし、再度の確認も行えるかと考えました。ひつじかいさんの内容を参照しているにも関わらず、作成に半日ほど要してしまいましたが、次のようなソフトを用いて確認を行いました。

確認を行う上で、今まで聞いたことのない、読み方が分からない島を見ながら「まだまだ、知らないことがいっぱいあるのだなぁ」と感じながら「機会あったら行ってみたい」と思いを馳せ楽しんでおりました。なお対馬は、ひつじかいさんのご指摘の通りのこと、確認しました。計算及び確認結果は次の通りです。誤差変化として、向上しているものを〇、そうでないものをー、ありませんでしたが同値のものを＝としています。〇は35、－は35でした。

そして、この確認を行ったうえで次のSQLを実行させフラグを付与します。

[SQL]
UPDATE t_PrefectureBorderArea
SET ActiveFlag = 1
WHERE Area > 20000000 # 単位：平米
[/SQL]

なお次の点は手作業で行うそうです。

　・対馬、小豆島、倉橋島、伯方島の面積20km2未満の領域→フラグON
　・出力地図のスペースの都合上、以下の領域はフラグOFFに
　　・東京都：八丈島・小笠原諸島　　　　　　
　　・鹿児島県：トカラ列島・奄美諸島（屋久島よりも南の諸島）
　　・沖縄県：沖縄本島以外

ん～、ここまで、お一人で実施されたのかと思うと感服いたします。