これまで、カーネル密度推定IDW(逆距離加重補間)では、東京都の輪郭を覆うように長方形のエリアをとり、その範囲全てにメッシュ状に予測対象地点を配置していました(下図)。

予測対象地点メッシュ(東京都)

図1 予測対象地点メッシュ(東京都)


 予測結果を地図上に出力する際には、東京都の輪郭に合わせてクリップしますので、上図の右上や右下の部分の結果は実際には使用していません。
 今回もこれと同じ方法を採用すると次のようになります(下図の例では点の数は100×75=7,500)。
予測対象地点メッシュ(1都3県)

図2 予測対象地点メッシュ(1都3県)


 やはり地図出力には使用されない輪郭外の箇所が結構あります。それでもIDWまではパフォーマンスには大して影響しなかったのですが(計算よりも地図描画に要する時間が問題となっていました)、クリギングでは計算が格段と複雑なる上に、範囲を東京都→1都3県に広げたばっかりに、観測点の数も5倍以上になってしまったので、なるべく予測地点の数を減らしたい。かといってメッシュを粗くすると、地図も粗くなってしまうのでそれは避けたい・・・ということで、今回は輪郭外の点は予測地点から除外することにしました。

 ある点がポリゴンの中にあるか否かを判別出来れば予測地点をセレクト出来るが、すごく複雑な計算になってしまうのだろうか?・・・なんて思いつつ寝ネット探索してみると、結構単純なアルゴリズムであることが分かりました。
 対象となる点から、一方向(どの方向でも良い)へ直線を引いた時に、ポリゴンの輪郭線とクロスする回数で内部か外部が判別出来るというものです。
     クロスする回数 : 奇数・・・ポリゴン内部   偶数・・・ポリゴン外部

対象点がポリゴン内か否かを判別するアルゴリズム

図3 対象点がポリゴン内か否かを判別するアルゴリズム


 上図の例では、点Aはクロス回数が4回なのでポリゴン外、点Bは3回なのでポリゴン内、点Cは1回なのでポリゴン内、点Dは2回なのでポリゴン外と判別出来ます。点Eのように全くクロスしない場合は0回=偶数なのでポリゴン外となります。
 どんなに複雑な形のポリゴンであっても結局は線分の集まりに過ぎないので、1つ1つの線分とクロスするか否かをチェックすれば、そのクロス回数からポリゴン内か外かが判別出来ます。
 このアルゴリズムをコーディングして、以下のような関数を作成しました。

