国土数値地理情報から白地図を作る。2


前の記事では市区町村(行政区域)の境界線を全て描画していました。そのままでは使いにくいので市区町村の境界線は消去して県の輪郭だけを描画するようにします。データにはすべての市区町村の境界線が入っています。そこですべてのデータについて重なったデータを消去します。

1 境界線のデータを消す

市区町村のデータと海岸線(県境)のデータの区別は次の様な仕組みになります。

上図は北海道の松前町、福島町、知内町、上ノ国町付近の地図です。欲しいデータは青い部分、海外線の部分。消したいデータは赤い部分で町と町の境界線になります。市区町村のデータを比較していき同じ座標のデータあればそれは境界線ということで削除するデータになります。すべてのデータを比較してデータを消せばよいので次のような処理をすることになります。

$imax = count($coordinate[$i]);
for ($i=0 ;$i<$imax ; ++$i){
  for ($m = $i+1 ; $m<$imax ; ++$m){
    $jmax = count($coordinate[$i]);
    $nmax = count($coordinate[$m]);

    for ($j=0 ;$j<$jmax ; ++$j){//各配列のjとnを比較するために2重ループ
      for ($n=0 ;$n<$nmax ; ++$n){
        if ($coordinate[$i][$j] === $coordinate[$m][$n]){
          $coordinate[$i][$j] = "";
          $coordinate[$m][$n] = "";
        }//if
      }//for $n ;
    }//for $j
  }
}

上記のコードではすべての行政区域を順番に調べています。そして座標データをすべて順番に調べ重なった(同一)の物があれば空欄にしています。この後に空欄のデータを削除して、といった処理を行います。この方法でも良いのですが非常に時間がかかります。明らかに重なっていないデータについてもチェックをするので無駄な処理が多数あります。そこで配列をまとめて処理をする方法を使います。配列に重なった部分があれば(同一の座標データがあれば)削除する処理を行い、なければ何もしない。という処理になります。

//行政区域の重なりを調べ、重なっている場合にはかさなっている座標を削除する。
  $figMax = count($coordinate);
  for ($fig1=0 ;$fig1<$figMax ; $fig1++){
    print("  ($fig1/$figMax)\r");
    for ($fig2=$fig1+1 ;$fig2<$figMax ; $fig2++){
      $checkArray = array();
      $checkArray = array_intersect($coordinate[$fig1],$coordinate[$fig2]);//2つの配列で重なった部分をしらべる
      if (count($checkArray) !== 0){//重なりがある場合
        $doublePoint = array();//重なっている座標を入れる配列。
        $doublePoint = array_intersect($coordinate[$fig1],$coordinate[$fig2]);//重なった値を返す
        $doublePoint = array_values($doublePoint);//キーをつめる

        //$coordinate[$fig1]と$doublePointで共通する値を削除する。
        $coordinate[$fig1] = array_diff($coordinate[$fig1],$doublePoint);//配列の差を求める。$coordinate[$fig1]から共通の項目を引く

        //$fig2について
        $coordinate[$fig2] = array_diff($coordinate[$fig2],$doublePoint);

        $coordinate[$fig1] = array_values($coordinate[$fig1]);//キーをつめる
        $coordinate[$fig2] = array_values($coordinate[$fig2]);//キーをつめる
      }
    }
  }

$fig1と$fig2が比較をする行政区域(または島)の座標データが入っている配列の番号になります。この$fig1と$fig2を二重のループで回します。その中で
array_intersect
で同一のデータを配列$doublePointにいれています。この配列が空なら重なった部分が無い。配列に中身があるならその部分は重なった座標なので、元の配列 $coordinate[$fig1]と $coordinate[$fig2]からそのデータを
array_diff
を使って削除します。
array_intersectとarray_diffについては以下のページに詳細があります。
array array_intersect ( array $array1 , array $array2 [, array $... ] )
array array_diff ( array $array1 , array $array2 [, array $... ] )

配列を削除した後には空の部分ができています。それをarray_valuesを使い隙間を埋めておきます。
array array_values ( array $array )

2. 余計な線を削除する

このままで描画すると市区町村の境界線の残りが直線として残っていしまいます。

完全に重ならない部分、線の始点、終点部分、そこいらへんが残ってしまいます。とにかく邪魔なものは消してしまいたいという事で座標データの個数が少ないデータは削除してしまいます。

  for ($fig1=0 ;$fig1<$figMax ; $fig1++){
    $dataCount = count($coordinate[$fig1]);//その行政区域の座標データ数を調べる
    if ($dataCount < 10){
      unset($coordinate[$fig1]);//座標データが10以下だったら削除
    }
  }
  $coordinate = array_values($coordinate);//キーをつめる
  $deleteCount = count($coordinate);

