地図について#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であるが、理科年表の式に当表記はない
    // ρ''=1として計算した所、正しいと思われる変換結果が得られた

    $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でした。

都道府県面積 計算値 計算値 
(今回) 
誤差 誤差 
(今回)
 
誤差変化
北海道83,457.4883,466.9083,437.450.01%-0.02%
青森県9,644.749,637.749,644.03-0.07%-0.01%
岩手県15,278.8915,276.3515,272.88-0.02%-0.04%
宮城県7,285.807,282.537,278.52-0.04%-0.10%
秋田県11,636.3211,638.5711,635.730.02%-0.01%
山形県9,323.469,326.069,327.710.03%0.05%
福島県13,782.7613,775.4613,782.02-0.05%-0.01%
茨城県6,095.846,096.006,096.440.00%0.01%
栃木県6,408.286,406.976,406.88-0.02%-0.02%
群馬県6,362.336,362.866,362.010.01%-0.01%
埼玉県3,798.083,796.293,797.15-0.05%-0.02%
千葉県5,156.625,157.795,157.810.02%0.02%
東京都2,188.672,193.982,193.180.24%0.21%
神奈川県2,416.052,422.272,415.820.26%-0.01%
新潟県12,583.8412,577.5212,579.07-0.05%-0.04%
富山県4,247.624,254.374,253.740.16%0.14%
石川県4,186.214,185.054,185.39-0.03%-0.02%
福井県4,189.894,188.764,189.80-0.03%0.00%
山梨県4,465.374,469.294,469.770.09%0.10%
長野県13,562.2313,552.7913,552.78-0.07%-0.07%
岐阜県10,621.1710,619.9010,619.89-0.01%-0.01%
静岡県7,780.607,773.137,770.81-0.10%-0.13%
愛知県5,165.165,171.065,171.870.11%0.13%
三重県5,777.355,774.975,773.63-0.04%-0.06%
滋賀県4,017.364,016.154,016.15-0.03%-0.03%
京都府4,613.264,604.764,611.65-0.18%-0.04%
大阪府1,901.421,911.551,905.010.53%0.19%
兵庫県8,396.478,397.248,399.820.01%0.04%
奈良県3,691.093,690.233,690.23-0.02%-0.02%
和歌山県4,726.324,724.194,723.98-0.05%-0.05%
鳥取県3,507.313,506.613,506.72-0.02%-0.02%
島根県6,707.986,707.486,707.43-0.01%-0.01%
岡山県7,113.247,109.897,114.09-0.05%0.01%
広島県8,479.818,475.048,477.92-0.06%-0.02%
山口県6,114.146,113.366,111.99-0.01%-0.04%
徳島県4,146.814,143.034,146.46-0.09%-0.01%
香川県1,876.581,874.851,875.83-0.09%-0.04%
愛媛県5,678.515,673.965,675.74-0.08%-0.05%
高知県7,105.207,101.637,102.52-0.05%-0.04%
福岡県4,979.424,982.404,985.790.06%0.13%
佐賀県2,439.672,441.992,440.630.10%0.04%
長崎県4,105.884,092.164,130.26-0.33%0.59%
熊本県7,404.897,407.857,408.340.04%0.05%
大分県6,339.826,337.886,339.64-0.03%0.00%
宮崎県7,736.087,734.397,733.78-0.02%-0.03%
鹿児島県9,188.999,182.869,185.82-0.07%-0.03%
沖縄県2,276.722,275.162,280.80-0.07%0.18%
全国377,961.73377,911.27377944.99-0.01%0.00%

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

島名(都道府県)面積 計算値 計算値 
(今回)
 
誤差 誤差 
(今回)
 