//
//  function innerPolygon
//     引数
//        xp   : 対象点のx座標
//        yp   : 対象点のy座標
//        poly : ポリゴンを構成する座標 [x,y]の配列
//
  function innerPolygon(xp,yp,poly) {
     var i,j,x1,y1,x2,y2,a,b;
     var cntcross = 0;
     for (i = 0; i <= poly.length - 1; i++) {
        if (i == poly.length - 1) {
           j = 0;
        } else {
           j = i + 1;
        }
        x1 = poly[i][0];
        y1 = poly[i][1];
        x2 = poly[j][0];
        y2 = poly[j][1];
        //対象点から水平方向(x→∞)に線を引いた時にポリゴンと交差する数をカウント
        //   交差数(cntcross)奇数→ポリゴン内
        //            偶数→ポリゴン外
        if ((xp == x1 && yp == y1) || (xp == x2 && yp == y2)
             || (xp == x1 && x1 == x2 && (yp > y1 != yp > y2)) || (yp == y1 && y1 == y2 && (xp > x1 != xp > x2))) {
           return true;  //ポリゴンを構成する点と座標が一致する場合、対象点がy=bまたはx=bの線分上にある場合→ポリゴン内(カウントする必要なし)
        } else if ((yp > y1 != yp > y2) && (xp < x1 || xp < x2)) {
           if (x1 == x2) {
              cntcross++;
           } else {
              a = (y2 - y1) / (x2 - x1); //傾き
              b = y1 - a * x1;           //Y切片
              if (((yp - b) / a) == xp) { //対象点が線分上にあればポリゴン内
                 return true;
              } else if (((yp - b) / a) > xp) { //直線上のy=yp時のxの値がxpより大きければ交差
                 cntcross++;
              }
           }
        }
     }
     if (cntcross % 2 == 1) {
        return true;
     } else {
        return false;
     }
  }

 参考にしたサイトのソース(C言語) 
 ・http://www.ecse.rpi.edu/Homepages/wrf/Research/Short_Notes/pnpoly.html
 は、非常にシンプルでたった10行なのですが、これだと対象点が輪郭線上にある場合の判別が点とポリゴンの位置関係(右か左か)によって異なるので、「輪郭線上の点は全てポリゴン内として判別する」というルールを追加した結果、上記のように若干長くなりました。地形のポリゴンでは、座標は大抵浮動小数点となり対象点が輪郭線と重なるケースは稀なので、ほとんど不要な処理ですが、まあそれでも十分シンプルなので・・・。
 なお、26,34行目の記述を「return true」から「return false」に変更すれば、輪郭線上はポリゴン外として判別されます。

 このアルゴリズムを実装した結果、図2の対象点は下図のように選別されました(赤い点の部分)。ただし、地図出力の際には「予測対象点を中心とした四角形」を描画するので、ギリギリでポリゴン外となる点を除外してしまうと一部に描画されない空白箇所が出来てしまうので、ポリゴン内部と判別された点と隣接する点も予測地点に含めることにしました(青い点の部分)。

予測対象地点メッシュ(1都3県)エリア内および隣接する点のみ選別

図4 予測対象地点メッシュ(1都3県) エリア内および隣接する点のみ選別


 図2と比較すると、対象点を約半分に減らすことが出来ました(7,500→3,766)。

 今回は、この選別した予測地点の情報をCSV保存することにしました。
 Javascript
 (1) 選別した対象点の情報をオブジェクトに収納しJSONに変換
 (2) Ajax(juery)でPHPへJSONデータをPOST送信
 (3) PHP側で受け取ったJSON文字列をデコードして連想配列にし、CSV保存
 というプロセスを採用しました。
 
 全ソースは以下のとおり

・MinamiKantoAirPollutionMakeCalpoint.html

<!DOCTYPE html>
<html>
<head>
<title>南関東(1都3県)の大気汚染濃度 クリギング予測地点</title>
<meta http-equiv="content-type" content="text/html; charset=UTF-8">
<meta name="viewport" content="initial-scale=1.0, user-scalable=no" />
<meta http-equiv="Content-Style-Type" content="text/css">
<script type="text/javascript" src="jquery.js"></script>
<script type="text/javascript" src="js/d3.js"></script>
<script type="text/javascript" src="js/queue.v1.js"></script>
<script type="text/javascript" src="js/MinamiKantoAirPollutionMakeCalpoint.js"></script>
<style>
TD input,textarea, select {
  font-size: 12px;
  ime-mode: disabled;
}

.estimatedcurve {
  fill: none;
  stroke-width: 2.0px;
}

.gridline {
  stroke: #cccccc;
  stroke-width: 1.0px;
  stroke-dasharray: 1,1;
  shape-rendering: crispEdges;
}

.axis path,
.axis line {
  fill: none;
  stroke: #333333;
  shape-rendering: crispEdges;
}
</style>
</head>
<body style="margin:0px;" onload="initialize()">
  <div style="background-color:#006666;padding:5px;">
  <table style="padding-left:5px;">
  <tr>
    <td style="font-size: 16px;font-weight: bold;color: #FFFFFF;">南関東(1都3県)の大気汚染濃度 クリギング予測地点</td>
  </tr>
  </table>
  </div>
</body>
</html>