上のコードでは10個以下の点は削除しています。小さい島(岩)とかも同時に消えてしまいます。逆にある程度以下の島(岩)を消したい場合には

if ($dataCount < 400){
      unset($coordinate[$fig1]);
    }

とすると削除できます。これで余計な線が消えて次の様なものになります。

まだ残っている余計な線があります。これは市区町村の座標データの書きはじめが海岸線の区切りの良いところであれば良いのですが、そうなっていません。そこを切って直線でつなぐ部分がでてきるので上図のようになってしまいます。この線を消すために座標データを見て距離が離れている場合にはそのデータを削除するという処理を加えます。

$m = 0;
  $imax = count($coordinate);
  for ($i=0 ;$i<$imax ; $i++){
    $jmax = count($coordinate[$i]);
    $n = 0;//新しい配列に入れるためにnを0にする。
    for ($j=0 ;$j<$jmax-1 ; $j++){
      $latilongArray = explode(" ",$coordinate[$i][$j]);//空白区切りでデータを取り出す
      $latitude1String = $latilongArray[0];//各要素を取り出す
      $longitude1String = $latilongArray[1];

      $latilongArray = explode(" ",$coordinate[$i][$j+1]);//空白区切りでデータを取り出す
      $latitude2String = $latilongArray[0];//各要素を取り出す
      $longitude2String = $latilongArray[1];

      $x1 = floatval($latitude1String);//文字列のデータをfloatに変換する。
      $y1 = floatval($longitude1String);
      $x2 = floatval($latitude2String);
      $y2 = floatval($longitude2String);
      $dx = $x2 - $x1;
      $dy = $y2 - $y1;
      $s = $dx * $dx  + $dy * $dy;//$jと$j+1の距離の二乗 この値で判定する。

      $temp[$m][$n] = $coordinate[$i][$j];//$jのデータを配列に入れる。tempは作業用の仮の配列
      $n++;//次のためにnをひとつ進める。

      if (($s > 1.0E-4) || ($n > 30000)){//$j+1が遠くにあるときの処理 もしくは一つの図形が大きくなりすぎた時に分割する。
        $m++;//新規の図形 新しい配列に入れる。
        $n = 0;//新しい配列の先頭から
        //$temp[$m][$n] = $coordinate[$i][$j+1];
      }
    }// for $j
    $m++;//元の配列が一つ終わったら次の配列にする。
  }// for $i

  $coordinate = $temp;//元の配列に戻す。

$j番目の座標点と$j+1番目の距離を求めて距離が離れている場合には違う図形として新規に配列に入れています。距離はX座標とY座標の差を二乗して足しています。距離にするには最後にルートにいれる必要がありますが、厳密な距離を求める必要がないので、そのまま処理をしています。緯度によって実際の距離も変わるのでそこいらをちゃんと計算すると面倒なので適当に数値を決めています。

 if (($s > 1.0E-4) || ($n > 30000)){//$j+1が遠くにあるときの処理 もしくは一つの図形が大きくなりすぎた時に分割する。
        $m++;//新規の図形 新しい配列に入れる。
        $n = 0;//新しい配列の先頭から
        //$temp[$m][$n] = $coordinate[$i][$j+1];
      }

if文の1.0E-4というのがその適当に決めた距離になる部分です。
同時に1つの図形で点が30,000個以上の場合にも新規の図形として登録し直します。これはIllustratorで1つの図形が30,000個以上のアンカーポイントがあると表示されなくなるので、その対策のためです。

3.プログラム全文

以上のことでプログラムの全文は次のようになります。
plotPrefData2.php

<?php

ini_set('memory_limit', '256M');

mb_language("Japanese");//文字コードの設定
mb_internal_encoding("UTF-8");

$loadXML = "./国土数値情報行政区域/N03-170101_01_GML/N03-17_01_170101.xml";
$saveSVG= "hokkaido2.svg";//SVGファイルのファイル名

loadXMLsaveSVGXML($loadXML,$saveSVG);