誤差変化
択捉島 (北海道)3,182.653,188.693172.480.19%-0.32%
国後島 (北海道)1,498.561,498.861489.670.02%-0.59%
佐渡島 (新潟県)854.53854854.62-0.06%0.01%
奄美大島 (鹿児島県)712.52711.36712.22-0.16%-0.04%
対馬 (長崎県)696.64444.09440.92-36.25%-36.71%
淡路島 (兵庫県)592.3590.36592.43-0.33%0.02%
天草下島 (熊本県)574.25574.13574.97-0.02%0.12%
屋久島 (鹿児島県)504.89504.13504.22-0.15%-0.13%
種子島 (鹿児島県)445.05444.72444.21-0.08%-0.19%
福江島 (長崎県)326.48325.74326.32-0.23%-0.05%
西表島 (沖縄県)289.3289.52289.570.08%0.09%
色丹島 (北海道)250.16251.09247.840.37%-0.93%
徳之島 (鹿児島県)247.77247.57247.82-0.08%0.02%
島後 (島根県)241.64241.4241.54-0.10%-0.04%
天草上島 (熊本県)225.38225.48225.930.04%0.24%
石垣島 (沖縄県)222.63221.98222.20-0.29%-0.20%
利尻島 (北海道)182.16181.99182.08-0.09%-0.04%
中通島 (長崎県)168.42167.87168.36-0.33%-0.03%
平戸島 (長崎県)163.68163.22163.37-0.28%-0.19%
宮古島 (沖縄県)159.26159158.97-0.16%-0.18%
小豆島 (香川県)153.35144.38153.25-5.85%-0.06%
奥尻島 (北海道)142.75142.6142.69-0.11%-0.04%
壱岐島 (長崎県)133.93133.85134.61-0.06%0.50%
屋代島 (山口県)128.43128.32128.46-0.09%0.02%
沖永良部島 (鹿児島県)93.6793.5793.65-0.10%-0.02%
江田島・能美島 (広島県)91.5491.1591.31-0.42%-0.25%
大島 (東京都)91.059190.71-0.05%-0.37%
長島 (鹿児島県)90.6390.5290.66-0.12%0.03%
礼文島 (北海道)80.9581.1181.260.20%0.38%
加計呂麻島 (鹿児島県)77.3977.1977.24-0.26%-0.20%
倉橋島 (広島県)69.5968.4168.54-1.69%-1.51%
八丈島 (東京都)69.4869.3669.10-0.18%-0.55%
下甑島 (鹿児島県)66.1265.9965.55-0.20%-0.87%
大三島 (愛媛県)64.5764.3564.57-0.33%0.00%
志発島 (北海道)59.559.6258.350.20%-1.94%
久米島 (沖縄県)59.1159.3559.530.40%0.71%
喜界島 (鹿児島県)56.9356.8756.75-0.11%-0.32%
西ノ島 (島根県)55.8655.755.76-0.28%-0.17%
三宅島 (東京都)55.4455.4255.19-0.04%-0.46%
能登島 (石川県)46.6946.4946.58-0.43%-0.25%
上甑島 (鹿児島県)44.1444.4444.190.68%0.12%
大島 (愛媛県)42.0741.7141.87-0.86%-0.47%
大崎上島 (広島県)38.3938.1938.26-0.52%-0.33%
久賀島 (長崎県)37.3537.2237.23-0.36%-0.31%
口永良部島 (鹿児島県)35.7735.835.810.09%0.12%
因島 (広島県)34.9834.8535.04-0.38%0.19%
中之島 (鹿児島県)34.4834.3934.41-0.25%-0.20%
針尾島 (長崎県)33.3333.2533.16-0.25%-0.52%
中ノ島 (島根県)32.3732.2332.29-0.43%-0.26%
生口島 (広島県)31.0631.0231.21-0.12%0.47%
若松島 (長崎県)3131.0931.130.28%0.44%
南大東島 (沖縄県)30.5730.5530.52-0.05%-0.17%
厳島(宮島) (広島県)30.3930.3530.33-0.14%-0.21%
大矢野島 (熊本県)29.9830.0630.200.28%0.72%
伊良部島 (沖縄県)29.129.1229.070.08%-0.11%
与那国島 (沖縄県)28.9128.8828.95-0.10%0.13%
諏訪之瀬島 (鹿児島県)27.6627.5927.61-0.24%-0.20%
宇久島 (長崎県)24.9324.8824.93-0.19%0.01%
奈留島 (長崎県)23.8223.6123.68-0.87%-0.59%
父島 (東京都)23.823.7823.44-0.09%-1.50%
新島 (東京都)23.2223.1923.44-0.14%0.96%
硫黄島 (東京都)23.1623.1523.73-0.06%2.47%
伊江島 (沖縄県)22.7722.7522.76-0.09%-0.06%
向島 (広島県)22.2322.1722.31-0.29%0.38%
中島 (愛媛県)21.1721.0921.27-0.39%0.49%
伯方島 (愛媛県)20.8820.1920.92-3.29%0.20%
伊平屋島 (沖縄県)20.5920.5920.650.00%0.30%
御蔵島 (東京都)20.5520.520.51-0.24%-0.20%
与論島 (鹿児島県)20.4720.520.560.13%0.43%
母島 (東京都)20.2120.1919.88-0.11%-1.63%

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

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

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

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

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

コメント