・MinamiKantoAirPollutionMakeCalpoint.js

  var svg;
  var mapScale = 22000;  //マップスケール
  var mapWidth = 1000;   //マップ幅
  var mapHeight = 700;   //マップ高さ
  var mapTop = 20;
  var mapPath, mapBounds;
  var meshCountX=100,meshCountY=75;  //メッシュ数(横・縦)
  var cborder,pborder,outline;
  var jsonOutlineFileName = "json/OutlineMinamiKanto.json?=" + new Date().getTime();       //南関東(1都3県)全体の輪郭データ
  var jsonPrefFileName = "json/PrefectureBorderMinamiKanto.json?=" + new Date().getTime(); //南関東各都県の輪郭データ
  var jsonCityFileName = "json/CityBorderMinamiKanto.json?=" + new Date().getTime();       //市区町村境データ
  var calpointFileName = "csv/CalPoint.csv";                                               //出力ファイル(クリギング予測対象点情報)
  var cpmatrix = new Array();
  var cpselect = new Array();

  function initialize() {
     queue()
        .defer(d3.json, jsonCityFileName)
        .defer(d3.json, jsonPrefFileName)
        .defer(d3.json, jsonOutlineFileName)
        .await(ready);
  }

  function ready(error, data1, data2, data3) {
     cborder = data1;
     pborder = data2;
     outline = data3;
     setCalPoint();
     drawMap();
  }

  function setCalPoint() {

     var i,j;
     var outlineXY = new Array();

     mapBounds = d3.geo.bounds(outline);
     var mapcenter = [(mapBounds[0][0]+mapBounds[1][0])/2, (mapBounds[0][1]+mapBounds[1][1])/2];

     mapPath = d3.geo.path()
        .projection(d3.geo.mercator()
        .center(mapcenter)
        .translate([mapWidth/2, mapHeight/2])
        .scale(mapScale)
     );

     outline.features[0]["geometry"]["coordinates"][0].forEach(function(d,i) {
        outlineXY.push(mapPath.projection()([parseFloat(d[0]),parseFloat(d[1])]));
     });

     var xRange = d3.extent(outlineXY, function (d) { return d[0]; });
     var yRange = d3.extent(outlineXY, function (d) { return d[1]; });
     var gx = d3.range(xRange[0],xRange[1]+Math.pow(10,-10),(xRange[1]-xRange[0])/(meshCountX-1));
     var gy = d3.range(yRange[0],yRange[1]+Math.pow(10,-10),(yRange[1]-yRange[0])/(meshCountY-1));

     //地図範囲全体を覆う長方形上にメッシュ状に点を配置(cpmatrix)
     for (i = 0; i < gx.length; i++) {
        cpmatrix[i] = new Array();
        for (j = 0; j < gy.length; j++) {
           cpmatrix[i][j] = new Array();
           cpmatrix[i][j].x = gx[i];
           cpmatrix[i][j].y = gy[j];
           if (innerPolygon(cpmatrix[i][j].x,cpmatrix[i][j].y,outlineXY)) { //領域内か否かをチェック
              cpmatrix[i][j].inner = 1;
           } else {
              cpmatrix[i][j].inner = 0;
           }
        }
     }

     //cpmatrixのうち、1都3県の領域内の点及びそれに隣接する点のみをセレクト(cpmatrix)
     var p = 0;
     for (i = 0; i < gx.length; i++) {
        for (j = 0; j < gy.length; j++) {
           if (cpmatrix[i][j].inner == 0  
              && ((i > 0 && cpmatrix[i-1][j].inner == 1) || (i < gx.length - 1 && cpmatrix[i+1][j].inner == 1)
                 || (j > 0 && cpmatrix[i][j-1].inner == 1) || (j < gy.length - 1 && cpmatrix[i][j+1].inner == 1)
                 || (i > 0 && j > 0 && cpmatrix[i-1][j-1].inner == 1) || (i > 0 && j < gy.length - 1 && cpmatrix[i-1][j+1].inner == 1)
                 || (i < gx.length - 1 && j > 0 && cpmatrix[i+1][j-1].inner == 1) || (i < gx.length - 1 && j < gy.length - 1 && cpmatrix[i+1][j+1].inner == 1))) {
              cpmatrix[i][j].inner = 2; //隣接点
           }
           if (cpmatrix[i][j].inner > 0) {
              cpselect[p] = new Object();  //JSONで通信するためにオブジェクトにする
              cpselect[p].Lat = mapPath.projection().invert([cpmatrix[i][j].x,cpmatrix[i][j].y])[1];
              cpselect[p].Lng = mapPath.projection().invert([cpmatrix[i][j].x,cpmatrix[i][j].y])[0];
              cpselect[p].x = cpmatrix[i][j].x;
              cpselect[p].y = cpmatrix[i][j].y;
              p++;
           }
        }
     }

     saveCalPoint(cpselect);  //CSV保存

  }

  function drawMap() {

     svg = d3.select("body").append("svg")
        .attr("width", mapWidth)
        .attr("height", mapHeight)
        .attr("id", "svg")
        .style("position", "absolute")
        .style("top", mapTop+"px")
        .style("left", "0px")

     svg.selectAll(".outline")
        .data(outline.features)
        .enter()
        .append("path")
        .attr("d", mapPath)
        .attr("class", "outline")
        .style("fill","#c0ffc0")
        .style("opacity",0.3)

     svg.selectAll(".prefborder")
        .data(pborder.features)
        .enter()
        .append("path")
        .attr("d", mapPath)
        .attr("class", "prefborder")
        .style("fill", "none")
        .style("stroke", "#808080")
        .style("stroke-width", "1")

     svg.selectAll(".cityborder")
        .data(cborder.features)
        .enter()
        .append("path")
        .attr("d", mapPath)
        .style("fill", "none")
        .attr("class", "cityborder")
        .style("stroke", "#808080")
        .style("stroke-width", "0.3")

     var calpointline = svg.selectAll("g")
        .data(d3.map(cpmatrix).entries())
        .enter()
        .append("g")
        .attr("id",function(d,i) { 
           return "calpointline" + i;
        });

     calpointline.selectAll("calpoint")
        .data(function(d) {
            return d.value;
        })
        .enter()
        .append("circle")
        .attr("cx", function(d) {
           return d.x;
        })
        .attr("cy", function(d) {
           return d.y;
        })
        .attr("r", function(d) {
           if (d.inner > 0) {
              return 1.5;
           } else {
              return 0.8;
           }
        })
        .style("fill", function(d) {
           if (d.inner == 1) {
              return "#ff0000";
           } else if (d.inner == 2) {
              return "#0000ff";
           } else {
              return "#a0a0a0";
           }
        })

  }