function loadXMLsaveSVGXML($loadXML,$saveSVG){
  print("XML file loading...");
  $dom = new DOMDocument('1.0', 'UTF-8');
  $dom->preserveWhiteSpace = false;
  $dom->formatOutput = true;
  $dom->load($loadXML);

  $root = $dom->getElementsByTagName("*")->item(0);//最初の要素の中身を取り出す

  $i = 0;
  foreach( $root->getElementsByTagNameNS("*","Curve") as $itemGML ){//これでgml:Curveのところが取れるらしい。
    $dataArray = explode("\n",$itemGML->nodeValue);//改行コードで配列にわける

    $j=0;
    foreach ($dataArray as $item) {//各要素ごとに繰り返す
      $item = trim($item);
      if($item != ""){//空白文字列の時には処理をしない。
        $coordinate[$i][$j]=$item;
        $j++;
      }
    }
    $i++;
  }

  print("\n");
  //行政区域の重なりを調べ、重なっている場合にはかさなっている座標を削除する。
  $figMax = count($coordinate);
  for ($fig1=0 ;$fig1<$figMax ; $fig1++){
    print("  ($fig1/$figMax)\r");
    for ($fig2=$fig1+1 ;$fig2<$figMax ; $fig2++){
      $checkArray = array();
      $checkArray = array_intersect($coordinate[$fig1],$coordinate[$fig2]);//2つの配列で重なった部分をしらべる
      if (count($checkArray) !== 0){//重なりがある場合
        $doublePoint = array();//重なっている座標を入れる配列。
        $doublePoint = array_intersect($coordinate[$fig1],$coordinate[$fig2]);//重なった値を返す
        $doublePoint = array_values($doublePoint);//キーをつめる

        //$coordinate[$fig1]と$doublePointで共通する値を削除する。
        $coordinate[$fig1] = array_diff($coordinate[$fig1],$doublePoint);//配列の差を求める。$coordinate[$fig1]から共通の項目を引く
        //$fig2について
        $coordinate[$fig2] = array_diff($coordinate[$fig2],$doublePoint);

        $coordinate[$fig1] = array_values($coordinate[$fig1]);//キーをつめる
        $coordinate[$fig2] = array_values($coordinate[$fig2]);//キーをつめる
      }
    }
  }


  //削除してなくなった行政区域を調べる
  for ($fig1=0 ;$fig1<$figMax ; $fig1++){
    $dataCount = count($coordinate[$fig1]);//その行政区域の座標データ数を調べる
    if ($dataCount < 10){
      unset($coordinate[$fig1]);//座標データが0だったらその配列を削除
    }
  }
  $coordinate = array_values($coordinate);//キーをつめる
  $deleteCount = count($coordinate);
  print("deleted($deleteCount/$figMax)\n");


  //======================距離を調べて点と点が離れている時には別の配列データにする。
  $m = 0;
  $figMax = count($coordinate);
  for ($i=0 ;$i<$figMax ; $i++){
    $jmax = count($coordinate[$i]);
    $n = 0;//新しい配列に入れるためにnを0にする。
    for ($j=0 ;$j<$jmax-1 ; $j++){
      $latilongArray = explode(" ",$coordinate[$i][$j]);//空白区切りでデータを取り出す
      $latitude1String = $latilongArray[0];//各要素を取り出す
      $longitude1String = $latilongArray[1];

      $latilongArray = explode(" ",$coordinate[$i][$j+1]);//空白区切りでデータを取り出す
      $latitude2String = $latilongArray[0];//各要素を取り出す
      $longitude2String = $latilongArray[1];

      $x1 = floatval($latitude1String);//文字列のデータをfloatに変換する。
      $y1 = floatval($longitude1String);
      $x2 = floatval($latitude2String);
      $y2 = floatval($longitude2String);
      $dx = $x2 - $x1;
      $dy = $y2 - $y1;
      $s = $dx * $dx  + $dy * $dy;//$jと$j+1の距離の二乗 この値で判定する。

      $temp[$m][$n] = $coordinate[$i][$j];//$jのデータを配列に入れる。tempは作業用の仮の配列
      $n++;//次のためにnをひとつ進める。

      if (($s > 1.0E-4) || ($n > 30000)){//$j+1が遠くにあるときの処理 もしくは一つの図形が大きくなりすぎた時に分割する。
        $m++;//新規の図形 新しい配列に入れる。
        $n = 0;//新しい配列の先頭から
        //$temp[$m][$n] = $coordinate[$i][$j+1];
      }
    }// for $j
    $m++;//元の配列が一つ終わったら次の配列にする。
  }// for $i
  $temp = array_values($temp);//キーをつめる
  $coordinate = $temp;//元の配列に戻す。

  //
  $figMax = count($coordinate);
  for ($fig1=0 ;$fig1<$figMax ; $fig1++){
    $dataCount = count($coordinate[$fig1]);//その行政区域の座標データ数を調べる
    if ($dataCount < 10){
      unset($coordinate[$fig1]);//座標データが0だったらその配列を削除
    }
  }
  $coordinate = array_values($coordinate);//キーをつめる
  $deleteCount = count($coordinate);
  print("deleted($deleteCount/$figMax)\n");



  //=========================SVGで書き出す
  //スクリーンサイズ
  $screenWidth = 4000;
  $screenHeight = 4000;

  $north_max = 46.0;//最北端:北緯45度33分28秒 東経148度45分14秒  46.0
  $south_max = 20.0;//最南端:北緯20度25分30.6585秒 東経136度4分11.1766秒  20.0
  $east_max = 154.0;//最東端:北緯24度17分12秒 東経153度58分50秒 154.0
  $west_max = 122.0;//  最西端:北緯24度26分58秒 東経122度56分01秒 122.0

  $y_max = log(tan(pi()/4.0 + deg2rad($north_max)/2.0));
  $y_min = log(tan(pi()/4.0 + deg2rad($south_max)/2.0));
  $x_max = deg2rad($east_max);//横方向
  $x_min= deg2rad($west_max);

  $cLo = ($x_max - $x_min)/2 + $x_min;//描画中央になる東経

  //描画するための縮尺を求める
  $dx = $x_max - $x_min;
  $dy = $y_max - $y_min;
  $mX = $screenWidth/$dx;
  $mY = $screenHeight/$dy;

  //描画するための縮尺で縦方向と横方向で大きい方を使う。
  if ($mX < $mY){
    $multiple = $mX;
  } else {
    $multiple = $mY;
  }

  //y軸の描画用の定数
  $yOffset = $y_max;

  //var_dump($coordinate);

  //==========================SVGファイルを書き出す
  $contents='<?xml version="1.0" standalone="no"?>';
  $contents .="\n";
  $contents .='<!DOCTYPE svg PUBLIC "-//W3C//DTD SVG 1.1//EN" "http://www.w3.org/Graphics/SVG/1.1/DTD/svg11.dtd">';
  $contents .="\n";
  $contents .='<svg width="40cm" height="40cm" viewBox="0 0 4000 4000" xmlns="http://www.w3.org/2000/svg" version="1.1">';
  $contents .="\n";
  $contents .='<desc>都道府県の境界線データ</desc>';
  $contents .="\n";
  $contents .='<polyline fill="none" stroke="black" stroke-width="1" points="0,0 4000,0 4000,4000,0,4000" />';//図形の書き出し
  //境界線データの書き出し polylineでデータで書き出す。

  $imax = count($coordinate);
  print("imax=$imax\n");
  for ($i=0 ;$i<$imax ; $i++){
    $jmax = count($coordinate[$i]);
    //print("jmax=$jmax\n");
    $svgData = "";
    for ($j= 0;$j<$jmax ; $j++){
        //メルカトル図法による計算。

        $latilongArray = explode(" ",$coordinate[$i][$j]);//空白で区切る
        //print("strring=" . $coordinate[$i][$j] . "\n");
        $latitude1String = $latilongArray[0];//各要素を取り出す
        $longitude1String = $latilongArray[1];

        $radian = deg2rad(floatval($latitude1String));
        $x = deg2rad(floatval($longitude1String)) - $cLo;//横方向
        $y = log(tan(pi()/4.0 + $radian/2.0));//縦方向

        $x = $multiple * $x + $screenWidth/2;
        $y = $multiple * (1.0 - $y) - ($multiple * (1.0-$yOffset));
        $svgData .= $x . "," . $y . " ";

    }
    $contents .='<polyline fill="none" stroke="blue" stroke-width="1" points="';//図形の書き出し
    $contents .= $svgData . '" />' . "\n";
  }

  $contents .="</svg>";

  file_put_contents($saveSVG, $contents);

}

完成したSVGをIllustratorで開くと次の様に表示されます。(アウトライン表示)

ベクターデータになっているので県毎に色をつけたりしたくなりますがそれは出来ません。塗りを付けるためには閉じた曲線にする必要がありますが、下図の様に途切れた線になっておりこれをつながないと塗りを入れることはできません。

SVGでpolylineではなくpolygonにして閉じた曲線にすれば良いのですが、このデータ量でそれをやってしまうとIllustratorで開く事のできないデータになってしまいます。データを間引くなどの処理が必要になってきます。それはまた後ほど。

これで都道府県の輪郭をXMLから取り出す事はできました。次に各県で同様の処理をして日本地図を作成します。