逆ジオコーディングサーバ(リバースジオコーディング)を自作しよう
概要
緯度経度から住所を割り出す逆ジオコーディング(reverse geocoder)はYahooやGoogle等から各種APIサービスが提供されています。
参考
しかしながら、件数制限があったり、当然ながらインターネット接続が必要だったりするため、クローズド環境で動作する使いたい放題の逆ジオコーディングサーバを自作しました。
（日本国内の住所のみ対応です。）
- 2022/09/18更新
Githubに公開しました。https://github.com/wrwrh/reverse-geocoder-japan/tree/main
注意点
逆ジオロケーションサービスで単純に最近傍を探すことについての問題点は「Goで超高速かつスケーラブルな逆ジオコーディングAPIサーバーを作ってみた」で詳しく説明しています。
以下のコーディングは単純に街区レベルで最近傍を探しているため、この点をクリアできていません。
あくまで使用用途で要求される精度が出ているかを確認して、使うをことを前提としています。
（自身の用途では、多少間違ったいたところで手動修正すればヨシ！という前提）
住所データの入手
住所と緯度経度の紐づけデータには国土地理院の位置参照情報ダウンロードサービスを利用しました。
本来であれば、街区レベル(「○○町△丁目□番」)の情報のみで逆ジオコーディングしたいところですが、まだまだ住居表示がなされていない場所も多いので、大字・町丁目レベルの情報も併せてダウンロードします。
必要な都道府県のデータ(zipファイル)をすべてダウンロードして、CSVファイルを取り出します。
(街区レベルも大字・町丁目レベルもCSVのファイル名が同一のため別ディレクトリに展開)
今回ダウンロードした「東京」のデータの場合、ファイル名は「13_2021.csv」でした。
展開後はこんなディレクトリ構成にします。
└─csv
├─gaiku（街区レベル）
│ 13_2021.csv
│ 14_2021.csv
│ ...
└─oaza（大字・丁目レベル）
13_2021.csv
14_2021.csv
...
参考事項
今回の逆ジオコーダーのアルゴリズムは、検索対象の緯度経度から一番近い住所情報の代表点の緯度経度を探すものです。
そのため、飛び地や三日月の地域は、実際の住所情報情報とは異なる可能性がほかより高くなります。
※参考
本アルゴリズムでは逆ジオコーダーのみのため、なるべる街区レベルで最近傍の、住所情報を表示しています。
データベースの作成
ダウンロードしたCSVファイルからSQLiteのデータベースを作成します。
データベース仕様
- データベース名はgis.db
- 街区レベルと大字・町丁目レベルはテーブルを分ける(detail と simple)
- 街区レベルは、複数地点の代表点が同一の場合があるため、同一地点はまとめて1レコードにする
createDB.py
from pathlib import Path
import pandas as pd
import sqlite3
dir = Path(r'./csv/gaiku')#街区レベルのCSVファイル
sdir = Path(r'./csv/oaza')#大字レベルのCSVファイル
db_path = r'gis.db'
if __name__ == '__main__':
#まず大字・町丁目レベルの処理
sfiles = sdir.glob('*.csv')
sdf_concat = pd.DataFrame()
for file in sfiles:
with open(file) as f:
df = pd.read_csv(f,
encoding='shift_jis',
#必要な項目のみ読み取り
usecols=['都道府県名','市区町村名','大字町丁目名','緯度','経度'],
dtype={'都道府県名': str,'市区町村名': str,'大字町丁目名':str,'緯度':float,'経度':float}
)
df.fillna("",inplace=True)#nanだとgroupbyでエラーがおきるため置換
#複数都道府県を処理する場合、一つにまとめる
sdf_concat = pd.concat([sdf_concat,df])
sdf_concat = sdf_concat.rename(columns={
'緯度': 'Latitude',
'経度': 'Longitude',
'都道府県名': 'Prefecture',
'市区町村名': 'City',
'大字町丁目名': 'District',
})
#続いて、街区レベルを処理
files = dir.glob('*.csv')
df_concat = pd.DataFrame()
for file in files:
with open(file) as f:
df = pd.read_csv(f,
encoding='shift_jis',
usecols=['都道府県名','市区町村名','大字・丁目名','小字・通称名','街区符号・地番','緯度','経度'],
dtype={'都道府県名': str,'市区町村名': str,'大字・丁目名': str,'小字・通称名': str,'緯度':float,'経度':float, '街区符号・地番': str}
)
#sqlite用に並び替え
df = df.loc[:, ['都道府県名','市区町村名','大字・丁目名','小字・通称名','街区符号・地番','緯度','経度']]
df.fillna("",inplace=True)
#地番以外が同一のものをまとめる
grp = df.groupby(['都道府県名','市区町村名','大字・丁目名','小字・通称名','緯度','経度'])
print(grp.size())
#街区符号は-でつないで1レコードに
num = grp.agg({'街区符号・地番': lambda x: '-'.join(x)}).reset_index()
df_concat = pd.concat([df_concat,num])
#不要な文字を削る
df_concat = df_concat.replace('（大字なし）', '').replace('-.*-', '-', regex=True)
#クエリ用にカラム名変更
df_concat = df_concat.rename(columns={
'緯度': 'Latitude',
'経度': 'Longitude',
'都道府県名': 'Prefecture',
'市区町村名': 'City',
'大字・丁目名': 'District',
'小字・通称名': 'Street',
'街区符号・地番': 'Numbers'
})
with sqlite3.connect(db_path) as conn:
sdf_concat.to_sql('simple', con=conn,index=False)
#高速化のため緯度経度にindexを忘れずに
conn.execute('CREATE INDEX sidx ON simple(Latitude ASC,Longitude ASC)')
df_concat.to_sql('detail', con=conn,index=False)
conn.execute('CREATE INDEX didx ON detail(Latitude ASC,Longitude ASC)')
これでgis.dbが出来上がります。（全国分のデータで約400MB）
検索用PHPの作成
リクエストを受け付けるPHPファイルを作成します
仕様
- ファイル名はgis.php
- gis.dbとgis.phpは同一ディレクトリに置く
- 渡すパラメータはlat(緯度)とlng(経度)
- 緯度経度は十進数
- 結果はJSONで返す
- ウェブサーバにはCORS設定を忘れずに！
gis.php?lat=35&lng=136
gis.php
<?php
$lat = (float)filter_input(INPUT_GET,'lat',FILTER_SANITIZE_SPECIAL_CHARS);
$lng = (float)filter_input(INPUT_GET,'lng',FILTER_SANITIZE_SPECIAL_CHARS);
//日本以外除外
if($lat < 20 || $lat > 46 ||$lng < 122 || $lng > 154) return;
list($result,$max) = mainLoop($lat,$lng);
//出力のJSONを構築
if($result){
$max = ceil($max);
$arr = array("accuracy"=>$max,"geo"=>$result);
$json = json_encode($arr);
echo($json . "\n");
}
function mainLoop($lat,$lng){
$db = new SQLite3('gis.db',SQLITE3_OPEN_READONLY);
//街区レベル
$stmt = $db->prepare('SELECT Latitude,Longitude,Prefecture,City,District,Street,Numbers FROM detail WHERE (Latitude BETWEEN ? AND ?) AND (Longitude BETWEEN ? AND ?)');
//大字・町丁目レベル
$stmt2 = $db->prepare('SELECT Latitude,Longitude,Prefecture,City,District FROM simple WHERE (Latitude BETWEEN ? AND ?) AND (Longitude BETWEEN ? AND ?)');
//対象の緯度経度を中心として、東西と南北の四角形内に該当するデータを抽出する
$arr = array(500, 3000, 10000); //10km四方で日本国内はだいたい網羅可能（山奥ではヒットしないかも）
foreach($arr as $MaxLength){
//範囲内の住所情報を取得
$sqlResult = FetchAddr($lat, $lng, $stmt, $MaxLength);
$max = INF;
$minIndex = INF;
//ある程度の範囲内で街区情報が見つからなかった場合、大字・町丁目レベルも検索
if(count($sqlResult) == 0 && $MaxLength >= 3000){
$sqlResult = FetchAddr($lat, $lng, $stmt2, $MaxLength);
}
//見つかった住所情報から、検索地点に一番近い住所情報を取得
for($i = 0; $i < count($sqlResult); $i++) {
$tmpLength = calcLength($lat, $lng, $sqlResult[$i]->lat, $sqlResult[$i]->lng);
if ($max >= $tmpLength)
{
$max = $tmpLength;
$minIndex = $i;
}
}
//1ループで見つからなかった場合、さらに広い範囲で検索
if(is_infinite($minIndex)) {
continue;
}
else{
$db->close();
return [$sqlResult[$minIndex],$max];
}
}
$db->close();
return null;
}
//SQLの回答
class SQLResult
{
public $lat;
public $lng;
public $pref;
public $city;
public $district;
public $street;
public $numbers;
}
//２点間の緯度経度から距離を求める
function calcLength($startLat, $startLng, $endLat, $endLng)
{
//近似式だけど、遠近の比較に用いるのみなのでOK
$pi2 = 0.018;// Pi÷180
$Rx = 6371000.0;//地球平均半径
$startlatPi = $startLat*$pi2;
$endlatPi = $endLat*$pi2;
$d = $Rx*acos(cos($startlatPi)*cos($endlatPi)*cos($endLng*$pi2 - $startLng*$pi2) + sin($startlatPi)*sin($endlatPi));
return $d;
}
//ある緯度経度から、指定した距離、角度の緯度経度を取得する(ちょっとマジメにヒュベニの式を記述)
function calcLatLng($startLat, $startLng, $length, $angle)
{
//角度は北を0度とし、東は90度となる。
$Rx = 6378137; //赤道半径
$ecc2 = 0.006694; //第2離心率
$pi2 = 0.0174532; // Pi÷180
$pi3 = 57.29578; //180÷Pi
//ヒュベニの式で対象の緯度経度を算出
$tmpLat = $startLat * $pi2 + $length * cos($angle * $pi2) / $Rx * (1 - $ecc2) / pow(sqrt(1 - $ecc2 * pow(sin(($startLat * $pi2)), 2)), 3) / 2;
$w = sqrt(1 - $ecc2 * pow(sin($tmpLat), 2));
$endLat = $startLat + $length * cos($angle * $pi2) / $Rx * (1 - $ecc2) / pow($w, 3) * $pi3; //計算した緯度
$endLng = $startLng + $length * sin(($angle * $pi2)) / ($Rx / $w * cos($tmpLat)) * $pi3; //計算した経度
$latlng = array($endLat, $endLng);
return $latlng;
}
//最大と最小の緯度経度を算出、その範囲内のデータをSQLiteファイルから抽出
function FetchAddr($lat,$lng,$stmt,$length)
{
//$stmt->reset()//PHP7.2以前用
//検索地点を中心とした矩形内にある地点を抽出
$maxLat = calcLatLng($lat, $lng, $length, 0)[0];
$maxLng = calcLatLng($lat, $lng, $length, 90)[1];
$minLat = calcLatLng($lat, $lng, $length, 180)[0];
$minLng = calcLatLng($lat, $lng, $length, 270)[1];
$stmt->bindValue(1, $minLat ,SQLITE3_FLOAT);
$stmt->bindValue(2, $maxLat ,SQLITE3_FLOAT);
$stmt->bindValue(3, $minLng ,SQLITE3_FLOAT);
$stmt->bindValue(4, $maxLng ,SQLITE3_FLOAT);
$result = $stmt->execute();
$sqlResult = array();
$cnt = 0;
//クエリ結果を配列へ変換
while ($sdr = $result->fetchArray(SQLITE3_NUM)) {
$sqlResult[$cnt] = new SQLResult();
$sqlResult[$cnt]->lat = $sdr[0];
$sqlResult[$cnt]->lng = $sdr[1];
$sqlResult[$cnt]->pref = $sdr[2];
$sqlResult[$cnt]->city = $sdr[3];
$sqlResult[$cnt]->district = $sdr[4];
if(count($sdr)==7){
$sqlResult[$cnt]->street = $sdr[5];
$sqlResult[$cnt]->numbers = $sdr[6];
}
$cnt++;
}
return $sqlResult;
}
?>
実行結果
東京駅の緯度経度で検索
gis.php?lat=35.681363707720784&lng=139.7672604332142
レスポンス(accuracyは、検索地点と住所の代表地点との距離)
{"accuracy":14,"geo":{"lat":35.681252,"lng":139.767235,"pref":"\u6771\u4eac\u90fd","city":"\u5343\u4ee3\u7530\u533a","district":"\u4e38\u306e\u5185\u4e00\u4e01\u76ee","street":"","numbers":"9"}}
整形
{
"accuracy": 14
"geo": {
"lat": 35.681252,
"lng": 139.767235,
"pref": "東京都",
"city": "千代田区",
"district": "丸の内一丁目",
"street": "",
"numbers": "9"
}
}
後述
一応、PHPを使わずクライアントサイドのみで完結する方法はあった → WebWorker + indexedDB
(やってみたけど、初回の数百メガバイトのテキスト読込からDB変換のタイムラグがきつすぎる)
備忘録代わりに次の執筆ネタに・・。