//
//  function innerPolygon
//     引数
//        xp   : 対象点のx座標
//        yp   : 対象点のy座標
//        poly : ポリゴンを構成する座標 [x,y]の配列
//
  function innerPolygon(xp,yp,poly) {
     var i,j,x1,y1,x2,y2,a,b;
     var cntcross = 0;
     for (i = 0; i <= poly.length - 1; i++) {
        if (i == poly.length - 1) {
           j = 0;
        } else {
           j = i + 1;
        }
        x1 = poly[i][0];
        y1 = poly[i][1];
        x2 = poly[j][0];
        y2 = poly[j][1];
        //対象点から水平方向(x→∞)に線を引いた時にポリゴンと交差する数をカウント
        //   交差数(cntcross)奇数→ポリゴン内
        //            偶数→ポリゴン外
        if ((xp == x1 && yp == y1) || (xp == x2 && yp == y2)
             || (xp == x1 && x1 == x2 && (yp > y1 != yp > y2)) || (yp == y1 && y1 == y2 && (xp > x1 != xp > x2))) {
           return true;  //ポリゴンを構成する点と座標が一致する場合、対象点がy=bまたはx=bの線分上にある場合→ポリゴン内(カウントする必要なし)
        } else if ((yp > y1 != yp > y2) && (xp < x1 || xp < x2)) {
           if (x1 == x2) {
              cntcross++;
           } else {
              a = (y2 - y1) / (x2 - x1); //傾き
              b = y1 - a * x1;           //Y切片
              if (((yp - b) / a) == xp) { //対象点が線分上にあればポリゴン内
                 return true;
              } else if (((yp - b) / a) > xp) { //直線上のy=yp時のxの値がxpより大きければ交差
                 cntcross++;
              }
           }
        }
     }
     if (cntcross % 2 == 1) {
        return true;
     } else {
        return false;
     }
  }

