地図について#4

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

Web上で操作可能な日本の白地図(都道府県別)を作る(4) ― レコード数削減(データ抽出) ―

前回に引き続き、データを間引く作業です。ひつじかいさんのブログには、データが約260万件との記載がありますが、今回試した平成31年データでは約960万件と3倍以上のようです。

データを間引く方法は、点と点との間隔を計算し、加えて点を緯度・経度でグループ化して処理していきます。まず次のSQL文を実行させフィールドを加えます。

ALTER TABLE t_PrefectureBorder
ADD ExtractFlag TINYINT(2) NOT NULL DEFAULT 0

そして解説頂いたPHPのソースコードは次の通りです。ひつじかいさんの記載されている内容と変わりありません。

<?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>";
    echo "  ----- Reading Data -----<br>";
    ob_flush();
    flush();

    $border = array();
    $strSQL = "SELECT t_PrefectureBorder.*";
    $strSQL .= " FROM t_PrefectureBorder";
    $strSQL .= " INNER JOIN t_PrefectureBorderArea";
    $strSQL .= " ON t_PrefectureBorder.PrefectureCode = t_PrefectureBorderArea.PrefectureCode";
    $strSQL .= " AND t_PrefectureBorder.BorderCode = t_PrefectureBorderArea.BorderCode";
    $strSQL .= " WHERE t_PrefectureBorderArea.PrefectureCode = " . $PrefCode;
    $strSQL .= " AND t_PrefectureBorderArea.ActiveFlag = 1";
    $strSQL .= " ORDER BY t_PrefectureBorder.BorderCode, t_PrefectureBorder.PointNo";
    $rst = $mysqli->query($strSQL);
    while ($col = $rst->fetch_array(MYSQLI_ASSOC)) {
        array_push($border, $col);
    }
    $rst->close();

    for ($b = 0; $b < count($border); $b++) {
        $border[$b]["PrefGrp"] = GroupPointLatLng($mysqli, $border[$b]);
        //ブラウザがタイムアウトしないよう1万件毎に出力
        if ($b % 10000 == 0) {
            echo "   PrefCode=" . $border[$b]["PrefectureCode"] . "   BorderCode=" . $border[$b]["BorderCode"] . "  PointNo=" . $border[$b]["PointNo"]
                . "  " . $b . "<br>";
            ob_flush();
            flush();
        }
    }

    $intervals = 1000;  //点と点の間隔
    echo "  ----- Active Point -----<br>";
    $formerborder = array();
    for ($b = 0; $b < count($border); $b++) {
        if (count($border[$b]["PrefGrp"]) > 0 && $border[$b]["PrefGrp"][0] == $PrefCode) {
            if (
                !isset($formerborder["BorderCode"]) || $border[$b]["BorderCode"] <> $formerborder["BorderCode"]
                || count($border[$b]["PrefGrp"]) > count($formerborder["PrefGrp"])
                || ($b + 1 < count($border) && count($border[$b]["PrefGrp"]) > count($border[$b + 1]["PrefGrp"]))
                || getDistance($formerborder, $border[$b]) > $intervals
            ) {
                $formerborder = $border[$b];

                $strSQL = "UPDATE t_PrefectureBorder";
                $strSQL .= " SET ExtractFlag = 1";
                $strSQL .= " WHERE Lat = " . $border[$b]["Lat"];
                $strSQL .= " AND Lng = " . $border[$b]["Lng"];
                if (!$mysqli->query($strSQL)) {
                    echo "失敗!" . $strSQL . "<br>";
                    exit();
                }
                //抽出された点の情報をブラウザ出力(タイムアウト対策)
                echo "   PrefCode=" . $border[$b]["PrefectureCode"] . "   BorderCode=" . $border[$b]["BorderCode"] . "  PointNo=" . $border[$b]["PointNo"]
                    . "  Grp=" . implode("-", $border[$b]["PrefGrp"]) . "<br>";
                ob_flush();
                flush();
            }
        }
    }
}

if ($mysqli) {
    $mysqli->close();
}

echo "終了しました。<br>";
ini_restore('error_reporting');
ini_restore('memory_limit');

//対象の点($point)と同じ緯度・経度の点を探索して配列で返す
function GroupPointLatLng($mysqli, $point)
{
    $strSQL = "SELECT DISTINCT PrefectureCode FROM t_PrefectureBorder";
    $strSQL .= " WHERE Lat = " . $point["Lat"] . " AND Lng = " . $point["Lng"];
    $strSQL .= " ORDER BY PrefectureCode";

    $rst = $mysqli->query($strSQL);
    $prefgrp = array();
    while ($col = $rst->fetch_array(MYSQLI_ASSOC)) {
        array_push($prefgrp, $col["PrefectureCode"]);
    }
    $rst->close();

    return $prefgrp;
}

//2点間の距離を計算(GRS80)
//  引数:$p0,$p1 --- 緯度・経度収納の配列   array("Lat"=>latitude,"Lng"=>longitude)
function getDistance($p0, $p1)
{

    $a = 6378137.0;             //長半径
    $f = 1.0 / 298.257222101;   //扁平率
    $b = $a * (1 - $f);           //短半径
    $e2 = (pow($a, 2) - pow($b, 2)) / pow($a, 2);    //第一離心率^2
    $dy = deg2rad($p0["Lat"] - $p1["Lat"]);          //2点の緯度差(ラジアン)
    $dx = deg2rad($p0["Lng"] - $p1["Lng"]);          //2点の経度差(ラジアン)
    $my = deg2rad(($p0["Lat"] + $p1["Lat"]) / 2.0);  //2点の緯度の中間
    $w = sqrt(1.0 - $e2 * pow(sin($my), 2));
    $m = ($a * (1 - $e2)) / pow($w, 3);              //子午線曲率半径
    $n = $a / $w;                                    //卯酉線曲率半径

    return sqrt(pow($dy * $m, 2) + pow($dx * $n * cos($my), 2));
}

