Zenn
😀

逆ジオコーディングのためのデータ整備 #3 位置参照情報を取り込む

2020/10/27に公開
#
postgis
tech

はじめに

街区レベル位置参照情報をPostGISデータベースに格納します。ただし、これまでに格納している国土数値情報（行政区域）とリンクするようにしておきましょう。

リンクを考えなければ非常に楽なんですが、今回はかなり泥沼化しています。

今回は「汚い」

また、今回ご紹介するスクリプトは、とても汚いものですので、こういうかんじでやってますよ、これでもだいたいうまくいきますよ、という情報提供程度にとらえていただければと思います。

簡単でない点

街区レベル位置参照情報と国土数値情報（行政区域）とをつなげる際には、次のポイントに注意しなければなりません。

  • 街区レベル位置参照情報には市区町村コードがないので市区町村名でチェックしなければならない
  • そして表記が異なる場合があり、たとえば「ケ」「ヶ」の違いがある
  • 国土数値情報と街区レベル位置参照情報との作成タイミングが異なり、市町村合併について一方が反映され、もう一方が反映されていないことがある
  • そして合併後は大字・町丁目名が変わることが結構ある

必要なPerlモジュール

標準のPerlにないモジュール "DBD::Pg" が必要ですので、引っ張ってきておいて下さい。

位置参照情報を取得して抽出する

位置参照情報ダウンロードサービスからダウンロードして下さい。結構面倒くさいかも知れませんが、全国です。

zip/ に落とす

% mkdir zip

ここに思う存分落としてください。

次にcsvディレクトリを作り、zipからファイルを思う存分抽出して下さい。

% mkdir csv
% cd csv
% foreach F ( ../zip/*.zip )
  unziup $F
  end

今回は"CSV"と"csv"が混在している

csvファイルで配られていますが、今回はファイル名のサフィックスが"CSV"になっているものがあります。ファイル数が多いので見過ごしがちですが、気をつけてください。

% ls csv/*/*.CSV
...
% ls csv/*/*.csv
...

cshなら次のようなスクリプトでファイル名を変更しちゃいましょう。終わったら、トップディレクトリに戻っておきます。

% foreach F ( csv/*/*.CSV )
  set D = `echo $F|sed 's/\.CSV$/\.csv/'`
  echo "$F to $D"
  mv "$F" "$D"
end
% cd ..

文字コード変換、タブ区切り化を行う

文字コードをcp932からutf-8に変換し、タブ区切り化するスクリプトmktxt.perlを示します。

#!/usr/bin/perl

use utf8;
use Encode;


# 県名、市区町村名、字・町丁目名、街区符号、系番号、Ｘ、Ｙ、緯度、経度、住居表示フラグ、代表フラグ
# Ｘ、ＹはＸが南北、Ｙが東西
#

$line = 0;
while(<>) {
  chomp;
  $_ = decode('cp932',$_);
  @arr = split(/,/);
  @a = ();
  for($n = 0; $n < 11; $n++ ) {
    $arr[$n] =~ s/^\s*\"([^\"]*)\"\s*$/\1/;
    push(@a,$arr[$n]);
  }
  if( $line > 0 ) {
    $_ = join("\t",@a);
    $_ = encode('utf8',$_);
    print "$_\n";
  }
  $line++;
}

これを使って、いっきにタブ区切りファイルを作っちゃいましょう。ファイル名がcsvのままだけど気にしないことにします。

mkdir utf

foreach F ( `ls csv/*/*.csv | sort ` )
  set body = `echo $F | awk -F / '{print $NF;}'`
  set out = "utf/$body"
  echo "$F to $out"
  ./mktxt.perl $F > $out
end

g_mncplにない市区町村名を抽出する

g_mncplは、市区町村コードで管理していますが、位置参照情報には市区町村コードがありません。このため、市区町村名から市区町村コードを得る必要があります。

国土数値情報（行政区域）と位置参照情報とが完全にリンクしていればいいのですが、必ずしもそうではありません。リンクしきれない原因は、最初に言及していますが、表記ゆれと市町村合併です。

リンクしきれない市区町村については、手動でどうにかする必要があります。

まずは、g_mncplにない市区町村名を抽出します。

mkerrmcode.perlというファイル名のスクリプトを、次のように作ります。

#!/usr/bin/perl

use strict;
use utf8;

use Encode;

use DBI;

my $DB = "db";
my $C_PNAME = 0;
my $C_MNAME = 1;

binmode STDOUT, ":utf8";
binmode STDERR, ":utf8";

my %mcodehash = ();
my $conn = DBI->connect(
  "dbi:Pg:dbname=$DB", "", "",
  {pg_enable_utf8 => 1}
) || die "$!\n";