//
//  予測地点情報をCSV保存
//
  function saveCalPoint(cp) { 

     var jsonObj = JSON.stringify(cp);
     var param = "filename=" + calpointFileName;
     param += "&json=" + jsonObj;

     jQuery.ajax({
        url : "savecalpoint.php?",
        type : "post",
        data: param,
        async : false,
        success: function(request){
           var res = request;
           if (res.length > 0){
              alert("エラー:"+res);
           } else {
              alert("出力しました");
           }
        },
        error: function() {
           alert('失敗しました');
        }
     });
  }

 ※29行目で呼び出しているdrawMap(97~173行目)は、確認用として地図上に対象点を出力するものなので(図4のMAPが出力されます)、CSV保存自体には不要です。

・savecalpoint.php

<?php
   $filename = $_REQUEST['filename'];
   $json = $_REQUEST['json'];
   $dataArray = json_decode($json,true);

   $cp_item = array('x','y','Lng','Lat','EstimatedValue');

   $file = fopen($filename, "w");
   fputcsv($file, $cp_item);

   foreach($dataArray as $data) {
      $cp_row = array();
      foreach ($cp_item as $item_name) {
         if (isset($data[$item_name])) {
            array_push($cp_row,$data[$item_name]);
         } else {
            array_push($cp_row,null);
         }
      }
      fputcsv($file, $cp_row);
   }
   fclose($file);
?>

 ※「EstimatedValue」は予測結果収納用のカラムで、本ソース実行時点では空です。

 輪郭データ(1都3県、都道府県、市区町村:JSON形式)はこちら(ZIP圧縮)

 この他、別途jquery.jsd3.jsqueue.jsをダウンロード、設置する必要があります。
 また、PHPのバージョンが古い場合(5.2未満)は、デフォルトでJSON拡張モジュールJSONが入っていないため追加が必要です。

 上記ソースの実行で出来上がった予測対象地点情報ファイル(CalPoint.csv)はこちら
 
 これでようやく準備が整いました。次回はいよいよクリギング計算・地図出力です。

コメント   

 前回のD3.jsで空間分析にチャレンジ(4) IDW(逆距離加重補間)から8か月以上。時間がかかってしまいましたが、関数最適化の手法の1つレーベンバーグ・マーカート法のアルゴリズムを実装出来たことで、ようやく空間統計の本命と思われるクリギング法に目途が立ちました。

 これまでは東京都のSPM濃度を対象にしてきましたが、いつも同じではつまらないので、今回は南関東の1都3県(東京都・埼玉県・千葉県・神奈川県)に範囲を広げ、対象データも二酸化窒素(NO2)の年平均値に変更しました(データ年度も2013年度から2014年度に更新)。対象測定局の位置は下図のとおり。

南関東(1都3県)大気汚染常時監視測定局

図1 南関東(1都3県)大気汚染常時監視測定局
このMAPを表示

  

 まずは観測値のセミバリオグラム(セミバリアランス)を求めます。観測値のばらつきを距離の関数として定義したものをバリオグラムと呼び、セミバリオグラムはバリオグラムを半分にしたもので、次のような式で定義されます。