//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;
    }
}

これで設定値である1000mごとのポイント、かつ隣接する県境に対する対応ができています。SQL文でデータ数を確認すると3万3千で、ひつじかいさんの計算結果とここで一致します。これを白地図で表示してみます。

白地図は描けており問題ありませんが、データ抽出において時間を要しているようなので、新たなテーブルを作成します。

[SQL]
CREATE TABLE t_JapanPrefectureBorder (
    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] 

そして次のPHP文を実行し、生成した[t_JapanPrefectureBorder]にデータをコピーしていきます。

<?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 "PrefCode=" . $PrefCode . "<br>";
    ob_flush();
    flush();

    $oldBorderCode = array();
    $strSQL = "SELECT DISTINCT BorderCode";
    $strSQL .= " FROM t_PrefectureBorder";
    $strSQL .= " WHERE t_PrefectureBorder.PrefectureCode = " . $PrefCode;
    $strSQL .= " AND t_PrefectureBorder.ExtractFlag = 1";
    $strSQL .= " ORDER BY t_PrefectureBorder.BorderCode";
    $rst = $mysqli->query($strSQL);
    while ($col = $rst->fetch_array(MYSQLI_ASSOC)) {
        array_push($oldBorderCode, $col["BorderCode"]);
    }
    $rst->close();

    for ($i = 0; $i < count($oldBorderCode); $i++) {
        echo "old bordercode=" . $oldBorderCode[$i] . "--->  new bordercode=" . $i . "<br>";
        $strSQL = "SELECT * FROM t_PrefectureBorder";
        $strSQL .= " WHERE PrefectureCode = " . $PrefCode;
        $strSQL .= " AND BorderCode = " . $oldBorderCode[$i];
        $strSQL .= " AND ExtractFlag = 1";
        $strSQL .= " ORDER BY BorderCode";
        $rst = $mysqli->query($strSQL);
        $point = array();
        while ($col = $rst->fetch_array(MYSQLI_ASSOC)) {
            array_push($point, $col);
        }
        $rst->close();

        for ($j = 0; $j < count($point); $j++) {
            $strSQL = "INSERT INTO t_JapanPrefectureBorder (PrefectureCode,BorderCode,PointNo,Lat,Lng)";
            $strSQL .= " VALUES (" . $PrefCode . "," . $i . "," . $j . "," . $point[$j]["Lat"] . "," . $point[$j]["Lng"] . ")";
            if (!$mysqli->query($strSQL)) {
                echo "失敗!" . $strSQL . "<br>";
                exit();
            }
        }
    }
}

if ($mysqli) {
    $mysqli->close();
}

echo "終了しました。<br>";
ini_restore('error_reporting');
ini_restore('memory_limit');

//対象の点($point)と同じ緯度・経度の点を探索して配列で返す
function GroupPointLatLng($mysqli, $point)
{

    $strSQL = "SELECT DISTINCT PrefectureCode FROM t_PrefectureBorder";
    $strSQL .= " WHERE Lat = " . $point["Lat"] . " AND Lng = " . $point["Lng"];
    $strSQL .= " ORDER BY PrefectureCode";

    $rst = $mysqli->query($strSQL);
    $prefgrp = array();
    while ($col = $rst->fetch_array(MYSQLI_ASSOC)) {
        array_push($prefgrp, $col["PrefectureCode"]);
    }
    $rst->close();

    return $prefgrp;
}

//2点間の距離を計算(GRS80)
//  引数:$p0,$p1 --- 緯度・経度収納の配列   array("Lat"=>latitude,"Lng"=>longitude)
function getDistance($p0, $p1)
{

    $a = 6378137.0;             //長半径
    $f = 1.0 / 298.257222101;   //扁平率
    $b = $a * (1 - $f);           //短半径
    $e2 = (pow($a, 2) - pow($b, 2)) / pow($a, 2);    //第一離心率^2
    $dy = deg2rad($p0["Lat"] - $p1["Lat"]);          //2点の緯度差(ラジアン)
    $dx = deg2rad($p0["Lng"] - $p1["Lng"]);          //2点の経度差(ラジアン)
    $my = deg2rad(($p0["Lat"] + $p1["Lat"]) / 2.0);  //2点の緯度の中間
    $w = sqrt(1.0 - $e2 * pow(sin($my), 2));
    $m = ($a * (1 - $e2)) / pow($w, 3);              //子午線曲率半径
    $n = $a / $w;                                    //卯酉線曲率半径

    return sqrt(pow($dy * $m, 2) + pow($dx * $n * cos($my), 2));
}

//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;
    }
}

ここで対象テーブルを[t_prefectureborder]から[t_JapanPrefectureBorder]に変更し、白地図を描いてみます。結果は前述と同じですが、途中の計算スピードは圧倒的にアップしています。

それにしても、隣接する県境においてグルーピングする方法を教示されていますが、私では考えられないこと、勉強になります。