# "(県名) (郡名)(市区町村名)" の文字列を作成する。
# where条件
# xx000 - 県別未定地
# xx100, xx300 - 区別未定地 (百の位が奇数)
# 27140 - 堺市未定地(奇数じゃない)
my $sql = "SELECT mcode, pname || ' ' || (CASE WHEN gname IS NULL THEN '' ELSE gname END)||(CASE WHEN mname IS NULL THEN '' ELSE mname END) AS mname FROM g_mncpl LEFT JOIN g_pref ON g_mncpl.pcode=g_pref.pcode ".
  "WHERE mod(mcode,1000) <> 0 AND mod(mcode,1000) <> 100 AND mod(mcode,1000) <> 130 AND mod(mcode,1000) <> 150 AND mcode <> 27140;";

my $statement = $conn->prepare($sql);
$statement->execute();
while (my ($mcode, $mname) = $statement->fetchrow_array()) {
  $mcodehash{$mname} = $mcode;
}
print STDERR "All data are extracted.\n";

open(WE,">:utf8", "mcode-error.txt");
my @already = ();
while(<>) {
  chomp;
  # 入力は内部表現utf8でないので、デコードが必要
  # STDINに対してbinmode指定しているがダイヤモンド演算子は影響を受けない
  $_ = decode('utf8', $_);
  my @arr = split(/\t/);
  my $mname = $arr[$C_PNAME].' '.$arr[$C_MNAME];
  if( !grep(/^$mname$/,@already) ) {
    print STDERR " \"$mname\" ... ";
    my $mcode = $mcodehash{$mname};
    if( !$mcode ) {
      my $mname1 = $mname;
      $mname1 =~ s/ヶ/ケ/g;
      $mcode =$mcodehash{$mname1};
      if( !$mcode ) {
        print STDERR "Not found.\n";
        $mcode = '';
        print WE $arr[$C_PNAME]."\t".$arr[$C_MNAME]."\t\n";
      }
      else {
        print STDERR "\"ga\" problem.\n";
        print WE $arr[$C_PNAME]."\t".$arr[$C_MNAME]."\t$mcode\n";
      }
    }
    else {
      print STDERR "Found.\n";
    }
    push(@already, $mname);
  }
}
close(WE);

このスクリプトに(csvとなっていますが)タブ区切りファイルを突っ込みます。