セミバリオグラム定義式


 この計算をする前に、各観測点間の距離と分散(観測値:NO2濃度差の2乗)の関係をプロットすると以下のようになります(バリオグラムクラウド)。
バリオグラムクラウド

図2 バリオグラムクラウド


 点の数が多すぎて見づらいですが、上図のように各点間の距離はそれぞれ異なります。そのため、上のセミバリオグラム定義式では「2点間の距離がhとなる点(i,j)」とありますが、実際にはそんな組み合わせは通常1組しか存在せず、「集合」にはなり得ません。そこで、距離(h)を範囲指定することでグループ分けし、各グループについて定義式に基づいてセミバリオグラムを求めます。
 今回は、2点間の距離(最短:約1km、最長:約162km)を50のグループに分けました。その結果、上のバリオグラムクラウドは下図のように整理されました。
セミバリオグラム

図3 セミバリオグラム


 次は、このセミバリオグラムの近似曲線を求めます。近似曲線は下図のような形状となります。
セミバリオグラム 近似曲線の例

図4 セミバリオグラム 近似曲線の例


 上図のように、ある程度の距離に達すると、セミバリオグラムの値は一定となります。その値のことをシルと呼びます。また、シルに達する距離のことをレンジと呼びます(※シル・レンジの存在しないモデルもあります)。
 また、距離がゼロの時は、セミバリオグラムもゼロとなるのが通常ですが(同じ点なので)、測定誤差等を考慮してゼロでない値(>0)を設定する場合もあります(ナゲットと呼びます)。
 
 曲線のモデルには次のようなものがあります。

・円形モデル(spherical model)
      円形モデル

・指数モデル(exponential model)
      指数モデル

・ガウスモデル(gaussian model)
      ガウスモデル

・べき乗モデル(power model)
      べき乗モデル

・線形モデル(linear model)
      線形モデル

・ナゲット効果モデル(nugget effect)
      ナゲット効果モデル

    各モデル式のパラメータ

 ここで、いよいよ苦心の末に実装したレーベンバーグ・マーカート法モジュールの登場です。結果の返し方などをクリギング仕様へと若干変更していますが、基本的には以前作成したものと同じです。
 今回はナゲットは考慮せずに(=0)、上記それぞれのモデルについて、近似曲線を求めました(線形モデルとナゲット効果モデルは直線ですが)。この程度であれば、レーベンバーグ・マーカート法を使用するまでもなくガウス・ニュートン法で十分だったような気もしますが・・・。
 下図のような結果になりました。

各モデルのセミバリオグラム近似曲線

図5 各モデルのセミバリオグラム近似曲線
このグラフを表示

  
 円形モデルとガウスモデルが、残差平方和(residual sum of squares)が小さく、目で見た感じでも良くフィットしている印象です。
 次回は予測対象地点の設定をします。
 
 
 

 パナソニック汐留ミュージアムにて開催中の特別展「アール・ヌーヴォーのガラス展」を観覧してきました。
 都内で1時間少々の空き時間が出来たので、時間潰しをかねて・・・。


アール・ヌーヴォーのガラス展 ちらし

アール・ヌーヴォーのガラス展 ちらし

 ガラス工芸家といえばエミール・ガレとドーム兄弟くらいしか思いつかないような素人ですが、その2組の作品が展示の大半を占めているので、安心(?)して観覧出来ました。
 これだけ精巧な作品をどのようにして作成するのか、さっぱり分かりません。「複数の色のガラスを重ねている」等といった解説がありましたが、重ねるといってもそうそう出来るものではないでしょう・・・。
 照明の当て方によって見え方が変わる作品が何点かありましたが、実際にそれが分かるよう、時間経過とともに照明を変える演出も施されていました・・・さすがに家電メーカー傘下の美術館、こういった演出はお手のものですね。

 全体を通して「わあ、綺麗だなあ」と感心しながら眺めるだけでしたが、十分楽しめました。
 9月6日(日)まで。

コメント