% ./mkerrmcode.perl utf/*.csv
All data are extracted.
 "北海道 札幌市中央区" ... Found.
 "北海道 札幌市北区" ... Found.
...
 "岩手県 岩手郡滝沢村" ... Not found.
...
 "宮城県 宮城郡七ヶ浜町" ... "ga" problem.
...

市区町村コードと"(都道府県名) (市区町村名)"で整えた文字列との対応表を作成して、それにあわないものを出力するものです。

ここで、国土数値情報では「ケ」となっている箇所は位置参照情報では「ヶ」となっているので、配慮しています。

./mcode-error.txtが生成されます。

これは、市区町村名がdb.g_mncplに存在しなかったものの一覧です。

タブ区切りで左から県名、市区町村名、市区町村コードとなっています。「ケ」問題と推測される市区町村名については、推定される市区町村コードが付与されていますが、それ以外は市区町村コードが付与されていません。

青森県	上北郡六ヶ所村	2411
岩手県	岩手郡滝沢村	
宮城県	宮城郡七ヶ浜町	4404
茨城県	龍ヶ崎市	8208
埼玉県	南埼玉郡白岡町	
千葉県	鎌ヶ谷市	12224
千葉県	袖ヶ浦市	12229
千葉県	山武郡大網白里町	
神奈川県	茅ヶ崎市	14207

これに市区町村コード等を追加して対応表を上書きするためのデータを作成します。

市制施行や市区町村合併の際は、変更先市区町村コードを第3カラムに記入します。

市制施行や市区町村合併の際は、大字名に変更が発生することが多く有ります。たとえば○○町が合併して××市になった場合、「○○町 △△」が「××市 ○○町△△」となったり、「××市 ○○△△町」になったり「××市 △△町」とかになったりと、さまざまです。そこで、第4カラム以降に大字名変換規則を追加します。mksqls.perlで変換を行います。

変換式次の通りです。

  • '-○○': ○○を削除(正規表現可)
  • '○○$': $を字名に置換
  • '!○○!××': ○○を××に置換

大字名変換規則は、各市町村のウェブサイトに示されています。

今回は、市区町村合併はなく、埼玉県南埼玉郡白岡町および岩手県岩手郡滝沢村の市制施行がありました。
白岡市は大字名変更はないので、第3カラムに新市コードを記入するだけで、第4カラム以降は記入しません。
滝沢市の場合は第3カラムはいいとして、大字名については、対応不能なレベルで問題があるのが分かったので、とりあえず「滝沢字○○」を「○○」に変更し、「字××」を「××」に変更することとし、それ以上は考えないことにしました。

これで、次のような内容になります。

青森県	上北郡六ヶ所村	2411
岩手県	岩手郡滝沢村	3216	-滝沢字	-字
宮城県	宮城郡七ヶ浜町	4404
茨城県	龍ヶ崎市	8208
埼玉県	南埼玉郡白岡町	11246
千葉県	鎌ヶ谷市	12224
千葉県	袖ヶ浦市	12229
千葉県	山武郡大網白里町	12239
神奈川県	茅ヶ崎市	14207

SQLを作ってPostGISに格納させる

./mksqls.perlというファイル名のスクリプトを作ります。

なお、``$YEAR''はデータ作成年度にあわせて変更する必要があります。

また、EWKB作成にあたって、浮動小数点数にpackを使っているので、ビッグエンディアンの場合には、末尾付近で、エンディアンに関する記述変更する必要があります。

#!/usr/bin/perl

use utf8;
binmode STDIN,  ":utf8";
binmode STDOUT, ":utf8";
binmode STDERR, ":utf8";

use DBI;

my $DB = "db";
my $INDIR = "./utf/";
my $OUTDIR = "./sql/";

# READ directory
opendir(D, $INDIR) || die "$INDIR: $!\n";
my(@files) = ();
foreach my $f (readdir(D)) {
  if( $f =~ /\.csv$/ ) {
    push(@files,$f);
  }
}
@files = sort @files;

# g_mncpl table
my %mcodehash = ();
my $conn = DBI->connect("dbi:Pg:dbname=$DB","") || die "$!\n";
my $statement = $conn->prepare("SELECT mcode, pname || ' ' || (CASE WHEN gname IS NULL THEN '' ELSE gname END)||(CASE WHEN mname IS NULL THEN '' ELSE mname END) AS mname FROM g_mncpl LEFT JOIN g_pref ON g_mncpl.pcode=g_pref.pcode;");
$statement->execute();
while (my ($mcode, $mname) = $statement->fetchrow_array()) {
  $mcodehash{$mname} = $mcode;
}
print STDERR "DB data are extracted.\n";

# mcode-error.txt
my %err = ();
open(R,"<:utf8", "mcode-error.txt");
while(<R>) {
  chomp;
  my @a = split(/\t/);
  my ($pn,$mn,$mcode);
  $pn = shift(@a);
  $mn = shift(@a);
  $mcode = shift(@a) * 1;
  my $filter = [];
  foreach my $f ( @a ) {
    push(@$filter, $f);
  }
  $err{$pn.' '.$mn} = [$mcode,$filter];
}
close(R);
print STDERR "mcode-error.txt data are read.\n";

# mkdir
if( ! -d $OUTDIR ) {
  print STDERR "$OUTDIR not found. I will try to create.\n";
  mkdir($OUTDIR) || die "$OUTDIR: $!\n";
}

# created_at field value
my @todayarr = localtime();
my $today = sprintf('%04d-%02d-%02d', $todayarr[5]+1900, $todayarr[4]+1, $todayarr[3]);

@cnts = ();

foreach my $f (@files) {
  my $inpath = $INDIR.$f;
  my $outpath = $OUTDIR.$f;
  $outpath =~ s/\.csv$/.sql/;
  print STDERR "$inpath -> $outpath\n";
  open(R,"<:utf8",$inpath) || die "$inpath: $!";
  open(W,">:utf8",$outpath) || die "$outpath: $!";
  print W "COPY g_section (mcode, sernum, section, blockcode, syscode, lastat, geom, jpr_geom, rep) FROM stdin;\n";
  while(<R>) {
    chomp;
    my($pn,$mn,$section,$blockcode,$syscode,$north,$east,$lat,$lon,$f1,$rep,$f2,$f3) = split(/\t/);
    my $mname = $pn.' '.$mn;
    my $mcode;
    $mcode = $mcodehash{$mname};
    if( !$mcode ) {
      my $e = $err{$mname};
      $e || die "\"$mname\" not found !";
      my($mc,$filter) = @$e;
      $mcode = $mc;
      if( $filter ) {
        foreach my $f ( @$filter ) {
          if( !$f || $f eq '' ) {
          }
          elsif( substr($f,0,1) eq '!' ) {
            my ($ptn,$rep) = split(/\!/, substr($f,1));
            if( $section =~ /$ptn/ ) {
              $section =~ s/$ptn/$rep/;
              break;
            }
          }
          elsif( substr($f,0,1) eq '-' ) {
            # 削除
            my $f1 = substr($f,1);
            $section =~ s/$f1//;
            break;
          }
          else {
            # 追加
            my $sc = $f;
            $sc =~ s/\$/$section/;
            $section = $sc;
            break;
          }
        }
      }
    }
    $cnts[$mcode] = $cnts[$mcode] + 1;
    my $sernum = $cnts[$mcode];
    my $geod = &mkewkb(4612,$lon,$lat);
    my $jpr = &mkewkb(-1,$east,$north);
    my @outa = ($mcode, $sernum, $section, $blockcode, $syscode, $today, $geod, $jpr, $rep);
    print W join("\t",@outa)."\n";
  }
  print W "\\.\n";
  close(W);
  close(R);
}

sub mkewkb {
  my($s, $x,$y) = @_;
  my $hc = '';
  if( $s > 0 ) {
    my $tc = unpack('H8',pack('i',0x20000001));
    my $sc = unpack('H8',pack('i',$s));
    $hc = $tc.$sc;
  }
  else {
    $hc = unpack('H8',pack('i',0x00000001));
  }
  my $xc = unpack('H16',pack('d',$x));
  my $yc = unpack('H16',pack('d',$y));
  # ビッグエンディアンの場合は、'01'でなく'00'として下さい。
  my $wkb = '01'.$hc.$xc.$yc;
  return $wkb;
}

これを使って、sqlデータを作成します。

% ./mksqls.perl

これで、都道府県ごとに ./sql/??.sql というファイルが作られます。

テーブル作成

次の通りテーブルを作成します。

CREATE TABLE g_section (
  gid SERIAL PRIMARY KEY,
  mcode int REFERENCES g_mncpl(mcode),
  sernum int,
  section text,
  blockcode text,
  syscode int,
  lastat timestamp,
  geom GEOMETRY(POINT, 4612),
  jpr_geom GEOMETRY(POINT, 0),
  rep int,
  UNIQUE(mcode,sernum)
);
CREATE INDEX ix_g_section_syscode ON g_section (syscode);
CREATE INDEX ix_g_section_the_geom ON g_section USING GiST (geom);
CREATE INDEX ix_g_section_jpr_geom ON g_section USING GiST (jpr_geom);
CREATE INDEX ix_g_section_rep ON g_section (rep);
  • mcode - 市区町村コード
  • section - 町丁目名、大字名
  • blockcode - 街区符号、番地
  • syscode - 平面直角座標系の系番号
  • geom - JGD2000地理座標系による位置
  • jpr_geom - 平面直角座標系による位置

住居表示がされている場合の住所は、「××市 △△町×丁目 □番 ○号」というふうになっていると思いますが、「△△町×丁目」が町丁目名、「□番」が街区符号、に対応します。「○号」に対応する住居番号は、街区レベル位置参照情報にはありません。

また、repカラムは、街区レベル位置参照情報にある「代表点フラグ」ですが、逆ジオコーディングでは使用しません。

geomとjpr_geom

よく見ていただくと、JGD2000地理座標系(geom)と、平面直角座標系(jpr_geom)の、ふたつのジオメトリカラムにポイントデータを入れているのが分かると思います。

jpr_geomは平面直角座標系ですが、SRID指定をしていません(SRID=0)。

jpr_geomには複数の系が入り乱れるため、こうせざるを得ません。

逆ジオコーディングのクエリを渡す際には、計算対象位置のSRIDを0にあわせないとエラーとされます。SRIDを0にすること自体は、そういうEWKTを生成すればすみます(SRID=<数字>;を書かなければいいだけ)し、既に生成されたジオメトリに対してはST_SetSRIDを使えば問題ありません。ただし、空間参照系を間違っていないかは自分で管理しておかなければならないので、注意が必要です。

データ格納

% foreach F ( `ls sql/*.sql|sort` )
  echo $F
  psql db -f $F
end

おわりに

今回は街区レベル位置参照情報(ポイント)をPostGISに格納しました。これでとりあえず、都市計画区域内では街区符号までの住所を計算できるようになりました。都市計画区域外においても、少なくとも市区町村レベルの逆ジオコーディングは可能となりました。

今回は、PostGISの機能は全く使わず、主役はPerlでした。街区レベル位置参照情報には市区町村名しかないので、PostGISに入れていた市区町村テーブルを使って自動で変換し、かつ自動変換が不能な市区町村名については手動で訂正するためのスクリプト群を示しました（汚いスクリプトですが）。

ただ、滝沢市の大字名に問題があるのが発覚し、今回は回復不能であることが分かったので、まずい対応をせざるを得ませんでした。

次回は、とりあえずいくつかのクエリを投げてみようなかと思います。

本記事のライセンス

クリエイティブ・コモンズ・ライセンス
この記事は クリエイティブ・コモンズ 表示 4.0 国際 ライセンス の下に提供されています。

Discussion