(Category:システム関連)

 前回までで、円の大きさが地球を1周するような大きなものにならない限りは等距離の円を同心円状に描くことが出来るようになりました。
 今回は、
 ・円の半径を地球を1周するような規模に設定した場合、本来とは逆の領域が選択されてしまう場合がある
 この不具合に対応したいと思います。

 ポリゴンを描く際に選択する領域をA、選択していない方の領域をBとした時、Google Maps APIではAとBの面積の差があまりないような場合には、A・Bのうち北側に位置する方の領域を選択する傾向にあるようです。そのため、南半球にある地点を円の中心点にする場合はAの方が南側になりますからBが選択されてしまいます。また、例えば東京から15000kmの等距離円などBの面積が明らかにAよりも小さくなる場合にはBが選択されます。

 google.maps.Polygonメソッドに「指定した緯度・経度がポリゴン内部に含まれるように領域選択する」オプションがあればいいのになあ・・・と思うのですが、残念ながら備わっていません。
 でも、そのオプションの有無を探している際に、google.maps.geometryライブラリのポリゴン関数にcontainsLocationというものがあるのを見つけました。

 ・Maps Javescript API ジオメトリ ライブラリより抜粋
 containsLocation(point:LatLng, polygon:Polygon) 
 この関数は、指定の地点がポリゴンの内部または境界に存在すれば、true を返します。
 
 これを利用すれば、次のようなロジックでポリゴン描画すれば、上記ケースにも対応出来るのではないかと考えました。

 1. 半径の大きい円のポリゴンをダミーで描く(visible=falseとして表示させない)
 2. 円の中心点がダミーポリゴンの内部に含まれているかどうかをcontainsLocation関数で調べる
 3. 内部に含まれている場合はそのままでOKなのでダミーと同じものをもう一度描き表示させる。
 4. 内部に含まれていない場合は逆の領域Bがポリゴンとして選択されてしまったことになる。
  領域Aを選択させるには、一旦地球(地図)の全領域を覆うポリゴンを描いた後、領域Bを中抜け
  領域としてくり抜いてしまえば良い。
 
 コードにすると次のようになります

   //地図全体を覆うポリゴン用座標データ(反時計回り)
   var polyWholeWorld = [
      new google.maps.LatLng(0, 180),
      new google.maps.LatLng(-90, 180),
      new google.maps.LatLng(-90, 90),
      new google.maps.LatLng(-90, 0),
      new google.maps.LatLng(-90, -90),
      new google.maps.LatLng(-90, -180),
      new google.maps.LatLng(0, -180),
      new google.maps.LatLng(90, -180),
      new google.maps.LatLng(90, -90),
      new google.maps.LatLng(90, 0),
      new google.maps.LatLng(90, 90),
      new google.maps.LatLng(90, 180),
      new google.maps.LatLng(0, 180)
   ];
   //ダミーのポリゴンを作成   
   var dummyCircle = new google.maps.Polygon({
      paths: polyCoords,  //ポリゴン領域の位置情報(緯度・経度)の配列(時計回り)
      visible: false
   });
   //中心点(myPoint)がダミーポリゴン内部に含まれているかどうかをチェック
   var innerFlag = google.maps.geometry.poly.containsLocation(myPoint, dummyCircle);

   if (innerFlag) {
      //内部にある場合は、そのままのデータでOK
      polyPath = [polyCoords];  
   } else {
      //内部にない場合は、地球全領域からダミーポリゴンの領域をくり抜く
      polyPath = [polyWholeWorld, polyCoords];  //polyWholeWorldとpolyCoordsは回りが逆
   }
   //等距離円の描画
   myCircle = new google.maps.Polygon({
      paths: polyPath,
      strokeWeight: 0,
      fillColor: red,
      fillOpacity: 0.4
   });
   myCircle.setMap(map);

 後は前回同様に内側の円領域を切り抜く処理を加える必要があります。通常のケースでは単に1つ内側の円領域データ(配列)を逆順にすれば良いのですが、「内側の円の段階で本来とは逆の領域が選択されてしまっている場合にはどうすればいいのだろうか?」という疑問が浮かびましたが、このケースでも正しい領域を描くために使用したデータ(地球全体を選択して領域Bをくり抜くパターン)をそのまま全て逆順にする次のような手法で上手くいきました。

 地球全体領域-外円の逆領域(領域B1)-{地球全体領域-内円の逆領域(領域B2)}

 この式を展開すれば 

 内円の逆領域(領域B1)-外円の逆領域(領域B1)

 これで良かったのか・・・でも考え方としては上式の方が分かりやすい気がするので、結果は同じですが今回は上式を採用しました。

計算・描画結果

計算・描画結果 ボツワナの首都(ハボロネ)から3000km,5000km,8000km,10000kmの同心円

 上手く出来たようです。しかし、ボツワナから10,000kmでは日本に全然届かない・・・やっぱり遠いなあ。

コメント   

 前回、一応それっぽい同心円を描いてみたものの、等距離円の距離が世界地図規模(千km単位)となると計算誤差が大きくなるため、今回はもう少し正確性の高い計算方法を探索します。

 以前に位置情報(緯度・経度)の分かっている2点間の距離を知りたいと思った時には、その求め方が記載されたサイトはWeb上に多数存在し、結構簡単に見つけられたのですが、今回必要としている「ある1点からの距離がrである点の緯度・経度を求める方法」はなかなか見つけられず、一番頼りになる国土地理院のサイトにも記載なし。「もうここはおとなしくgoogle.maps.Circle利用で妥協するか(同心円の色が重なってしまいますが)」と諦めかけたところで漸く見つけました・・・Vincenty法(順解法)。一旦見つけてしまえば「何で今まで見つけられなかったんだろう?」と思うのですが、探し物にはよくあることですね。
 始点の緯度・経度、始点からの距離・方位角を、このVincenty法(順解法)の式に代入することで到達点の緯度・経度が求まります。
 同心円を描くには始点の周り360度の点を全て求める必要があるので、以下のように関数化しました。

  //Vincenty法(順解法)
  function VincentyDirect(p1,deg1,s) {
/*
     	引数
        p1    : 中心となる点(点1)の緯度・経度(google.maps.LatLngクラス)
        deg1  : 求める点の方角(真北からの時計回りの角度。北:0、東:90、南:180、西:270)
        s     : 点1からの距離(m)
*/
     var p2                         //求める点(点2)の緯度・経度(google.maps.LatLngクラス)
     var a = 6378137.0;             //長半径(赤道半径)
     var f = 1.0 / 298.257223563;   //扁平率(WGS84)
     var b = a * (1 - f);           //短半径(極半径)
     var phi1,phi2;                 //点1、点2の緯度
     var L1,L2;                     //点1、点2の経度
     var L;                         //点1、点2の経度差
     var alpha1,alpha2;
     var alpha;
     var U1,U2;
     var lambda
     var sigma;
     var sigma1,formersigma;
     var sigmam,deltasigma;
     var epsilon = 1.0e-12;         //計算精度(sigmaの許容誤差)
     var sinalpha,cosalpha2;
     var u2,A,B,C;
     var cnt=0;

     phi1 = calRadian(p1.lat());
     L1 = calRadian(p1.lng());
     alpha1 = calRadian(deg1);

     U1 = Math.atan((1 - f)*Math.tan(phi1));
     sigma1 = Math.atan2(Math.tan(U1),Math.cos(alpha1));
     sinalpha = Math.cos(U1)*Math.sin(alpha1);
     cosalpha2 = 1-Math.pow(sinalpha,2);
     u2 = cosalpha2*((Math.pow(a,2)-Math.pow(b,2))/Math.pow(b,2));
     A = 1+u2/16384*(4096+u2*(-768+u2*(320-175*u2)));
     B = u2/1024*(256+u2*(-128+u2*(74-47*u2)));

     sigma = s/(b*A);
     do {
        sigmam = (2*sigma1+sigma)/2;
        deltasigma = B*Math.sin(sigma)
                       *(Math.cos(2*sigmam)
                         +1/4*B*(Math.cos(sigma)*(-1+2*Math.pow(Math.cos(2*sigmam),2))
                                 -1/6*B*Math.cos(2*sigmam)*(-3+4*Math.pow(Math.sin(sigma),2))*(-3+4*Math.pow(Math.cos(2*sigmam),2))));
        formersigma = sigma;
        sigma = s/(b*A)+deltasigma;
        cnt++;
        if (cnt>100) {
           alert("error");
           return null;
        }
     } while (Math.abs(sigma-formersigma)>epsilon);

     phi2 = Math.atan2((Math.sin(U1)*Math.cos(sigma)+Math.cos(U1)*Math.sin(sigma)*Math.cos(alpha1))
                     ,((1-f)*Math.sqrt(Math.pow(sinalpha,2)+Math.pow((Math.sin(U1)*Math.sin(sigma)-Math.cos(U1)*Math.cos(sigma)*Math.cos(alpha1)),2))));
     lambda = Math.atan2(Math.sin(sigma)*Math.sin(alpha1),(Math.cos(U1)*Math.cos(sigma)-Math.sin(U1)*Math.sin(sigma)*Math.cos(alpha1)));
     C = f/16*cosalpha2*(4+f*(4-3*cosalpha2));
     L = lambda-(1-C)*f*sinalpha*(sigma+C*Math.sin(sigma)*(Math.cos(2*sigmam)+C*Math.cos(sigma)*(-1+2*Math.pow(Math.cos(2*sigmam),2))));
     L2 = L + L1;
     alpha2 = Math.atan(sinalpha/(-1*Math.sin(U1)*Math.sin(sigma)+Math.cos(U1)*Math.cos(sigma)*Math.cos(alpha1)));
     p2 = new google.maps.LatLng(calDegree(phi2),calDegree(L2));
     return p2;
  }

 アークタンジェントの関数はMath.atan()とMath.atan2()の2種類あり、これらを使い分ける必要がありますが、Wikipediaではその辺も親切に解説されています。
 Vincenty法での緯度・経度を計算するサイト(数箇所)での計算結果と比較することで検証しましたが、誤差はほとんどなかったので上記ソースで多分大丈夫かと思います。なお、扁平率はWGS84のものを採用しましたが、GRS80の扁平率を採用しても計算結果はほとんど変わりません。

 と、ここまで作成した段階で残念な発見・・・同様の計算をする関数がgoogle.maps.geometryライブラリに用意されていました。

 ・Maps Javescript API ジオメトリ ライブラリより抜粋
 computeOffset() を使用すると、特定の角度、出発地点、移動距離(メートル単位)から、目的地の座標を計算できます。

 Google Maps APIを読み込む際に、

<script type="text/javascript" src="http://maps.google.com/maps/api/js?language=jp&key=取得したAPIキー&libraries=geometry"></script>

 といった具合にlibraries=geometryのパラメータを追記するだけで利用可能です。google.maps.Circleも当関数を利用しているようです。
 computeOffset()と自作の上記関数とで誤差があるか検証したところ、低緯度の鉛直方向(北半球なら南側)の誤差が大きくなるようで、ハボロネ(ボツワナの首都)から5000km地点計算時には約31kmのズレが出ました。5000kmで31kmならば誤差は1%未満、地図上ではほとんど分からないレベルですが、原因が気になるのでいろいろ調べてみるとgoogle.maps.geometryライブラリでは地球を完全な球体とみなして計算していることが分かりました(自作関数の扁平率をゼロにして再度検証した結果は両者一致しました)。
 大差ないとはいえ、自作関数の方が扁平率を考慮している分だけより正確ということになるのか・・・わざわざ作った甲斐がありました。computeOffset()ではなく、こちらを採用します。
 
 ここまで出来てしまえば、地図に落とす処理は簡単です。
 今回のテスト用html、Javascriptのソースを以下に記載します。 

・concentriccircles_test2.html

<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN" "http://www.w3.org/TR/html4/loose.dtd">
<html>
<head>
<title>concentric circles on the map</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="http://maps.google.com/maps/api/js?language=jp&libraries=geometry"></script>
<style type="text/css">
<!--
.circletitle {
   font-size:12px;
   color:#FFFF00;
   padding-left:20px;
   white-space: nowrap; 
}
.circletitle input {
   width:40px;
}
.circletitle input[type="button"] {
   width:60px;
   font-size:14px;
}

-->
</style>
<script type="text/javascript" src="js/concentriccircles_test2.js"></script>
</head>

<body style="margin:0px;" onload="initialize()">
  <div id="top_bar" style="background-color:#333399;padding:5px;">
  <table>
  <tr>
    <th style="color:#FFFFFF;text-align:left;padding-left:20px;" colspan="7" style="font-size: 16px;">中心となる地点から一定距離にある点を求めて地図上に同心円を描く(テスト2)~Vincentyの順解法利用~</th>
  </tr>
  <tr>
    <td class="circletitle">circle1</td>
    <td class="circletitle">circle2</td>
    <td class="circletitle">circle3</td>
    <td class="circletitle">circle4</td>
    <td class="circletitle">circle5</td>
    <td colspan="2">&nbsp;</td>
  </tr>
  <tr>
    <td class="circletitle">
        <input type='text' id="f_distance0">km
        <input type='text' id="f_color0">
    </td>
    <td class="circletitle">
        <input type='text' id="f_distance1">km
        <input type='text' id="f_color1">
    </td>
    <td class="circletitle">
        <input type='text' id="f_distance2">km
        <input type='text' id="f_color2">
    </td>
    <td class="circletitle">
        <input type='text' id="f_distance3">km
        <input type='text' id="f_color3">
    </td>
    <td class="circletitle">
        <input type='text' id="f_distance4">km
        <input type='text' id="f_color4">
    </td>
    <td class="circletitle">
       <input type="button" id="btnDrawCircle" name="btnDrawCircle" onclick="drawCircle()" value="draw">
    </td>
    <td class="circletitle">
       <input type="button" id="btnResetCircle" name="btnResetCircle" style="color:#ff0000;" onclick="resetCircle()" value="reset">
    </td>
  </tr>
  </table>
  </div>
  <div id="map_canvas" style="margin:5px;height:600px"></div>
</body>
</html>

・concentriccircles_test2.js

  var map;
  var svg;
  var Gaborone = new google.maps.LatLng(-24.6091799,25.7900739);
  var iniZoom = 4;
  var pointMarker = new google.maps.Marker();
  var circleInfo = new Array();
  var myCircle = new Array();
  var polygons = new Array();
  var drawFlag;

  function initialize() {

     drawFlag = false;
     var dColor = ["blue","green","red","yellow","magenta"];
     var dDistance = [500,1000,1500];

     for (var i = 0; i <= 4; i++) {
        if (dDistance[i]) {
           document.getElementById("f_distance"+i).value = dDistance[i];
        }
        if (dColor[i]) {
           document.getElementById("f_color"+i).value = dColor[i];
        }
     }

     var myOptions = {
        zoom: iniZoom,
        center: Gaborone,
        mapTypeControl: false,
        mapTypeId: google.maps.MapTypeId.ROADMAP
     };
     map = new google.maps.Map(document.getElementById("map_canvas"), myOptions);

     google.maps.event.addListener(map, "click", function(clickevent) {
        var myPoint = new google.maps.LatLng(clickevent.latLng.lat(), clickevent.latLng.lng());
        setPointMarker(myPoint);
     });

  }

  function setPointMarker(location, pile) {

    if(!pointMarker.getMap()){
       pointMarker = new google.maps.Marker({
          position: location,
          map: map,
          draggable: true,
       });
    }else{
       pointMarker.setPosition(location);
    }

  }

  function setCircleInfo() {
     
     var curDistance,curColor,curPplText;

     circleInfo = new Array();

     for (var i = 0; i <= 4; i++) {
        curDistance = parseFloat(document.getElementById("f_distance"+i).value)*1000;
        curColor = document.getElementById("f_color"+i).value;
        if (!curColor) {
           curColor = "#ffffff";
        }
        if (isFinite(curDistance)) {
           circleInfo.push({"distance":curDistance,"color":curColor});
        }
     }

     circleInfo.sort(function(a,b){
        if(a.distance<b.distance) return -1;
        if(a.distance > b.distance) return 1;
        return 0;
     });

  }

  function resetCircle() {

     var i;

     for (var i = 0; i < myCircle.length; i++) {
        if (myCircle[i]) {
           myCircle[i].setMap(null);
        }
     }
     myCircle.length = 0;
     
     if (pointMarker.getMap()) {
        pointMarker.setMap(null);
     }
     drawFlag = false;
  }

  function drawCircle() {

     var polyCoords = new Array();
     var polyPath = new Array();
     var degpitch = 0.5;//0.5度ピッチで計算

     if (!pointMarker.getMap()) {
        alert("中心となる点にマーカーを設置して下さい。");
        return 0;
     }

     if (drawFlag) {
        return 0;
     }

     var myPoint = pointMarker.getPosition();

     setCircleInfo();

     for (var i = 0; i < circleInfo.length; i++) {
        var polyCoord = new Array();
        //中心点から指定の距離だけ離れた地点の座標を計算し、配列に収納(ポリゴン用)
        for (var theta = 0; theta < 360; theta += degpitch) {  //thetaは真北からの角度(時計回り)
           polyCoord.push(VincentyDirect(myPoint,theta,circleInfo[i].distance));
        }
        polyCoord.push(polyCoord[0]);
        polyCoords[i] = polyCoord;

        polyPath = [polyCoords[i]];
        //最内以外の円描画時は、1つ内側の円の座標を反対周りにした配列を付加することでドーナツ型ポリゴンにする
        if (i > 0) {
           polyPath.push(polyCoords[i-1].slice().reverse());
        }

        myCircle[i] = new google.maps.Polygon({
           paths: polyPath,
           strokeWeight: 0,
           fillColor: circleInfo[i].color,
           fillOpacity: 0.4
        });
        myCircle[i].setMap(map);

     }

     drawFlag = true;

  }

  function calRadian(deg){
     return deg * Math.PI / 180.0;
  }

  function calDegree(rad){
     return rad * 180 / Math.PI;
  }

  //Vincenty法(順解法)
  function VincentyDirect(p1,deg1,s) {
/*
     	引数
        p1    : 中心となる点(点1)の緯度・経度(google.maps.LatLngクラス)
        deg1  : 求める点の方角(真北からの時計回りの角度。北:0、東:90、南:180、西:270)
        s     : 点1からの距離(m)
*/
     var p2                         //求める点(点2)の緯度・経度(google.maps.LatLngクラス)
     var a = 6378137.0;             //長半径(赤道半径)
     var f = 1.0 / 298.257223563;   //扁平率(WGS84)
     var b = a * (1 - f);           //短半径(極半径)
     var phi1,phi2;                 //点1、点2の緯度
     var L1,L2;                     //点1、点2の経度
     var L;                         //点1、点2の経度差
     var alpha1,alpha2;
     var alpha;
     var U1,U2;
     var lambda
     var sigma;
     var sigma1,formersigma;
     var sigmam,deltasigma;
     var epsilon = 1.0e-12;         //計算精度(sigmaの許容誤差)
     var sinalpha,cosalpha2;
     var u2,A,B,C;
     var cnt=0;

     phi1 = calRadian(p1.lat());
     L1 = calRadian(p1.lng());
     alpha1 = calRadian(deg1);

     U1 = Math.atan((1 - f)*Math.tan(phi1));
     sigma1 = Math.atan2(Math.tan(U1),Math.cos(alpha1));
     sinalpha = Math.cos(U1)*Math.sin(alpha1);
     cosalpha2 = 1-Math.pow(sinalpha,2);
     u2 = cosalpha2*((Math.pow(a,2)-Math.pow(b,2))/Math.pow(b,2));
     A = 1+u2/16384*(4096+u2*(-768+u2*(320-175*u2)));
     B = u2/1024*(256+u2*(-128+u2*(74-47*u2)));

     sigma = s/(b*A);
     do {
        sigmam = (2*sigma1+sigma)/2;
        deltasigma = B*Math.sin(sigma)
                       *(Math.cos(2*sigmam)
                         +1/4*B*(Math.cos(sigma)*(-1+2*Math.pow(Math.cos(2*sigmam),2))
                                 -1/6*B*Math.cos(2*sigmam)*(-3+4*Math.pow(Math.sin(sigma),2))*(-3+4*Math.pow(Math.cos(2*sigmam),2))));
        formersigma = sigma;
        sigma = s/(b*A)+deltasigma;
        cnt++;
        if (cnt>100) {
           alert("error");
           return null;
        }
     } while (Math.abs(sigma-formersigma)>epsilon);

     phi2 = Math.atan2((Math.sin(U1)*Math.cos(sigma)+Math.cos(U1)*Math.sin(sigma)*Math.cos(alpha1))
                     ,((1-f)*Math.sqrt(Math.pow(sinalpha,2)+Math.pow((Math.sin(U1)*Math.sin(sigma)-Math.cos(U1)*Math.cos(sigma)*Math.cos(alpha1)),2))));
     lambda = Math.atan2(Math.sin(sigma)*Math.sin(alpha1),(Math.cos(U1)*Math.cos(sigma)-Math.sin(U1)*Math.sin(sigma)*Math.cos(alpha1)));
     C = f/16*cosalpha2*(4+f*(4-3*cosalpha2));
     L = lambda-(1-C)*f*sinalpha*(sigma+C*Math.sin(sigma)*(Math.cos(2*sigmam)+C*Math.cos(sigma)*(-1+2*Math.pow(Math.cos(2*sigmam),2))));
     L2 = L + L1;
     alpha2 = Math.atan(sinalpha/(-1*Math.sin(U1)*Math.sin(sigma)+Math.cos(U1)*Math.cos(sigma)*Math.cos(alpha1)));
     p2 = new google.maps.LatLng(calDegree(phi2),calDegree(L2));
     return p2;
  }

 各距離の円の色が重ならないよう、外側の円のポリゴンはドーナツ状にしています。Google Maps APIのポリゴンでは、内側領域のデータを外側とは逆回りの配列で追記することで中抜けの形を描画することが出来るので、外側の円の描画時にはその1つ内側のポリゴン配列を逆回りにしてPolygonクラスにデータを渡しています(上記ソース125~129行目)。

 そして、上のHTMLを開いて描画させた結果が以下の図です。

計算・描画結果

計算・描画結果 ボツワナの首都(ハボロネ)から1000km,2000km,3000kmの同心円


 
 同心円の色も重ならず、また半透明なので下の地図が見えなくなってしまうこともなくバッチリです。
 「これで完成!」と思ったのですが・・・・・・。

 距離を10000kmなど大きくして試してみると・・・

計算・描画結果(長距離)

計算・描画結果 ボツワナの首都(ハボロネ)から3000km,5000km,10000kmの同心円

 このように一番外側の距離円(もはや円ではありませんが)のポリゴンが本来とは逆の領域を選択してしまっています。その関係で内側の円には外側の円の色(今回の場合は赤)が重なってしまいました。
 中心点や距離をいろいろ変えて描画させてみたところ、どうやらGoogle Maps APIのポリゴンは地球を1周して南北に2分割するような大きな領域を描く際には北半球側を選択するようです。

 ところが悔しいことに、google.maps.Circleで同条件の同心円を描いた場合は以下の図のように南半球でもきちんと選択されています

google.maps.Circle

同条件でにてgoogle.maps.Circleを利用した結果


 「色が重なる程度の小さなことは気にせず、おとなしくgoogle.maps.Circleを使いなさい!」ということなんでしょうか?・・・いえいえ、まだ諦めません。次回へつづく

コメント   

 久々に本職のシステム関連の投稿を・・・。

 ある日、配属先の同僚から次のような相談を受けました。
「今、ここに首都(ハボロネ)からの距離が500km、1000km、1500kmの同心円になっている地図があるんだけど、この中心をハボロネに限らず任意の点から作図出来るようなアプリ作れない?」
 地図上の距離に関しては以前に「緯度・経度情報のある2点間の距離」を計算する関数を作成しているので「ああ、出来るよ!」と気軽に答えたものの・・・細部にこだわったこともあり(相手はそこまで求めていないのですが)、結構苦労しました。

 まず、地図に関しては中心とする地点の選びやすさを考慮してGoogle Maps APIを利用することにしました。Google Mapにはgoogle.maps.Circleという便利なクラスがあり、これを使えば簡単に同心円が作れます。なので、円を描くだけならそれでいいのです。
 しかし、見本には各円に色が塗ってある。もちろんCircleクラスでも円を塗りつぶすことが出来ますが、複数の円を描く同心円では色が重なってしまいます。透過度をゼロにして内側から描けば色は重なりませんが、今度は地図が見えなくなってしまい場所が何処だか分からなくなってしまうので・・・。結局自分で計算することにしました。

 「点Aからの距離がrである点の座標(緯度・経度)」を求める計算(点Bは無数に存在)。最初考えたのはAから真東にだけ離れた点(同緯度)と真北にrだけ離れた点(同経度)を最初に求めて、その2点を通る楕円上の点の座標を求める方法です(下のイメージ図)。

点Aからの距離がrである点の座標(緯度・経度)を求める計算手法(1)~イメージ~

点Aからの距離がrである点の座標(緯度・経度)を求める計算手法(1)~イメージ~

 以前作成した2点間の距離を計算する関数を使えば、点Aから緯度1度分離れた点及び経度1度分離れた点の距離が計算出来るので、それを基準に上図の点B・点Cの位置情報を求め、後は三角関数を使用して楕円上の点を計算していくというものです。実際、ネットで調べてみると同じような理論で計算されている方もちらほらといらっしゃるので、これでOKかな?と目論んでいたものの・・・。
 以下のコードが当手法で計算しポリゴンとして地図上に落としたものです。

  var map;
  var mapCenter = new google.maps.LatLng(-23, 23);
  var pointMarker = new google.maps.Marker();
  var polygons = new Array();

  function initialize() {

     var myOptions = {
        zoom: 3,
        center: mapCenter,
        mapTypeControl: true,
        mapTypeId: google.maps.MapTypeId.ROADMAP
     };
     map = new google.maps.Map(document.getElementById("map_canvas"), myOptions);

     google.maps.event.addListener(map, "click", function(clickevent) {
        var myPoint = new google.maps.LatLng(clickevent.latLng.lat(), clickevent.latLng.lng());
        setPointMarker(myPoint);
     });


  }

  function setPointMarker(location, pile) {

    if(!pointMarker.getMap()){
       pointMarker = new google.maps.Marker({
          position: location,
          map: map,
          draggable: true,
       });
    }else{
       pointMarker.setPosition(location);
    }
  }

  function drawCircle() {

     if (!pointMarker.getMap()) {
        alert("中心となる点にマーカーを設置して下さい。");
        return 0;
     }

     var myPoint = pointMarker.getPosition();
     var pN = new google.maps.LatLng(myPoint.lat()-0.5,myPoint.lng()); //中心点から北へ0.5度
     var pS = new google.maps.LatLng(myPoint.lat()+0.5,myPoint.lng()); //中心点から南へ0.5度
     var pW = new google.maps.LatLng(myPoint.lat(),myPoint.lng()-0.5); //中心点から西へ0.5度
     var pE = new google.maps.LatLng(myPoint.lat(),myPoint.lng()+0.5); //中心点から東へ0.5度
     var dLat = getDistance(pN,pS); //緯度1度分の距離
     var dLng = getDistance(pW,pE); //距離1度分の距離
     var curLat,curLng;
     var cLat,cLng;
     var polyCircle = new Array;
     var r = 3000000;//今回は3000kmで計算
     var degpitch = 0.5; //今回は同心円上の点は0.5度ピッチで計算

     for (var theta = 0; theta < 360; theta += degpitch) {
        curLng = myPoint.lng() + (r / dLng) *  Math.cos(calRadian(theta));  //経度(X座標)
        curLat = myPoint.lat() + (r / dLat) *  Math.sin(calRadian(theta));  //緯度(Y座標)
        polyCircle.push(new google.maps.LatLng(curLat,curLng));
     }
     polyCircle.push(polyCircle[0]);

     var myCircle = new google.maps.Polygon({
        paths: polyCircle,
        fillColor: '#00FFFF',
        fillOpacity: 0.5,
        strokeColor: '#00FFFF',
        strokeOpacity: 1.0,
        strokeWeight: 1
     });     
     myCircle.setMap(map);

     //比較用 Google Maps API Circleクラスでの描画
     var myCircleB = new google.maps.Circle({
        map: map,
        center: myPoint,
        radius: r,
        fillOpacity: 0,
        strokeColor: "#FF0000",
        strokeOpacity: 1,
        strokeWeight: 2
     });

  }

  //度数→ラジアン変換
  function calRadian(deg){
     return deg * Math.PI / 180.0;
  }

  //2点間の距離を計算(GRS80)
  function getDistance(p1, p2) {
     var a = 6378137.0;             //長半径
     var f = 1.0 / 298.257222101;   //扁平率
     var b = a * (1 - f);           //短半径
     var e2 = (Math.pow(a, 2) - Math.pow(b, 2)) / Math.pow(a, 2);
     var dy = calRadian(p1.lat() - p2.lat());
     var dx = calRadian(p1.lng() - p2.lng());
     var my = calRadian((p1.lat() + p2.lat()) / 2.0);
     var w = Math.sqrt(1.0 - e2 * Math.pow(Math.sin(my), 2));
     var m = (a * (1 - e2)) / Math.pow(w, 3);
     var n = a / w;

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

 一見良さげなのですが、Circleクラスで描いた円(下図の赤い線)と比較すると結構ズレがあります。
 真東・真西の位置はぴったり合っているようなのですが、特に南東・南西方向のズレは大きく100km以上です(下の拡大図参照)。これだと都市地図等の近距離なら気にならない程度の誤差で収まるかもしれませんが、世界地図レベルとなると実用には耐えられそうもありません。

計算・描画結果

計算・描画結果 ボツワナの首都(ハボロネ)から3000kmの円
Google Maps APIのCircleクラス利用で描いたもの(赤い線)とは結構ズレがあります。
この地図を別窓で表示


誤差の大きい部分を拡大(中心から南東方向)

ズレの大きい部分を拡大(中心から南東方向)。100km以上の誤差があります。これではちょっと・・・


 考えてみたら、経度1度当たりの距離は高緯度になるほど短くなるのに、それを考慮せずに1度当たりの距離を過大評価してしまっている状態ですから、等距離の円も高緯度では実際よりも萎んだ形となってしまうのも当然です。
 また、地球は真球ではなく楕円体なので同じ「緯度1度分」であっても北へ1度と南へ1度とでは距離が違うのですね。つまり、緯度差と距離は比例しないので比率換算では緯度は求められない(上記の方法はそこが間違い)。
 参考までに、上記の2点間の距離を計算する関数を利用して、東京から北へ10度離れた点と南へ10度離れた点までの距離を調べてみたところ、

 ・東京(35.69,139.69)から北へ1度離れた点までの距離:1110479.2m
 ・東京(35.69,139.69)から南へ1度離れた点までの距離:1108641.5m

 2km近い距離差となりました。意外だったのは北側の方が長かったことです。地球は赤道半径の方が極半径より長いので、「低緯度側(北半球の場合は南側)の方が当然長くなるだろうと計算前には思っていたのですが結果は逆でした。プログラムが間違っているのかとも思いましたが、国土地理院のサイトでの計算でもほぼ同じ結果を得られたので間違いないようです。

 緯度は地球の中心からの角度で決めたもの(地心緯度)ではなく、楕円体面の法線と赤道の角度で決めたもの(地理緯度)だそうで、そのため高緯度ほど緯度1度当たりの距離が長くなるようで・・・、私にはなんだか難しいですが、いずれにしても上記の方法では駄目なので、考え直しが必要です。次回へつづく

コメント   

 前回までで準備は整ったので、後は次の式に従って予測計算し地図出力するだけです。

クリギング計算式

クリギング計算式


 上式の重み係数(w)は次の方程式を解く形で求めます。
重み係数を求める方程式

重み係数を求める方程式


 この式の右辺をBとおくと、重み係数(w)は w=A-1×B で求まります。
 逆行列の計算は既にレーベンバーグ・マーカート法の中で実装済み、セミバリオグラムを求める近似式も前々回で求めているので、計算は簡単なはず・・・なんですが、実際はかなり難航しました。

 セミバリオグラムのモデルはガウスモデルを選択してテストしていたのですが、予測結果はあり得ないような数値になってしまいます。
 参考書籍(地理空間データ分析(Rで学ぶデータサイエンス 7)共立出版)のデータを使用した際には正常な結果が得られるので、プログラムが根本的に間違っているというわけではなさそうなのですが・・・。

 細かく検証した所、ガウスモデル選択時には逆行列(A-1)の計算に失敗していることが判明しました。
 下図はセミバリオグラム近似曲線の一部を拡大したものですが、ガウスモデルは他モデルと比較して距離ゼロ付近ではセミバリオグラムの上がり方が緩やかになります。

ガウスモデル

ガウスモデル

円形モデル

円形モデル

 
図1 セミバリオグラム近似曲線の比較(原点付近を拡大)

 このため距離が非常にゼロに近いペアが存在する場合、そのセミバリオグラムは他のペアに対して極端に小さい数値となるため、逆行列の計算過程において有効桁数(16桁)では処理しきれなくなって異常な値となってしまうものと思われます。
 それならば「逆行列の計算をしない方法を採用すればどうか?」とLU分解にて解を求める方法も試みましたが、やはり芳しくない結果となってしまいます。逆行列・LU分解ともにいくつかの計算手法が存在するので、今回試した方法以外で上手く計算出来るものがあるかもしれませんが・・・。
 ガウスモデルを採用する場合は、今回のように逆行列の計算に失敗したり、重み係数が異常値(大きなマイナス数値など)になったりし、良い結果が得られないことが時々あるようで、ネット検索してみるとコミュニティサイトや論文(いずれも英語ですが)にそういった記述がいくつか見つかりました。
 おそらく観測点を絞るなどして再計算すれば上手くいくのでしょうが、他のモデルでは正常と思われる結果が得られたので、今回はガウスモデルを除外することにしました。
 
 計算に結構時間がかかるので地図出力とはプログラムを分け、予測結果を一旦csv出力した後に、そのcsvデータを読み込んで出力する形式にしました。最終的なプログラムソースは以下の通り。

・MinamiKantoAirPollutionCalKriging.html(計算用)

<!DOCTYPE html>
<html>
<head>
<title>南関東(1都3県)の大気汚染(NO2)濃度 クリギング法による予測計算</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="js/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/MinamiKantoAirPollutionCalKriging.js"></script>
<script type="text/javascript" src="js/LevenbergMarquardtKriging.js"></script>
</head>
<body style="margin:10px;" onload="initialize()">南関東(1都3県)の大気汚染(NO2)濃度 クリギング法による予測計算
</body>
</html>

 

・MinamiKantoAirPollutionCalKriging.js

// *************************************************************** 
//  D3.jsで空間分析にチャレンジ クリギング
// *************************************************************** 

  var obspointFileName = "csv/T_Place_AirPollutionMinamiKanto.csv";//大気汚染常時監視測定局情報データ(緯度・経度・NO2濃度)
  var calpointFileName = "csv/CalPoint.csv";                       //クリギング予測対象点情報データ
  var resultFileName = "csv/EstimatedResult.csv";                  //クリギング予測結果(出力ファイル)
  var calpoint;
  var samplepoint = new Array();
  var mathmodels = new Array();

  function initialize() {
     queue()
        .defer(d3.csv, obspointFileName)
        .defer(d3.csv, calpointFileName)
        .await(Kriging);
  }

  function Kriging(error, data1, data2) {

     var dRange = 1.0e+20;
     var uvec;
     var dInterval;
     var i,j,k,c;
     var x,dd;
     var dist = new Array();
     var svar = new Array();
     var disvar = new Array();
     var dsum = new Array();
     var vsum = new Array();
     var vbind = new Array();
     var cntdata = new Array();
     var cntbind = 50;
     var modelName;
     var mathmodel;

     i = 0;
     data1.forEach(function(d) {
        if (d.NO2 > 0) {
           samplepoint[i] = d;
           samplepoint[i].NO2ppb = parseFloat(d.NO2)*1000; //ppm*1000=ppb
           i++;
        }
     });

     calpoint = data2;

     for (i = 0; i < samplepoint.length - 1; i++) {
        for (j = i + 1; j < samplepoint.length; j++) {
           dd = getDistance(samplepoint[i],samplepoint[j])/1000; //m→km
           if (dd < dRange) {
              dist.push(dd);
              svar.push(Math.pow(samplepoint[i].NO2ppb - samplepoint[j].NO2ppb, 2) / 2);
           }
        }
     }

     dInterval = (d3.max(dist)-d3.min(dist)) / (cntbind - 1);
     uvec = d3.range(d3.min(dist),d3.max(dist)+1.0e-10,dInterval);

     for (c = 0; c < uvec.length; c++) {
        cntdata1 = 0;
        dsum1 = 0;
        vsum1 = 0;
     }

     for (i = 0; i < dist.length - 1; i++) {
        for (c = 0; c < uvec.length; c++) {
           if (dist[i] >= uvec1 - (dInterval / 2) && dist[i] < uvec1 + (dInterval / 2)) {
              vsum1 += svar[i];
              cntdata1 ++;
              disvar.push({'distance':dist[i], 'semivariance':svar[i], 'distanceclass':c});
           }
        }
     }

     for (var c = 0; c < uvec.length; c++) {
        if (cntdata1 > 0) {
           vbind1 = vsum1 / cntdata1;
        } else {
           vbind1 = 0;
        }
     }

     //モデルの選択
     modelName = "spherical";initialP = [20,100];    //球形モデル
//     modelName = "exponential";initialP = [20,100];  //指数モデル
//     modelName = "gaussian";initialP = [20,100];     //ガウスモデル
//     modelName = "power";initialP = [0.05,1.5];      //べき乗モデル
//     modelName = "linear";initialP = [1];            //線形モデル
//     modelName = "nugget";initialP = [20];           //ナゲット効果モデル
     mathmodel = setModel(modelName,initialP);

     LevenburgMarquardtKriging(uvec,vbind,mathmodel);

     var fct = mathmodel.fct
     var fctsub = mathmodel.fctsub
     var rangedist = 1.0e+20;

     if (mathmodel.distpNo > 0) {
        rangedist = mathmodel.parameters[mathmodel.distpNo];
     }
     for (i = 0; i < mathmodel.parameters.length; i++) {
        fct = fct.replace(new RegExp('p\\['+i+'\\]','g'), mathmodel.parameters[i]);
        fctsub = fctsub.replace(new RegExp('p\\['+i+'\\]','g'), mathmodel.parameters[i]);
     }
     fct = fct.replace(new RegExp(' ','g'), '');
     fct = fct.replace(new RegExp('--','g'), '+');
     fctsub = fctsub.replace(new RegExp(' ','g'), '');
     fctsub = fctsub.replace(new RegExp('--','g'), '+');

     var g = new Array();
     var gm = new Array();
     var w = new Array();

     for (var i = 0; i <= samplepoint.length; i++) {
        g[i] = new Array();
        for (var j = 0; j <= samplepoint.length; j++) {
           if (i == samplepoint.length || j == samplepoint.length) {
              g[i][j] = 1;
              if (i == samplepoint.length && j == samplepoint.length) {
                g[i][j] = 0;
              }
           } else {
              x = getDistance(samplepoint[i],samplepoint[j])/1000; //m→km
              if (x >= rangedist) {
                 g[i][j] = eval(fctsub);
              } else {
                 g[i][j] = eval(fct);
              }
           }
        }
     }

     var invg = InverseMatrix(g,true);

     if (!checkInverseMatrix(g,invg)) {
        alert("逆行列の計算に失敗しました。");
        return false;
     }

     for (i = 0; i < calpoint.length; i++) {
        calpoint[i].EstimatedValue = 0;
        gm = new Array();
        for (var j = 0; j <= samplepoint.length; j++) {
           if (j == samplepoint.length) {
              gm[j] = 1;
           } else {
              x = getDistance(calpoint[i],samplepoint[j])/1000;
              if (x >= rangedist) {
                 gm[j] = eval(fctsub);
              } else {
                 gm[j] = eval(fct);
              }
           }
           w[j] = 0;
        }

        for (var j = 0; j <= samplepoint.length; j++) {
           if (modelName == "nugget") {
              w[j] = 1 / samplepoint.length;
           } else {
              for (var k = 0; k <= samplepoint.length; k++) {
                 w[j] += invg[k][j] * gm[k];
              }
           }
        }

        for (var j = 0; j < samplepoint.length; j++) {
           calpoint[i].EstimatedValue += w[j] * samplepoint[j].NO2ppb;
        }

     }

     saveResult(modelName,calpoint);

  }

  //度数→ラジアン変換
  function calRadian(deg){
     return deg * Math.PI / 180.0;
  }

  //2点間の距離を計算(GRS80)
  function getDistance(p1, p2) {
     var a = 6378137.0;             //長半径
     var f = 1.0 / 298.257222101;   //扁平率
     var b = a * (1 - f);           //短半径
     var e2 = (Math.pow(a, 2) - Math.pow(b, 2)) / Math.pow(a, 2);
     var dy = calRadian(parseFloat(p1.Lat) - parseFloat(p2.Lat));
     var dx = calRadian(parseFloat(p1.Lng) - parseFloat(p2.Lng));
     var my = calRadian((parseFloat(p1.Lat) + parseFloat(p2.Lat)) / 2.0);
     var w = Math.sqrt(1.0 - e2 * Math.pow(Math.sin(my), 2));
     var m = (a * (1 - e2)) / Math.pow(w, 3);
     var n = a / w;

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

  //逆行列が正しいかをチェック
  function checkInverseMatrix(matrix,invmatrix) {

     var elem; // = new Array();
     var errRate = 0.01;
     var idm;
     var i,j,k;

     for (i = 0; i < matrix.length; i++) {
        for (j = 0; j < matrix.length; j++) {
           elem = 0;
           for (k = 0; k < matrix.length; k++) {
              elem += matrix[i][k]*invmatrix[k][j];
           }
           if (i == j) {
              idm = 1;
           } else {
              idm = 0;
           }
           if (Math.abs(elem - idm) > errRate) {
              return false;
           }
        }
     }
     
     return true;

  }

  //モデルの選択
  function setModel(modelname,initialP) {

     var modelName;   //モデル名
     var strFct;      //バリオグラム関数
     var strFctSub;   //バリオグラム関数2(パラメータ値によって式が切り替わる場合のみ設定。球型モデルなど)
     var distpNo;     //上記バリオグラム関数2を適用する場合、切替対象となるパラメータ番号。適用しない場合は-1

     //球形モデル
     if (modelname == "spherical") {
        modelName = "spherical model";
        strFct = "p[0] * (1.5 * (x / p[1]) - 0.5 * Math.pow((x / p[1]),3))";
        strFctSub = "p[0]";
        distpNo = 1;
     //指数モデル
     } else if (modelname == "exponential") {
        modelName = "exponential model";
        strFct = "p[0] * (1 - Math.exp(-x / p[1]))";
        strFctSub = "";
        distpNo = -1;
     //ガウスモデル
     } else if (modelname == "gaussian") {
        modelName = "gaussian model";
        strFct = "p[0] * (1 - Math.exp(-Math.pow((x / p[1]),2)))";
        strFctSub = "";
        distpNo = -1;
     //べき乗モデル
     } else if (modelname == "power") {
        modelName = "power model";
        strFct = "p[0] * Math.pow(x,p[1])";
        strFctSub = "";
        distpNo = -1;
     //線形モデル
     } else if (modelname == "linear") {
        modelName = "linear model";
        strFct = "p[0] * x";
        strFctSub = "";
        distpNo = -1;
     //ナゲット効果モデル
     } else if (modelname == "nugget") {
        modelName = "nugget effect";
        strFct = "p[0]";
        strFctSub = "";
        distpNo = -1;
     }

     return new ModelInfo(modelName,strFct,strFctSub,distpNo,initialP);

  }

  //予測結果をCSVに出力
  function saveResult(modelname,cp) { 

     var jsonObj = JSON.stringify(cp);
     var param = "filename=" + resultFileName.replace(".csv", "_" + modelname + ".csv");
     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('失敗しました');
        }
     });
  }

 
 checkInverseMatrix(201~227行目)は、逆行列が正しいかをチェックをする関数で、元の行列との積が単位行列になるかどうかで判断しています。今回のデータでは、ガウスモデル選択時にはこのチェックでFALSE判定され予測計算を中断します。

・savecalpoint.php(結果のcsv出力用) ※前回掲載したものと同じ

<?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);
?>

・MinamiKantoAirPollutionMapKriging.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="js/d3.js"></script>
<script type="text/javascript" src="js/queue.v1.js"></script>
<script type="text/javascript" src="js/MinamiKantoAirPollutionMapKriging.js"></script>
</head>
<body style="margin:0px;" onload="setMap()">
  <div style="background-color:#006666;padding:5px 0px 2px 5px;">
  <table style="padding-left:5px;">
  <tr>
    <td style="font-size: 16px;font-weight: bold;color: #FFFFFF;">南関東(1都3県)の大気汚染(NO<sub>2</sub>)濃度分布 </td>
    <td style="font-size: 12px;color: #FFFFFF;">※クリギング法にて予測 </td>
    <td style="font-size: 12px;color:#FFFFFF;padding-left:10px;">バリオグラムモデル:</td>
    <td>
      <select id="f_model" name="f_model" style="width:140px;font-size:12px;" onchange="switchModel(this.value);">
         <option value="spherical" selected>球形モデル</option>
         <option value="exponential">指数モデル</option>
         <option value="power">べき乗モデル</option>
         <option value="linear">線形モデル</option>
         <option value="nugget">ナゲット効果モデル</option>
      </select>
    </td>
  </tr>
  </table>
  </div>
</body>
</html>

・MinamiKantoAirPollutionMapKriging.js

// *************************************************************** 
//  D3.jsで空間分析にチャレンジ クリギング
// *************************************************************** 

  var svg;
  var canvas;
  var mapScale = 22000;  //マップスケール
  var mapWidth = 1000;   //マップ幅
  var mapHeight = 700;   //マップ高さ
  var mapTop = 20;
  var mapPath, mapBounds;
  var cborder,pborder,outline,obspoint,estimatepoint;
  var color = d3.interpolateHsl("#66ff00", "#ff0000");
  var alpha = 0.7;
  var pMin = 3, pMax = 25;  //濃度 Min,Max(ppb)
  var jsonOutlineFileName = "json/OutlineMinamiKanto.json";        //南関東(1都3県)全体の輪郭データ
  var jsonPrefFileName = "json/PrefectureBorderMinamiKanto.json";  //南関東各都県の輪郭データ
  var jsonCityFileName = "json/CityBorderMinamiKanto.json";        //市区町村境データ
  var obspointFileName = "csv/T_Place_AirPollutionMinamiKanto.csv";//大気汚染常時監視測定局情報データ(緯度・経度・NO2濃度)
  var resultFileName = "csv/EstimatedResult.csv";                  //クリギング予測結果(出力ファイル)
  var samplepoint = new Array();
  var outlineXY = new Array();
  var mathmodels = new Array();

  function setMap() {
     var modelname = window.location.search.substring(1);
     if (modelname.length == 0) {
        modelname = "spherical";
     }
     document.getElementById('f_model').value = modelname;

     queue()
        .defer(d3.json, jsonCityFileName)
        .defer(d3.json, jsonPrefFileName)
        .defer(d3.json, jsonOutlineFileName)
        .defer(d3.csv, obspointFileName)
        .defer(d3.csv, resultFileName.replace(".csv", "_" + modelname + ".csv"))
        .await(ready1);
  }

  function switchModel(modelname) {
     queue()
        .defer(d3.csv, resultFileName.replace(".csv", "_" + modelname + ".csv"))
        .await(ready2)
  }

  function ready1(error, data1, data2, data3, data4, data5) {
     cborder = data1;
     pborder = data2;
     outline = data3;
     obspoint = data4;
     estimatepoint = data5;
     drawMap();
  }

  function ready2(error, data1) {
     estimatepoint = data1;
     drawCanvas(false);
  }

  function drawMap() {

     var i,j;
     var x,y;

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

     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")

     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])]));
     });

     i = 0;
     obspoint.forEach(function(d) {
        if (d.NO2 > 0) {
           samplepoint[i] = d;
           samplepoint[i].x = mapPath.projection()([parseFloat(d.Lng),parseFloat(d.Lat)])[0];
           samplepoint[i].y = mapPath.projection()([parseFloat(d.Lng),parseFloat(d.Lat)])[1];
           i++;
        }
     });

     svg.selectAll(".prefborder")
        .data(pborder.features)
        .enter()
        .append("path")
        .attr("d", mapPath)
        .style("fill", "none")
        .attr("class", "outline")
        .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 ObservationPoint = svg.selectAll(".point")
        .data(samplepoint)
        .enter()
        .append("circle")
        .attr("cx", function(d) {
           return d.x;
        })
        .attr("cy", function(d) {
           return d.y;
        })
        .attr("r", 2)
        .style("fill", "#0000ff")
        .on('mouseover', function(d) {
           d3.select(this).style("fill","#ff0000");
           curText = d.PointName + " " + parseFloat(d.NO2).toFixed(3);
           SetTooltip(mapPath.projection()([d.Lng,d.Lat]),curText);
        })
        .on('mouseout', function() {
           d3.select(this).style("fill","#0000ff");
           var tooltip = svg.select(".tooltip")
           if(!tooltip.empty()) {
              tooltip.style("visibility", "hidden")
           }
        })

     drawCanvas(true);

  }

  function drawCanvas(initFlag) {

     var currentcolor;
     var canvas = document.getElementById('canvas')
     var ctx = canvas.getContext('2d');

     if (initFlag) {
        ctx.beginPath();
        outlineXY.forEach(function(d,i) {
           if (i==0) {
              ctx.moveTo(d[0], d[1]);
           } else {
              ctx.lineTo(d[0], d[1]);
           }
        });
        ctx.closePath();
        ctx.clip();
     }

     ctx.clearRect(0,0,mapWidth,mapHeight); 
     ctx.globalAlpha = alpha;

     var vrange = d3.extent(estimatepoint, function (d) { return parseFloat(d.EstimatedValue); });
     var meshWidth = 0;
     var meshHeight = 0;

     for (var i = 0; i < estimatepoint.length - 1; i++) {
        if (meshWidth == 0 && estimatepoint[i].x != estimatepoint[i+1].x) {
           meshWidth = Math.abs(estimatepoint[i+1].x - estimatepoint[i].x);
        }
        if (meshHeight == 0 && estimatepoint[i].y != estimatepoint[i+1].y) {
           meshHeight = Math.abs(estimatepoint[i+1].y - estimatepoint[i].y);
        }
        if (meshWidth > 0 && meshHeight > 0) {
           break;
        }
     }

     estimatepoint.forEach(function(d) {
        currentcolor = color((parseFloat(d.EstimatedValue)-pMin)/(pMax-pMin));
        ctx.fillStyle = currentcolor;
        ctx.fillRect(d.x-0.5*meshWidth,d.y-0.5*meshHeight,meshWidth,meshHeight);
     });

     setMapLegend(initFlag);

  }

  function SetTooltip(tipXY,tipText) {

     var tooltip = svg.select(".tooltip")
     var fontSize = 12;
     var rectWidth = tipText.length * fontSize + 5;
     var rectHeight = 20;

     if(tooltip.empty()) {
        tooltip = svg
           .append("g")
           .attr("class","tooltip")

        tooltip
           .append("rect")
           .attr("height",rectHeight)
           .style("stroke","none")
           .style("fill","#ffffff")
           .style("opacity","0.8")

        tooltip
           .append("text")
           .attr("text-anchor","left")
           .attr("id","tooltiptext")
           .style("font-size",fontSize+"px")
           .style("font-family","sans-serif")
     }

     var tipleft = tipXY[0] + 2;
     if (tipleft + rectWidth > mapWidth) {
        tipleft = mapWidth - rectWidth;
     }
     var tiptop = tipXY[1] - rectHeight - 2;

     tooltip
        .style("visibility", "visible")
        .attr("transform", "translate("+tipleft+","+tiptop+")")

     tooltip.select("text")
        .text(tipText + "ppm")
        .attr("transform", "translate(5,"+(fontSize+(rectHeight-fontSize)/2)+")")

     tooltip.select("rect")
        .attr("width",rectWidth)

  }

  function setMapLegend(initFlag) { 

    var fontsize = 14;
     var legendWidth = 200;
     var legendHeight = 20;
     var legendLeft = 50;
     var legendTop = mapPath.projection()(mapBounds[0])[1] - 80;
     var valuerange = d3.range(pMin,pMax+0.0001,(pMax-pMin)/(legendWidth-1));
     var pFormat = d3.format(".3f");

     if (!initFlag) {
        var v = [pMin,pMax];
        d3.selectAll(".axistext")
         .text(function(d,i) {
              return pFormat(v[i]/1000) +"ppm";
           })
           .append("tspan")
           .attr("class","super")
           .attr("dy","-2")
           .style("font-size",(fontsize-4)+"px")
           .text("3")
        return false;
     }

     var MapLegend = d3.select("svg")
        .append("g")
        .attr("id","maplegend")
        .attr("transform", "translate("+legendLeft+","+legendTop+")")

     MapLegend
        .append("circle")
        .attr("r", 3)
        .style("fill", "#0000ff")
        .on('mouseover', function(d) {
           d3.select(this).style("fill","#ff0000");
        })
        .on('mouseout', function() {
           d3.select(this).style("fill","#0000ff");
        })

     MapLegend
        .append("text")
        .attr("id","text1")
        .style("font-size",fontsize + "px")
        .style("font-family","sans-serif")
        .style("fill","#000000")
        .attr("transform", "translate("+10+","+ (fontsize * 0.5 - 2) +")")
      .text("大気汚染常時監視測定局(一般局)")

     MapLegend
        .append("text")
        .attr("id","text2")
        .style("font-size",(fontsize - 2) + "px")
        .style("font-family","sans-serif")
        .style("fill","#000000")
        .attr("transform", "translate("+10+","+ (fontsize * 1.5) +")")
      .text("※マウスオンすると測定局の名称とNO")
        .append("tspan")
        .attr("id","sub")
        .attr("dy","2")
        .text("2")
        .style("font-size",(fontsize - 4)+"px")
        .append("tspan")
        .attr("id","text2b")
        .style("font-size",(fontsize - 2) + "px")
        .attr("dy","-2")
      .text("濃度(2014年度 年平均値)が表示されます。")

     MapLegend.selectAll(".legend")
        .data(valuerange)
        .enter()
        .append("rect")
        .attr("class", "legend")
        .attr("transform", function(d, i) {
           return "translate("+((fontsize*10)+i)+",30)";
        })
        .attr("width",1)
        .attr("height",legendHeight)
        .style("fill", function(d, i){
           return color((d - pMin) / (pMax - pMin));
        })

     MapLegend
        .append("rect")
        .attr("id","legendframe")
        .attr("transform", "translate("+(fontsize*10)+",30)")
        .attr("width",legendWidth)
        .attr("height",legendHeight)
        .style("fill", "none")
        .style("stroke","#303030")
        .style("stroke-width",0.5);

     MapLegend
        .append("text")
        .attr("id","text3")
        .style("font-size",fontsize + "px")
        .style("font-family","sans-serif")
        .style("fill","#000000")
        .attr("transform", "translate(0,"+ (30 + (fontsize + legendHeight) * 0.5) +")")
      .text("NO")
        .append("tspan")
        .attr("id","sub")
        .attr("dy","2")
        .text("2")
        .style("font-size",(fontsize - 4)+"px")
        .append("tspan")
        .attr("id","text2b")
        .style("font-size",(fontsize - 2) + "px")
        .attr("dy","-2")
      .text("濃度(年平均値)")


     var MapLegendAxis = MapLegend
        .append("g")
        .attr("transform", "translate("+(fontsize*10)+","+ (30 + legendHeight) +")")
        .attr("class","axis")
        .style("stroke","none")
        .style("fill","#303030")

     MapLegendAxis
        .append("rect")
        .attr("transform", "translate(0,10)")
        .attr("width",legendWidth)
        .attr("height",1)

     var MapLegendAxisX = MapLegendAxis.selectAll(".axis")
        .data([d3.min(valuerange)/1000,d3.max(valuerange)/1000]) //ppb→ppm
        .enter()
        .append("g")

     MapLegendAxisX
        .append("rect")
        .attr("transform", function(d,i) {
           return "translate(" + (i * legendWidth) + ",5)";
        })
        .attr("width",1)
        .attr("height",10)

     MapLegendAxisX
        .append("text")
        .attr("transform", function(d,i) {
           return "translate("+(i * legendWidth)+","+(15+(fontsize+2))+")";
        })
        .attr("text-anchor","middle")
        .attr("class","axistext")
      .attr("id",function(d,i) {
           return "axistext"+i;
        })
        .style("font-size",(fontsize-2) + "px")
        .style("font-family","sans-serif")
      .text(function(d,i) {
           return pFormat(d) +"ppm";
        })

  }

上記ソースを実行した結果、なかなかいい感じのマップが出力されました。図2は円形モデル選択時のものです。秩父の辺り(地図左上部分)が若干イメージと異なりますが、この辺は東側にしか観測点が存在せず、内挿ではなく「外挿」といった形になってしまっているので仕方がないかと思います。
 指数モデル・線形モデルの結果は地図に落としてみると円形モデルとほとんど見分けがつきません(もちろん数値は異なるのですが)。
 べき乗モデルでは少しなだらかな分布となりました(図3)。
 また、参考までにナゲット効果モデルの結果も出力しましたが、距離との相関がないモデルなので当然ながら全点均一の値となります(図4)。

図 南関東(1都3県)の大気汚染(NO)濃度分布(クリギング法による)

円形モデル

図2 円形モデル
このMAPを表示


べき乗モデル

図3 べき乗モデル
このMAPを表示


ナゲット効果モデル

図4 ナゲット効果モデル(参考)
※距離との相関がゼロのモデルなので、このような均一な分布となります。
このMAPを表示

 クリギング・・・遠い道のりでしたがようやく計算・出力出来ました。今回の結果がどの程度現状に即したたものなのかは不明ですが、見た目的にはそれっぽい地図が出力されて嬉しいです。

 今回使用のプログラム・データ 一括ダウンロード : MinamiKantoAirPollutionKriging.zip
 【収納ファイル】
 ・MinamiKantoAirPollutionCalKriging.html(計算用)
 ・MinamiKantoAirPollutionMapKriging.html(地図出力)
 ・savecalpoint.php(結果のcsv出力用)
 ・js/MinamiKantoAirPollutionCalKriging.js(計算用)
 ・js/MinamiKantoAirPollutionMapKriging.js(地図出力)
  ※以上の5つは上にソース記載
 ・js/LevenbergMarquardtKriging.js(レーベンバーグ・マーカート法による最適化計算、逆行列の計算)
 ・json/OutlineMinamiKanto.json(1都3県 全体の輪郭)
 ・json/PrefectureBorderMinamiKanto.json(1都3県 都境・県境輪郭)
 ・json/CityBorderMinamiKanto.json(1都3県 エリア内市区町村輪郭)
 ・csv/T_Place_AirPollutionMinamiKanto.csv(1都3県 大気汚染常時監視測定局観測データ 位置情報・2014年度NO2年平均値)
 ・csv/CalPoint.csv(予測地点位置情報)
 ・csv/EstimatedResult_spherical.csv(クリギング予測結果 円形モデル)
 ・csv/EstimatedResult_exponential.csv(クリギング予測結果 指数モデル)
 ・csv/EstimatedResult_power.csv(クリギング予測結果 べき乗モデル)
 ・csv/EstimatedResult_linear.csv(クリギング予測結果 線形モデル)
 ・csv/EstimatedResult_nugget.csv(クリギング予測結果 ナゲット効果モデル)
  ※各クリギング予測結果ファイルは上のプログラム(MinamiKantoAirPollutionCalKriging.js)を実行することでも得られます。

 以上のファイルは全てGNUライセンスです。ご自由に利用・加工していただいて結構ですが、計算結果の正確性は保証出来ません(不具合報告歓迎致します)。

 この他、上記プログラムの実行の際には、以下ライブラリを別途入手し、jsディレクトリに置く必要があります。
 ・jquery.js 
 ・d3.js
 ・queue.js(今回のプログラムではqueue.v1.jsという名称で保存・使用しています) 
 

コメント   

 これまで、カーネル密度推定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)が小さく、目で見た感じでも良くフィットしている印象です。
 次回は予測対象地点の設定をします。
 
 
 

 これまで、「D3.jsで空間分析にチャレンジ」の名目でボロノイ図カーネル密度推定IDW(逆距離加重補間)を作成してきましたが、次はいよいよ本丸のクリギングです。しかし、ここで完全に躓き、前回から随分と時間が開いてしまいました。クリギングはIDWとは比較にならないくらいに複雑で・・・。

 空間座標のばらつきを距離の関数として定義する「バリオグラム」という手法を用いるのですが、その関数を決定するパラメータ推定が難しい。R言語ではそれらを計算する関数が既に存在しているので、どの書籍を読んでもその手順には特に触れることもなく、当たり前のように計算しているのですが、Javascript(D3.js)ではそうはいきません。当初の目的(空間分析)からは完全に脱線しますが、最適化数学なんていう、馴染みのない分野を学習することにしました。

 遠い昔に学んだ(はずの)行列計算や偏微分の復習から始め、
 ・これなら分かる最適化数学 基礎原理から計算手法 金谷健一著 共立出版
 この本を参考に、何とか「ガウス・ニュートン法」という最適化手法までは比較的簡単にプログラム化出来たのですが、どうもこの手法は万能ではないようで関数最適化(パラメータ推定)に失敗することが多々。上記参考書にも紹介されている「レーベンバーグ・マーカート法」はガウス・ニュートン法の改良版でこれが良さそう・・・なのですが、この実装が難解を極めました。理論的にはさほど難しくないのですが、関数最適化(データと目的関数の誤差の最小化)へのプロセスが収束しているのか否かの判断が難しく試行錯誤の繰り返し。いろいろ調べてみると、RやC、Pythonなどにはレーベンバーグ・マーカート法のライブラリが存在しますが、これらのほとんどは、minpackという今から30年くらい前にFORTRANで書かれたプログラムを移植したものです。なので、これをJavascriptに移植してみることにしました。

 とりあえず移植版は完成し、数種類かのテストモデルにて試したところでは正しい結果を得られたのですが、このminpackはどういう手順で計算しているのかがさっぱりわかりません。ヤコビアン行列の計算あたりまでは、参考書にある手順と同じなのですが、その後でQR分解(?)などの処理が出てきます。その代わりに当手法(ガウス・ニュートン法もですが)では必須のはずの逆行列の計算が出てきません。逆行列を求めるのは複雑な処理なので、「ひょっとすると昔のコンピュータでは時間がかかり過ぎるためにQR分解で代用しているのかな?」なんて想像しているのですが、何しろ私には数学的バックボーンがないのでその考えが正しいのかどうかも判断分からず・・・。計算の趣旨が分からないプログラムをそのまま使用するのは納得いかないので、自作をチャレンジしてみました。

 前述しましたが、レーベンバーグ・マーカート法の実装で最も難しいのはパラメータ収束の判断です。minpackはこの部分も理解不能・・・ただ、ソース中にあるratio、actual、predictedといった言葉が関係していることは分かったので、これらのワードで検索をかけて文献を読みまくり(ほとんどが理論的なことしか記述されていません)、ついに実装出来そうな記述のある論文を見つけました。
 METHODS FOR NON-LINEAR LEAST SQUARES PROBLEMS 2nd Edition, April 2004
    K. Madsen, H.B. Nielsen, O. Tingleff
    Informatics and Mathematical Modelling Technical University of Denmark

 当文献に記載されている収束条件を採用したところ、最適計算の成功率が高まりました。パラメータ初期値の設定が適切でないと局所解に到達してしまう等で失敗する場合もありますが、それはminpackでも同じこと。テストした限りではminpackと同等程度の成功率が得られています。計算部分のソースは以下の通り。

/*
    LevenbergMarquardt.js

	Copyright (C) 2015, Yoshiro Tamura (http://www.eclip.jp/)
	All rights reserved.

    This program is free software: you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation, either version 3 of the License, or any
    later version.

    This program is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.

    You should have received a copy of the GNU General Public License
    along with this program.  If not, see <http://www.gnu.org/licenses/>.
*/

function EstimatedResult() {
    this.parameters = new Array();
    this.iteration = 0;
    this.residualsumsqueres = 0;
}

var lmresult = new EstimatedResult();

/*
  function LevenburgMarquardt

     Input Parameters
        datax : array of independent variables
        datay : array of observed variables
        p0    : array of initial guess values for parameters
        strfcn: string of model function
            The string must obtain x and p[0,1,2,...]
            e.g. p[0]*x/(p[1] + x)      (Michaelis-Menten kinetics)

     Return value
        true  : success to find a solution
        false : fail to convergence
*/
function LevenburgMarquardt(datax,datay,p0,strfcn) {

   var m = datax.length;
   var n = p0.length;
   var p = p0.concat();
   var epsilon1 = 1.0e-6;
   var epsilon2 = 1.0e-10;
   var strfcnder = new Array();
   var x;
   var i,j,k;
   var ss_lm,ss_gn,ss0;
   var iter,cntinner,cntouter;
   var itermax = 10000
   var cntinnermax = 1000;
   var cntoutermax = 1000;
   var delta_lm = new Array();
   var delta_gn = new Array();
   var p_lm = new Array();
   var p_gn = new Array();
   var g = new Array();
   var r = new Array();
   var h = new Array();
   var jac,jacT,jacTjac_lm,jacTjac_gn;
   var invjacTjac_lm,invjacTjac_gn;
   var mu = new Array();
   var muMax;
   var v = 2;
   var lambda0 = 0.0001;
   var lambda = lambda0;
   var lambdaMax = 1.0e+20;
   var lambdaMin = 1.0e-16;
   var lv = 0.01;
   var actual,predicted,ratio;
   var convergence;
   var terminateouter = false;
   var terminateinner;

   for (j = 0; j < n; j++) {
      mu[j] = 0;
   }

   convergence = false;
   iter = 0;
   cntouter = 0;

   while (!terminateouter) {
      ss0 = 0;
      for (i = 0; i < m; i++) {
         x = datax[i];
         ss0 += Math.pow(datay[i] - eval(strfcn),2);
      }
      for (j = 0; j < n; j++) {
         h[j] = Math.max(Math.abs(p[j])*epsilon1,epsilon2);
         strfcnder[j] = strfcn.replace(new RegExp('p\\['+j+'\\]','g'), '\(p\['+j+'\]+'+h[j]+'\)');
      }

      jac = new Array();
      jacT = new Array();
      jacTjac_gn = new Array();
      r = new Array();

      //Jacobian matrix
      for (i = 0; i < m; i++) {
         x = datax[i];
         jac[i] = new Array();
         r[i] = datay[i] - eval(strfcn);
         for (j = 0; j < n; j++) {
            jac[i][j] = -(eval(strfcnder[j]) - eval(strfcn)) / h[j];
         }
      }

      //transpose of Jacobian matrix
      for (j = 0; j < n; j++) {
         jacT[j] = new Array();
         g[j] = 0;
         for (i = 0; i < m; i++) {
            jacT[j][i] = jac[i][j];
            g[j] += jacT[j][i] * r[i];
         }
      }

      //Jacobian matrix transpose * Jacobian matrix
      for (j = 0; j < n; j++) {
         jacTjac_gn[j] = new Array();
         for (k = 0; k < n; k++) {
            jacTjac_gn[j][k] = 0;
            for (i = 0; i < m; i++) {
               jacTjac_gn[j][k] += jacT[j][i] * jac[i][k];
            }
         }
      }

      //copy array jacTjac_gn -> jacTjac_lm
      jacTjac_lm = JSON.parse(JSON.stringify(jacTjac_gn));

      terminateinner = false;
      cntinner = 0;
      while (!terminateinner) {
         for (j = 0; j < n; j++) {
            if (mu[j] == 0) {
               mu[j] = lambda * Math.abs(jacTjac_lm[j][j]);
            }
         }

         for (j = 0; j < n; j++) {
            jacTjac_lm[j][j] = jacTjac_gn[j][j] + mu[j];
            delta_lm[j] = 0;
            delta_gn[j] = 0;
         }

         invjacTjac_lm = InverseMatrix(jacTjac_lm);
         invjacTjac_gn = InverseMatrix(jacTjac_gn);

         p_lm = p.concat();
         p_gn = p.concat();

         for (i = 0; i < m; i++) {
            for (j = 0; j < n; j++) {
               for (k = 0; k < n; k++) {
                  delta_lm[j] -= invjacTjac_lm[j][k] * jacT[k][i]*r[i];
                  delta_gn[j] -= invjacTjac_gn[j][k] * jacT[k][i]*r[i];
               }
            }
         }

         for (j = 0; j < n; j++) {
            p_lm[j] += delta_lm[j];
            p_gn[j] += delta_gn[j];
         }

         ss_lm = 0;
         ss_gn = 0;
         for (i = 0; i < m; i++) {
            x = datax[i];
            ss_lm += Math.pow(datay[i] - eval(strfcn.replace(new RegExp('p\\[', 'g'),'p_lm\[')),2);
            ss_gn += Math.pow(datay[i] - eval(strfcn.replace(new RegExp('p\\[', 'g'),'p_gn\[')),2);
         }

         actual = (ss0 - ss_lm) / 2;
         
         predicted = 0;
         for (j = 0; j < n; j++) {
            predicted += 0.5*delta_lm[j] * (mu[j] * delta_lm[j]-g[j]);
         }

         if (Math.abs(predicted) > 0) {
            ratio = actual / predicted;
         } else {
            ratio = 0;
         }

         //check convergence
         if (ss_lm == 0 || (Math.abs(ss0 - ss_lm) / ss_lm < epsilon2 && Math.abs(ss_gn - ss_lm) / ss_lm < epsilon2)) {
            lmresult.iteration = iter+1;
            lmresult.parameters = p_lm;
            lmresult.residualsumsqueres = ss_lm;
            terminateouter = true;
            convergence = true;
            break;
         }

         muMax = 0;

         //successful iteration
         if (ratio > 0) {
            p = p_lm;
            terminateinner = true;
            iter++;
            cntouter++;
            for (j = 0; j < n; j++) {
               mu[j] = mu[j] * Math.max(1/3, 1 - Math.pow(2 * ratio - 1,3));
            }
            v = 2;
         //unsuccessful
         } else {
            for (j = 0; j < n; j++) {
               mu[j] = mu[j] * v;
               muMax = Math.max(muMax,mu[j]);
            }
            v *= 2;
            cntinner++;
         }

         if (!isFinite(muMax) || cntouter > cntoutermax) {
            for (j = 0; j < n; j++) {
               mu[j] = 0;
            }
            v = 2;
            lambda *= lv;
            if (lambda < lambdaMin) {
               lv = 100;
               lambda = lambda0 * lv;
            } else if (lambda > lambdaMax) {
               terminateouter = true;
               break;
            }
            terminateinner = true;
            p = p0.concat();
            cntouter = 0;
         }

         if (cntinner > cntinnermax) {
            return false;
         }

      }

      if (iter > itermax) {
         terminateouter=true;
      }
   }

   return convergence;

}

function  InverseMatrix(matrix) {

   var idm = new Array();
   var mm = new Array();
   var invertm = new Array();

   for (var i = 0; i < matrix.length; i++) {
      idm[i] = new Array();
      for (var j = 0; j < matrix[i].length; j++) {
         if (i == j) {
            idm[i].push(1);
         } else {
            idm[i].push(0);
         }
      }
      mm[i] = matrix[i].concat(idm[i]);
   }

   for (var i = 0; i < mm.length; i++) {
      pivotMatrix(i,mm);
      GaussJordanElimination(i,mm);
   }

   for (var i = 0; i < mm.length; i++) {
      invertm[i] = mm[i].slice(matrix.length);
   }

   return invertm;

}

function pivotMatrix(lNo,matrix) {

   var vMax;
   var curlNo=lNo;
   var curV;
   var i,j;

   vMax = Math.abs(matrix[lNo][lNo]);
   for (i = lNo+1; i < matrix.length; i++) {
      if (vMax < Math.abs(matrix[i][lNo])) {
         curlNo = i;
         vMax = Math.abs(matrix[i][lNo]);
      }
   }

   if (curlNo > lNo) {
      for (i = 0; i < matrix[lNo].length; i++) {
         curV = matrix[lNo][i];
         matrix[lNo][i] = matrix[curlNo][i];
         matrix[curlNo][i] = curV;
      }
   }

}

function GaussJordanElimination(lNo,matrix){

   var v,vv;

   v = matrix[lNo][lNo];

   for (var j = 0; j < matrix[lNo].length; j++) {
      matrix[lNo][j] = matrix[lNo][j] / v;
   }

   for (var i = 0; i < matrix.length; i++) {
      if (i != lNo) {
         vv = matrix[i][lNo];
         for (var j = lNo; j < matrix[lNo].length; j++) {
            matrix[i][j] = matrix[i][j] - vv * matrix[lNo][j];
         }
      }
   }
}

 GNUライセンスですので、ご自由に利用・加工していただいて結構ですが、数学の素人が作成しているので「ある一定の条件では正しく計算されない」といった不具合があるかもしれませんのでご了承下さい。不具合報告歓迎します。
 
 上記プログラムを利用した最適化計算(レーベンバーグ・マーカート法)サンプルのページも作成しました。
 ・Nonlinear Regression (Levenberg–Marquardt algorithm)

 今回は完全に逸脱してますが、本来のプロジェクトは「D3.jsで空間分析にチャレンジ」ですので、データ点のプロット及び求めた関数の描画も出来るようにしました。ソースダウンロード(ZIP圧縮)はこちら

コメント   

 前回試したカーネル密度推定は、確率密度を求めるもので値そのものを補間計算で求めるものではありません(密度を実数値に換算出来るのかもしれませんが)。今回は実数値の推定を試みます。
 サンプル点の値から、平面上のデータを持たない点の値を補間計算する方法のうち、最もシンプルなものの1つがIDW(Inverse Distance Weighting:逆距離加重)法です。ある地点の値を推定する際には当然ながら周辺のサンプル点の値を参考にするわけですが、距離が近いサンプル点ほど相関が強いと考えて、距離に応じた加重平均を行い計算する方法がIDWです。計算式は以下のようになります。
  IDW 数式
  μ(s):地点sの予測値  n:サンプル点数  μ(i):i番目のサンプル点siの実測値   
  d(s,si):地点sとサンプル点siの距離  wi(s):地点sにおけるサンプル点siの重み  p:距離指数

   
 カーネル密度推定では何だか良く分からずに数式をそのままプログラムにしただけ・・・といった感じでしたが、今回のIDWは加重計算ですから私にも式の意味が理解出来ます。
 上記式のうち、p(距離指数)の値を大きくすると、より近い点の影響が強くなります。p=2とすることが多いようです。

 この計算を実行するプログラム(関数)は以下のようにしました。

function IDW(samplepoint, calpoint, p, filtermode, filterval) {
/* 引数
   samplepoint:サンプル点の情報 配列収納 [経度(x),緯度(y),実測値,0] ←4番目の要素(0)は未指定、Null等でも可
   calpoint:数値予測対象点の座標 配列収納 [経度(x),緯度(y)
   p:距離指数(省略時は2.0)
   filtermode:サンプル点のフィルタリング(省略時は全点を対象) 
     "PointCounts":より近い点からn点を対象とする "RadiusRange":半径r(m)以内にあるサンプル点を対象とする
   filterval:フィルタリング条件 
          filtermode="PointCounts"の時:対象点数  filtermode="PointCounts"の時:対象半径(m)
*/
   var idw = [];
   var i,j;
   var ivd;
   var ivdsum,calvalue;
   var x,y;
   var n = samplepoint.length;
   var r = Infinity;

   if (filtermode == "PointCounts" && !isNaN(filterval)) {
      n = filterval;
   }
   if (filtermode == "RadiusRange" && !isNaN(filterval)) {
      r = filterval;
   }
   if(!p || isNaN(p)) {
      p = 2.0;
   }

   for (i in calpoint) {
      ivdsum = 0;
      calvalue = 0;
      for (j in samplepoint) {
         samplepoint[j][3] = getDistance(calpoint[i],samplepoint[j]);
      }
      samplepoint.sort(function(a, b) {
         return d3.ascending(a[3],b[3]);
      })
      for (j in samplepoint) {
         if (j == 0 || (j < n && samplepoint[j][3] < r)) {
            ivd = 1 / Math.pow(samplepoint[j][3],p);
            calvalue += ivd * samplepoint[j][2];
            ivdsum += ivd;
         }
      }
      calvalue = calvalue / ivdsum;
      x = mapPath.projection()([calpoint[i][0],calpoint[i][1]])[0];
      y = mapPath.projection()([calpoint[i][0],calpoint[i][1]])[1];
      idw.push([x,y,calvalue])
   }
   return idw;
}

 全てのサンプル点の影響を考慮して計算する以外に、対象点の限定も良く行われているようなので、
  (1)より近いものからn点を選択
  (2)半径r以内にある点を選択
 の2種類の方法でフィルタリング出来るようにしました。

 全体のソースは以下の通り

・TokyoAirPollutionMapIDW.html

<!DOCTYPE html>
<html>
<head>
<title>東京都の大気汚染(SPM)濃度分布</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="js/d3.js"></script>
<script type="text/javascript" src="js/queue.v1.js"></script>
<script type="text/javascript" src="js/TokyoAirPollutionMapIDW.js"></script>
</style>
<script type="text/javascript">
  function switch_display(mode) {
     if (mode == "ALL") {
        document.getElementById('blankcell').style.display = "table-cell";
        document.getElementById('pcounts').style.display = "none";
        document.getElementById('rradius').style.display = "none";
     } else if (mode == "PointCounts") {
        document.getElementById('blankcell').style.display = "none";
        document.getElementById('pcounts').style.display = "table-cell";
        document.getElementById('rradius').style.display = "none";
     } else if (mode == "RadiusRange") {
        document.getElementById('blankcell').style.display = "none";
        document.getElementById('pcounts').style.display = "none";
        document.getElementById('rradius').style.display = "table-cell";
     }
  }
</script>
</head>
<body style="margin:0px;" onload="setMap(document.getElementById('f_mode').value,document.getElementById('f_pcounts').value,document.getElementById('f_radius').value,document.getElementById('f_power').value,true)">
  <div style="background-color:#006666;padding:5px;">
  <table style="padding-left:5px;">
  <tr>
    <td style="font-size: 16px;font-weight: bold;color: #FFFFFF;">東京都の大気汚染(SPM)濃度分布 </td>
    <td style="font-size: 12px;color: #FFFFFF;">※IDW法にて予測 </td>
    <td style="font-size:12px;color:#FFFFFF;padding-left:10px;">計算対象とする測定局:
      <select id="f_mode" name="f_mode" style="width:140px;font-size:12px;" onchange="switch_display(this.value);">
         <option value="ALL" selected>全て</option>
         <option value="PointCounts">近い点から指定数分</option>
         <option value="RadiusRange">指定の半径内</option>
      </select>
    </td>
    <td id="blankcell" style="padding-left:15px;width:150px;display:table-cell;">&nbsp;</td>
    <td id="pcounts" style="font-size:12px;color:#FFFFFF;padding-left:15px;width:150px;display:none;">対象測定局数:
      <select id="f_pcounts" name="f_pcounts" style="width:50px;font-size:12px;">
         <option value="1">1</option>
         <option value="2">2</option>
         <option value="3">3</option>
         <option value="4">4</option>
         <option value="5">5</option>
         <option value="7">7</option>
         <option value="10" selected>10</option>
         <option value="15">15</option>
         <option value="20">20</option>
      </select>
    </td>
    <td id="rradius" style="font-size:12px;color:#FFFFFF;padding-left:15px;width:150px;display:none;">対象半径:
      <select id="f_radius" name="f_radius" style="width:70px;font-size:12px;">
         <option value="3000">3km</option>
         <option value="5000">5km</option>
         <option value="10000" selected>10km</option>
         <option value="15000">15km</option>
         <option value="20000">20km</option>
         <option value="30000">30km</option>
      </select>
    </td>
    <td id="rradius" style="font-size:12px;color:#FFFFFF;width:110px;">距離指数:
      <select id="f_power" name="f_power" style="width:50px;font-size:12px;">
         <option value="0.5">0.5</option>
         <option value="1.0">1</option>
         <option value="1.5">1.5</option>
         <option value="2.0" selected>2</option>
         <option value="2.5">2.5</option>
         <option value="3.0">3</option>
         <option value="4.0">4</option>
         <option value="5.0">5</option>
      </select>
    </td>
    <td style="padding-left:10px;">
       <input type="button" style="width:80px;font-size:14px;" onclick="setMap(document.getElementById('f_mode').value,document.getElementById('f_pcounts').value,document.getElementById('f_radius').value,document.getElementById('f_power').value,false)" value="更新">
    </td>
  </tr>
  </table>
  </div>
</body>
</html>

・TokyoAirPollutionMapIDW.js

// *************************************************************** 
//  D3.jsで空間分析にチャレンジ
// *************************************************************** 

  var svg;
  var canvas;
  var mapScale = 50000;  //マップスケール
  var mapWidth = 1000;  //マップ幅
  var mapHeight = 500; //マップ高さ
  var mapTop = 40;
  var mapPath, mapBounds;
  var meshCountX=240,meshCountY=120;
  var Mode,pCounts,rRange,dPower;
  var border,outline,point;
  var color = d3.interpolateHsl("#66ff00", "#ff0000");
  var alpha = 0.7;
  var pMin, pMax;  //濃度 Min,Max
  var jsonCityFileName = "json/CityBorderTokyo.json?=" + new Date().getTime();
  var jsonOutlineFileName = "json/OutlineTokyo.json?=" + new Date().getTime();
  var pointFileName = "csv/T_Place_AirPollution.csv?" + new Date().getTime();
  var prefBorder = new Array();

  function setMap(mode,pcnt,radius,power,initFlag) {
     Mode = mode;
     pCounts = pcnt;
     rRange = radius;
     dPower = power;
     if (initFlag) {
        queue()
           .defer(d3.json, jsonCityFileName)
           .defer(d3.json, jsonOutlineFileName)
           .defer(d3.csv, pointFileName)
           .await(ready);
     } else {
        drawCanvas(false);
     }
  }

  function ready(error, data1, data2, data3) {
     border = data1;
     outline = data2;
     point = data3;
     drawMap();
  }

  function drawMap() {

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

     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")

     mapBounds = d3.geo.bounds(border);
     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.forEach(function(d,i) {	
        prefBorder[i] = new Array();
        d["geometry"]["coordinates"][0].forEach(function(dd,j) {
           prefBorder[i].push(mapPath.projection()([parseFloat(dd[0]),parseFloat(dd[1])]));
        });
     });

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

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

     var ObservationPoint = svg.selectAll(".point")
        .data(point)
        .enter()
        .append("circle")
        .attr("cx", function(d,i) {
           return mapPath.projection()([d.Lng,d.Lat])[0];
        })
        .attr("cy", function(d,i) {
           return mapPath.projection()([d.Lng,d.Lat])[1];
        })
        .attr("r", 3)
        .style("fill", "#0000ff")
        .on('mouseover', function(d) {
           d3.select(this).style("fill","#ff0000");
           curText = d.PointName + " " + parseFloat(d.SPM2p).toFixed(3);
           SetTooltip(mapPath.projection()([d.Lng,d.Lat]),curText);
        })
        .on('mouseout', function() {
           d3.select(this).style("fill","#0000ff");
           var tooltip = svg.select(".tooltip")
           if(!tooltip.empty()) {
              tooltip.style("visibility", "hidden")
           }
        })

     drawCanvas(true);

  }

  function drawCanvas(initFlag) {

     var currentcolor;
     var canvas = document.getElementById('canvas')
     var ctx = canvas.getContext('2d');
     var meshSize = 10;
     var xRange = d3.extent([mapPath.projection()(mapBounds[0])[0],mapPath.projection()(mapBounds[1])[0]]);
     var yRange = d3.extent([mapPath.projection()(mapBounds[0])[1],mapPath.projection()(mapBounds[1])[1]]);
     var meshWidth = (xRange[1]-xRange[0])/(meshCountX- 1);
     var meshHeight = (yRange[1]-yRange[0])/(meshCountY- 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));
     var samplepoint = new Array();
     var calpoint = new Array();

     point.forEach(function(d) {
        samplepoint.push([parseFloat(d.Lng),parseFloat(d.Lat),parseFloat(d.SPM2p),0])
     });

     gx.forEach(function(dx) {
        gy.forEach(function(dy) {
           calpoint.push(mapPath.projection().invert([dx,dy]));
        });
     });

     var val;
     if (Mode == "PointCounts") {
        val = pCounts;
     } else if (Mode == "RadiusRange") {
        val = rRange;
     }
     var idw = IDW(samplepoint,calpoint,dPower,Mode,val);
     var vrange = d3.extent(idw, function (d) { return d[2]; });
     pMin = Math.floor(vrange[0]*200)/200;
     pMax = Math.ceil(vrange[1]*200)/200;

     if (initFlag) {
        ctx.beginPath();
        var coordinate; 
        outline.features.forEach(function(d,i) {
           for (var i = 0; i < d["geometry"]["coordinates"][0].length; i++) {
              coordinate = mapPath.projection()([parseFloat(d["geometry"]["coordinates"][0][i][0]),parseFloat(d["geometry"]["coordinates"][0][i][1])])
              if (i==0) {
                 ctx.moveTo(coordinate[0], coordinate[1]);
              } else {
                 ctx.lineTo(coordinate[0], coordinate[1]);
              }
           }
        });
        ctx.closePath();
        ctx.clip();
     }
     ctx.clearRect(0,0,mapWidth,mapHeight); 
     ctx.globalAlpha = alpha;
     idw.forEach(function(d) {
        currentcolor = color((d[2]-pMin)/(pMax-pMin));
        ctx.fillStyle = currentcolor;
        ctx.fillRect(d[0]-0.5*meshWidth,d[1]-0.5*meshHeight,meshWidth,meshHeight);
     });

     setMapLegend(initFlag);
  }

  function IDW(samplepoint, calpoint, p, filtermode, filterval) {
  /* 引数
     samplepoint:サンプル点の情報 配列収納 [経度(x),緯度(y),実測値,0] ←4番目の要素(0)は未指定、Null等でも可
     calpoint:数値予測対象点の座標 配列収納 [経度(x),緯度(y)
     p:距離指数(省略時は2.0)
     filtermode:サンプル点のフィルタリング(省略時は全点を対象) 
      "PointCounts":より近い点からn点を対象とする "RadiusRange":半径r(m)以内にあるサンプル点を対象とする
     filterval:フィルタリング条件 
            filtermode="PointCounts"の時:対象点数  filtermode="PointCounts"の時:対象半径(m)
  */
     var idw = [];
     var i,j;
     var ivd;
     var ivdsum,calvalue;
     var x,y;
     var n = samplepoint.length;
     var r = Infinity;

     if (filtermode == "PointCounts" && !isNaN(filterval)) {
        n = filterval;
     }
     if (filtermode == "RadiusRange" && !isNaN(filterval)) {
        r = filterval;
     }
     if(!p || isNaN(p)) {
        p = 2.0;
     }

     for (i in calpoint) {
        ivdsum = 0;
        calvalue = 0;
        for (j in samplepoint) {
           samplepoint[j][3] = getDistance(calpoint[i],samplepoint[j]);
        }
        samplepoint.sort(function(a, b) {
           return d3.ascending(a[3],b[3]);
        })
        for (j in samplepoint) {
           if (j == 0 || (j < n && samplepoint[j][3] < r)) {
              ivd = 1 / Math.pow(samplepoint[j][3],p);
              calvalue += ivd * samplepoint[j][2];
              ivdsum += ivd;
           }
        }
        calvalue = calvalue / ivdsum;
        x = mapPath.projection()([calpoint[i][0],calpoint[i][1]])[0];
        y = mapPath.projection()([calpoint[i][0],calpoint[i][1]])[1];
        idw.push([x,y,calvalue])
     }
     return idw;
  }

  //度数→ラジアン変換
  function calRadian(deg){
     return deg * Math.PI / 180.0;
  }

  //2点間の距離を計算(GRS80)
  function getDistance(p1, p2) {
     var a = 6378137.0;             //長半径
     var f = 1.0 / 298.257222101;   //扁平率
     var b = a * (1 - f);           //短半径
     var e2 = (Math.pow(a, 2) - Math.pow(b, 2)) / Math.pow(a, 2);
     var dy = calRadian(p1[1] - p2[1]);
     var dx = calRadian(p1[0] - p2[0]);
     var my = calRadian((p1[1] + p2[1]) / 2.0);
     var w = Math.sqrt(1.0 - e2 * Math.pow(Math.sin(my), 2));
     var m = (a * (1 - e2)) / Math.pow(w, 3);
     var n = a / w;

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

  function SetTooltip(tipXY,tipText) {

     var tooltip = svg.select(".tooltip")
     var fontSize = 12;
     var rectWidth = tipText.length * fontSize + 5;
     var rectHeight = 20;

     if(tooltip.empty()) {
        tooltip = svg
           .append("g")
           .attr("class","tooltip")

        tooltip
           .append("rect")
           .attr("height",rectHeight)
           .style("stroke","none")
           .style("fill","#ffffff")
           .style("opacity","0.8")

        tooltip
           .append("text")
           .attr("text-anchor","left")
           .attr("id","tooltiptext")
           .style("font-size",fontSize+"px")
           .style("font-family","sans-serif")
     }

     var tipleft = tipXY[0] + 2;
     if (tipleft + rectWidth > mapWidth) {
        tipleft = mapWidth - rectWidth;
     }
     var tiptop = tipXY[1] - rectHeight - 2;

     tooltip
        .style("visibility", "visible")
        .attr("transform", "translate("+tipleft+","+tiptop+")")

     tooltip.select("text")
        .text(tipText + "mg/m")
        .attr("transform", "translate(5,"+(fontSize+(rectHeight-fontSize)/2)+")")
        .append("tspan")
        .attr("id","super")
        .attr("dy","-2")
        .style("font-size",(fontSize-2)+"px")
        .text("3")

     tooltip.select("rect")
        .attr("width",rectWidth)

  }

  function setMapLegend(initFlag) { 

    var fontsize = 14;
     var legendWidth = 200;
     var legendHeight = 20;
     var legendLeft = 30;
     var legendTop = mapPath.projection()(mapBounds[0])[1] - 60;
     var valuerange = d3.range(pMin,pMax+0.0001,(pMax-pMin)/(legendWidth-1));
     var pFormat = d3.format(".3f");

     if (!initFlag) {
        var v = [pMin,pMax];
        d3.selectAll(".axistext")
         .text(function(d,i) {
              return pFormat(v[i]) +"mg/m";
           })
           .append("tspan")
           .attr("class","super")
           .attr("dy","-2")
           .style("font-size",(fontsize-4)+"px")
           .text("3")
        return false;
     }

     var MapLegend = d3.select("svg")
        .append("g")
        .attr("id","maplegend")
        .attr("transform", "translate("+legendLeft+","+legendTop+")")

     MapLegend
        .append("circle")
        .attr("r", 3)
        .style("fill", "#0000ff")
        .on('mouseover', function(d) {
           d3.select(this).style("fill","#ff0000");
        })
        .on('mouseout', function() {
           d3.select(this).style("fill","#0000ff");
        })

     MapLegend
        .append("text")
        .attr("id","text1")
        .style("font-size",fontsize + "px")
        .style("font-family","sans-serif")
        .style("fill","#000000")
        .attr("transform", "translate("+10+","+ (fontsize * 0.5 - 2) +")")
      .text("大気汚染常時監視測定局(一般局)")

     MapLegend
        .append("text")
        .attr("id","text2")
        .style("font-size",(fontsize - 2) + "px")
        .style("font-family","sans-serif")
        .style("fill","#000000")
        .attr("transform", "translate("+10+","+ (fontsize * 1.5) +")")
      .text("※マウスオンすると測定局の名称とSPM濃度(2013年度2%除外値)が表示されます。")

     MapLegend.selectAll(".legend")
        .data(valuerange)
        .enter()
        .append("rect")
        .attr("class", "legend")
        .attr("transform", function(d, i) {
           return "translate("+((fontsize*10)+i)+",30)";
        })
        .attr("width",1)
        .attr("height",legendHeight)
        .style("fill", function(d, i){
           return color((d - pMin) / (pMax - pMin));
        })

     MapLegend
        .append("rect")
        .attr("id","legendframe")
        .attr("transform", "translate("+(fontsize*10)+",30)")
        .attr("width",legendWidth)
        .attr("height",legendHeight)
        .style("fill", "none")
        .style("stroke","#303030")
        .style("stroke-width",0.5);

     MapLegend
        .append("text")
        .attr("id","text3")
        .style("font-size",fontsize + "px")
        .style("font-family","sans-serif")
        .style("fill","#000000")
        .attr("transform", "translate(0,"+ (30 + (fontsize + legendHeight) * 0.5) +")")
      .text("SPM濃度(2%除外値)")

     var MapLegendAxis = MapLegend   //.selectAll(".axis")
        .append("g")
        .attr("transform", "translate("+(fontsize*10)+","+ (30 + legendHeight) +")")
        .attr("class","axis")
        .style("stroke","none")
        .style("fill","#303030")

     MapLegendAxis
        .append("rect")
        .attr("transform", "translate(0,10)")
        .attr("width",legendWidth)
        .attr("height",1)

     var MapLegendAxisX = MapLegendAxis.selectAll(".axis")
        .data([d3.min(valuerange),d3.max(valuerange)])
        .enter()
        .append("g")

     MapLegendAxisX
        .append("rect")
        .attr("transform", function(d,i) {
           return "translate(" + (i * legendWidth) + ",5)";
        })
        .attr("width",1)
        .attr("height",10)

     MapLegendAxisX
        .append("text")
        .attr("transform", function(d,i) {
           return "translate("+(i * legendWidth)+","+(15+(fontsize+2))+")";
        })
        .attr("text-anchor","middle")
        .attr("class","axistext")
      .attr("id",function(d,i) {
           return "axistext"+i;
        })
        .style("font-size",(fontsize-2) + "px")
        .style("font-family","sans-serif")
      .text(function(d,i) {
           return pFormat(d) +"mg/m";
        })

     MapLegendAxisX.selectAll("text")
        .append("tspan")
        .attr("class","super")
        .attr("dy","-2")
        .style("font-size",(fontsize-4)+"px")
        .text("3")

  }

 緯度・経度情報から平面座標に変換するに当たっては、投影法としてメルカトル図法を採用しています。
  ※以下のd3.geo.mercator()の箇所

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

 メルカトル図法では高緯度ほど歪みが大きくなるので、地図上の1目盛(1ピクセル)当たりの距離は高緯度になるほど短くなります。今回の地図描画の対象は、緯度範囲が小さい東京都(島しょ部を除く)なので誤差は無視出来るレベルですが、日本全図を描画する際等は誤差が大きくなり補間計算にも影響が出てきそうです。そのため、IDW法に用いる距離は平面座標(x,y)の値を緯度・経度に戻した上で求めることにしました(以下の部分)。

           calpoint.push(mapPath.projection().invert([dx,dy]));

 d3.jsには平面座標を緯度・経度に戻してくれる便利なメソッド「projection.invert」が用意されています。

 補間計算の結果はフィルタリング条件やp(距離指数)の設定によって異なりますので、実際には出力される地図を見ながら調整することになるかと思います(例を下図に示します)。全サンプル点を対象としp=2とした場合は濃度の高いサンプル点近隣が目立つ出力となりました(図1)。今回扱っているデータはサンプル点=観測点であって発生源ではないので、サンプル点周辺のみにピークがあるというのは実態とは合っていない気がします。いろいろと条件を変更して試した所、対象半径5km以内でp=1とした結果では幹線道路(環八・環七・1号線等)周辺の濃度が高いような印象で、比較的実態に合っているように思います(図2)。もう少し観測点の数が多ければもっと良い結果が得られるのでしょうが・・・。
 ちなみに、最近隣の1点のみを対象とした場合はボロノイ図となります(図3)。240×120のメッシュではギザギザした出力となりますが、メッシュサイズを小さくすれば前々回作成したボロノイ図と同じ形状になるかと思います。

図 東京都の大気汚染(SPM)濃度分布(IDW法による メッシュ数:240×120)

サンプル点のフィルタリングなし(全点対象) p=2

図1 サンプル点のフィルタリングなし(全点対象) p=2


サンプル点のフィルタリング:半径5km以内 p=1

図2 サンプル点のフィルタリング:半径5km以内 p=1


サンプル点のフィルタリング:最近隣1点のみを対象=ボロノイ図

図3 サンプル点のフィルタリング:最近隣1点のみを対象 ※ボロノイ図となります

 
 

コメント   

 前回はD3.jsのメソッドを利用してボロノイ図を作成しましたが、今の所D3.jsには他の空間分析系の関数は見当たりませんので、自ら作成するしかありません。

 基本的に空間分析は「R」言語とセットになっているようで、これがないと何も出来ないような状況ですので、まずはRをインストール。
 以下ページに丁寧な説明があり、特に問題なくインストール出来ました。 
   R とパッケージの簡単インストール

 同時にRと空間分析の書籍を購入。
   地理空間データ分析(Rで学ぶデータサイエンス 7)共立出版

 まだまだ学習中で全然読み進んでいないのですが、本書で紹介されている「カーネル密度推定法」というものが比較的導入し易そうなので、今回はこれを試すことにしました。カーネル密度推定法は犯罪発生マップ等に利用されているようです。
 利用例:大阪府警察犯罪発生マップ
       ※カーネル密度推定は「車上ねらい等」のマップに適用されています。
 
 カーネル密度推定とはイベント発生点の分布状況から全体の分布状況を推定する手法の1つで、次の計算式によって求めます。
  カーネル密度推定数式
  n:標本点数  Xi(i=1,2,…,n):各点の値(座標) x:密度推定点の値(座標)
  h:バンド幅(標本点からの広がりの範囲)  K:カーネル関数
  
 このうちKのカーネル関数に適用される関数にはガウス関数、イパネクニコフ関数、四次関数など数種類あり、どれを採用すれば良いのか判断が難しいですが、前記の書籍によれば「採用するカーネル関数による影響よりも、バンド幅(h)の選択による影響の方が大きい」とありますので、今回はガウス関数の1つである次の標準正規分布を採用することにしました。
  カーネル関数(標準正規分布)
 この式を上のカーネル密度推定数式に当てはめると、
  カーネル密度推定数式(ガウス関数を適用)
 となります。
 これをjavascript+D3.jsに実装。標本としたデータは、アメダス(東京都府中市)の日平均気温(2013年)です。

<!DOCTYPE html>
<meta charset="utf-8">
<title>カーネル密度推定 アメダスデータ(東京都府中市 日平均気温 2013年)</title>
<style>

   .axis text {
      font-family: sans-serif;
      font-size: 11px;
   }
   .axis path,
   .axis line {
      fill: none;
      stroke: #000000;
      stroke-width: 1px;
      shape-rendering: crispEdges;
   }
   .legend text {
      font-family: sans-serif;
      font-size: 11px;
   }

</style>
<script type="text/javascript" src="js/d3.js"></script>
<script type="text/javascript" src="js/queue.v1.js"></script>
<script>

   var svg;
   var width = 800;
   var height = 400;
   var marginLeft = 100;
   var marginTop = 20;
   var xInterval = 1;
   var xRange = [-5,35];

   var xScale = d3.scale.linear()
      .domain([xRange[0], xRange[1]])
      .range([0, width])

   var yScale = d3.scale.linear()
      .domain([0, 0.07])
      .range([height, 0])

   var xAxis = d3.svg.axis()
      .scale(xScale)
      .orient("bottom")
      .ticks(10)

   var yAxis = d3.svg.axis()
      .scale(yScale)
      .orient("left")
      .ticks(10)

   var line = d3.svg.line()
      .x(function(d) { return marginLeft + xScale(d[0]); })
      .y(function(d) { return marginTop + yScale(d[1]); })

   var histogram = d3.layout.histogram()
      .frequency(false)
      .bins(d3.range(xRange[0],xRange[1]+0.001,xInterval));

   function initialize() {

      svg = d3.select("body").append("svg")
         .attr("id", "svg")
         .attr("width", width+200)
         .attr("height", height+150)

      queue()
         .defer(d3.csv, "csv/temperature_fuchu2013.csv?=" + new Date().getTime())
         .await(ready);
   }

   function ready(error, data) {

      var x = new Array();
      data.forEach(function(d) {
         x.push(parseFloat(d.Temperature));
      });

      var hist = histogram(x);

      svg.selectAll(".hist")
         .data(hist)
         .enter()
         .append("rect")
         .attr("class", "hist")
         .attr("x", function(d) { return marginLeft + xScale(d.x); })
         .attr("y", function(d) { return marginTop + yScale(d.y); })
         .attr("width", width/(xRange[1]-xRange[0])*xInterval)
         .attr("height", function(d) { return height - yScale(d.y); })
         .attr("fill", "#E0E0E0")
         .attr("stroke", "#909090")
         .attr("stroke-width", "0.5px")

      var gx = d3.range(xRange[0],xRange[1]+0.001,xInterval/10);

      var kde = new Array();
      var bandwidth = [0.5, 2.0, 5.0, 10.0];
      bandwidth.forEach(function(d) {
         kde.push(GaussianKernelDensityEstimation(x,gx,d));
      });
      var lineColor = ["#0000C0","#00C000","#C00000","#F0F000"];

      var lines = svg.append("g")
         .attr("class","lines")
         .attr("fill", "none")
         .attr("stroke-width", "2px")

      lines.selectAll("path")
         .data(kde)
         .enter()
         .append("path")
         .attr("id", function(d,i) { return "kde"+i; })
         .attr("d", function(d,i) { return line(d); })
         .attr("stroke", function(d,i) { return lineColor[i]; })

      drawAxis();
      drawLegend(bandwidth,lineColor);

   }

   function drawAxis(){

      d3.select("#svg")
         .append("g")
         .attr("class", "axis")
         .attr("transform", "translate("+marginLeft+", "+marginTop+")")
         .call(yAxis)
         .append("text")
            .attr("class", "label")
            .attr("x", -50)
            .attr("y", height/2)
            .attr("fill", "#000000")
            .style("writing-mode", "tb")
            .style("text-anchor", "middle")
            .text("確率密度");

      d3.select("#svg")
         .append("g")
         .attr("class", "axis")
         .attr("transform", "translate("+marginLeft+", "+(marginTop+height)+")")
         .call(xAxis)
         .append("text")
            .attr("class", "label")
            .attr("x", width/2)
            .attr("y", 35)
            .attr("fill", "#000000")
            .style("text-anchor", "middle")
            .text("日平均気温(℃)");

   }

   function drawLegend(bw,cl){

      var offsetY = 40;

      var legendline = svg.append("g")
             .attr("fill", "none")
             .attr("stroke-width", "2px")
             .attr("class", "legend")

      legendline.selectAll("path")
         .data(bw)
         .enter()
         .append("path")
         .attr("d", function(d,i) { 
            return "M"+ (marginLeft+width-100)+","+(marginTop+height+offsetY+i*15)+"L"+ (marginLeft+width-50)+","+(marginTop+height+offsetY+i*15);
         })
         .attr("stroke", function(d,i) { return cl[i]; })

      var format = d3.format(".1f");

      legendline.selectAll("text")
         .data(bw)
         .enter()
         .append("text")
         .attr("x", marginLeft+width-25)
         .attr("y", function(d,i) { 
            return marginTop+height+offsetY+4+i*15;
         })
         .style("text-anchor", "end")
         .attr("fill", "#000000")
         .text(function(d,i) { return format(d); })

      legendline.append("text")
         .attr("x", marginLeft+width-115)
         .attr("y", function(d,i) { 
            return marginTop+height+offsetY+4;
         })
         .style("text-anchor", "end")
         .attr("fill", "#000000")
         .text("バンド幅")

      svg.append("text")
         .attr("x", marginLeft+50)
         .attr("y", marginTop+20)
         .style("font-size","12px")
         .text("東京都府中市 2013年 日平均気温の分布状況  ※出典:気象庁 アメダスデータ")

   }

   function GaussianKernelDensityEstimation(x, gx, h) {

      var kde = [];
      var dd;

      if(typeof h === 'undefined') {
         h = SilvermanBandWidth(x);
      }

      for (i in gx) {
         dd = 0;
         for (ii in x) {
            dd += (Math.exp(-0.5*Math.pow((gx[i]-x[ii])/h,2)))/Math.sqrt(2 * Math.PI);
         }
         kde.push([gx[i],dd/(x.length*h)])
      }
      return kde;

   }

   function SilvermanBandWidth(data,mode) {   
      var n = 0;
      var Sxx = 0;
      var Ex,Sd;
      var d = [].concat(data);
      var IQR = d3.quantile(d.sort(d3.ascending), 0.75) - d3.quantile(d.sort(d3.ascending), 0.25);

      Ex = d3.mean(d);
      for (var i in d) {
         Sxx += Math.pow(Ex - d[i],2);
      }
      Sd = Math.sqrt(Sxx/(d.length-1));

      if (mode == 0) {
         return 0.9 * d3.min([IQR/1.34,Sd]) * Math.pow(d.length,-0.2);  //Silverman's rule
      } else {
         return 1.06 * d3.min([IQR/1.34,Sd]) * Math.pow(d.length,-0.2); //Scott's rule
      }
   }

</script>
<body onload="initialize()">

 4通りのバンド幅(0.5,2.0,5.0,10.0)にてそれぞれ計算しています。また、ヒストグラムも作成しています。
 これをブラウザで表示させると、下図のようなグラフとなります。

東京都府中市の日平均気温分布(2013年)

東京都府中市の日平均気温分布(2013年)
このグラフを表示


 バンド幅の設定によって密度分布の形が大きく変わることがよく分かります。また、標本データの存在しない所でも、その周辺に標本点があれば密度推定値はゼロにはならないことも確認出来ます(上図の31~32℃の箇所)。

 と、ここまでは1次元の場合。今回は地図上に表現したいので2次元で計算する必要があるのですが、その計算式が分からない・・・Rや統計学の参考書を調べてみても、何の説明もなく用意された関数を使用しているものがほとんど。たまに説明されているものを見つけても専門的過ぎる上に計算例が記載されていないので理解出来ません。仕方がないのでRの2次元カーネル密度推定関数の構造を解析しました。
 Rには無数のライブラリ(パッケージと呼ぶようです)があり、2次元カーネル密度推定関数もいくつかあるのですが、その中でも最もシンプルな式を採用していると思われるkde2d(MASSパッケージに収納)を対象としました。当関数では次のような計算をしています。
  カーネル密度推定数式(2次元)
  n:標本点数  Xi,Yi(i=1,2,…,n):各点の値(X座標,Y座標)
  x,y:密度推定点の値(X座標,Y座標) hx,hy:バンド幅(X座標方向,Y座標方向)

 この式を見ると、X座標の値とY座標の値は互いに独立しているとみなし、共分散は考慮していないようです。それでいいのかどうかはよく分からないのですが、とりあえずはこのままの式で進めます。また、関数kde2dでは引数として与えたバンド幅の値を計算時に1/4にしているのですが(2次元の平面だから22ということでしょうか?)、今回は引数で与えた値をそのまま計算式に組み込むことにします。

 以下、2次元カーネル密度推定の実装プログラムソースです。 

function GaussianKernelDensityEstimation2D(x, y, nx, ny, h, rangex, rangey) {
/* 引数
     x,y:イベント発生点の座標(X座標,y座標) ※両変数とも配列形式でx,yの要素数は同数
   nx,ny:密度推定点数(X座標方向、Y座標方向)※計算対象の点数はnx×nyとなります。  
   h:バンド幅
    ※要素数=2の配列形式 [X座標方向のバンド幅,Y座標方向のバンド幅] 
    ※配列でなく数値を渡した場合は、その値がX,Y両方向共通のバンド幅として適用される
    ※Nullを渡した場合は、シルバーマンの方法(イベント発生点の分散や四分位範囲などから計算)により設定する
     rangex,rangey:密度推定点を配置する範囲(X座標方向、Y座標方向)要素数=2の配列形式
     例: rangex=[0,50]・・・Min=0,Max=50の範囲で配置
       nx=6であれば、推定点のX座標は0,10,20,30,40,50となる     
*/
   var hx,hy;
   var rx0,rx1,ry0,ry1;

   if(!h || isNaN(h) || ((h instanceof Array) && (h.length != 2))) {
      hx = SilvermanBandWidth(x);
      hy = SilvermanBandWidth(y);
   } else if (!(h instanceof Array)) {
      hx = h;
      hy = h;
   } else {
      hx = h[0];
      hy = h[1];
   }

   if((typeof rangex === 'undefined') || !(rangex instanceof Array) || (rangex.length != 2)) {
      rx0 = d3.min(x) - hx*1.5;
      rx1 = d3.max(x) + hx*1.5;
   } else {
      rx0 = rangex[0];
      rx1 = rangex[1];
   }
   if((typeof rangey === 'undefined') || !(rangey instanceof Array) || (rangey.length != 2)) {
      ry0 = d3.min(y) - hy*1.5;
      ry1 = d3.max(y) + hy*1.5;
   } else {
      ry0 = rangey[0];
      ry1 = rangey[1];
   }

   var gx = d3.range(rx0,rx1+Math.pow(10,-10),(rx1-rx0)/(nx-1));
   var gy = d3.range(ry0,ry1+Math.pow(10,-10),(ry1-ry0)/(ny-1));

   var kde = [];
   var dd;
   var i,j,ii;

   for (i in gx) {
      for (j in gy) {
         dd = 0;
         for (ii in x) {
            dd += (Math.exp(-0.5*Math.pow((gx[i]-x[ii])/hx,2))*Math.exp(-0.5*Math.pow((gy[j]-y[ii])/hy,2)))/(2 * Math.PI);
         }
         kde.push([gx[i],gy[j],(dd / (x.length * hx * hy))])
      }
   }
   return kde;
}

  
 前述したバンド幅の調整の部分を除けばkde2dと同じ結果が得られる・・・と思います。

 この計算を実際の地図上に反映させるに当たり、前回ボロノイ図作成時に使用した東京都大気汚染常時監視測定局は、「観測点」であって実際に大気汚染物質を排出している場所ではないので、カーネル密度推定に使用する点としては適当ではありません。そのため、今回は代わりに東京都内の清掃工場を使用することにしました。実際には大気汚染物質というのは自動車からの排出が半分以上で、しかも最近の清掃工場は環境にも配慮しているでしょうから都内の排出総量から考えれば微々たるものなのですが、一応排出源ではありますので標本点としてふさわしいかと思います。
 密度推定点は地図上の平面に格子状に配置します(メッシュ)。そして前回のボロノイ図と同様に、この計算結果の数値に応じて各メッシュを着色・・・すれば良いのですが、試したところメッシュを細かくすると描画に異常に時間がかかり固まってしまいます。svgのpath生成処理に時間がかかるのかとおもったのですが、色をつけなければ速いので描画自体に時間がかかる模様。ブラウザ(IE,Chrome,FireFox,Opera)による差異も大きく、速い順にChrome>Opera>IE>FireFoxとなりました。特にFireFoxの遅さが突出しています。
 描画速度はブラウザのバージョンにより改善される可能性はありますが、現時点ではメッシュ数=20,000(200×100)位にすると最も速いChromeでも相当時間がかかるので、実用に耐えられません。となると他の方法を考えなければなりませんが、他といってもCanvasくらいしか思いつきません。メッシュのみをCanvasに描画しSVGと重ねる方法を試してみました。以下が、そのソースです。

・WasteIncinerationPlantTokyo.html

<!DOCTYPE html>
<html>
<head>
<title>東京都の清掃工場分布とカーネル密度推定</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="js/d3.js"></script>
<script type="text/javascript" src="js/queue.v1.js"></script>
<script type="text/javascript" src="js/WasteIncinerationPlantTokyo.js"></script>
</style>
</head>

<body style="margin:0px;" onload="setMap(document.getElementById('f_mesh').value,document.getElementById('f_bandwidth').value,true)">
  <div style="background-color:#666600;padding:5px;">
  <table style="padding-left:5px;">
  <tr>
    <td style="font-size: 16px;font-weight: bold;color: #FFFFFF;">東京都の清掃工場分布とカーネル密度推定 </td>
    <td style="font-size:12px;color:#FFFFFF;padding-left:30px;">メッシュ:
      <select id="f_mesh" name="f_mesh" style="width:80px;font-size:12px;">
         <option value="20,10">20×10</option>
         <option value="50,25">50×25</option>
         <option value="100,50">100×50</option>
         <option value="160,80">160×80</option>
         <option value="240,120">240×120</option>
      </select>
    </td>
    <td style="font-size:12px;color:#FFFFFF;padding-left:20px;">バンド幅:
      <select id="f_bandwidth" name="f_bandwidth" style="width:180px;font-size:12px;">
         <option value="10">10</option>
         <option value="20" selected>20</option>
         <option value="30">30</option>
         <option value="40">40</option>
         <option value="50">50</option>
         <option value="">自動(SilvermanBandWidth)</option>
      </select>
    </td>
    <td style="padding-left:10px;">
       <input type="button" style="width:50px;font-size:14px;" onclick="setMap(document.getElementById('f_mesh').value,document.getElementById('f_bandwidth').value)" value="更新">
    </td>
  </tr>
  </table>
  </div>
</body>
</html>

・WasteIncinerationPlantTokyo.js

// *************************************************************** 
//  D3.jsで空間分析にチャレンジ
// *************************************************************** 

  var svg;
  var mapScale = 50000;
  var mapWidth = 1000;
  var mapHeight = 500;
  var mapTop = 40;
  var mapPath, mapBounds;
  var meshCountX,meshCountY,bandWidth;
  var jsonCityFileName = "json/CityBorderTokyo.json?=" + new Date().getTime();
  var jsonOutlineFileName = "json/OutlineTokyo.json?=" + new Date().getTime();
  var pointFileName = "csv/T_WasteIncinerationPlant.csv?" + new Date().getTime();
  var pointX = new Array();
  var pointY = new Array();
  var border,outline,point;
  var color = d3.interpolateHsl("#66ff00", "#ff0000");
  var alpha = 0.7;

  function setMap(meshxy,bw,initFlag) {
     bandWidth = bw;
     meshCount = meshxy.split(",");
     meshCountX = meshCount[0];
     meshCountY = meshCount[1];
     if (initFlag) {
        queue()
           .defer(d3.json, jsonCityFileName)
           .defer(d3.json, jsonOutlineFileName)
           .defer(d3.csv, pointFileName)
           .await(ready);
     } else {
        drawCanvas(false);
     }
  }

  function ready(error, data1, data2, data3) {
     border = data1;
     outline = data2;
     point = data3;
     drawMap();
  }

  function drawMap() {

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

     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")

     mapBounds = d3.geo.bounds(border);
     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)
        );

     point.forEach(function(d) {
        pointX.push(mapPath.projection()([d.Lng,d.Lat])[0]);
        pointY.push(mapPath.projection()([d.Lng,d.Lat])[1]);
     });

     drawCanvas(true);

     svg.selectAll(".outline")
        .data(outline.features)
        .enter()
        .append("path")
        .attr("d", mapPath)
        .style("fill", "none")
        .attr("class", "outline")
        .style("stroke", "#C0C0C0")
        .style("stroke-width", "2")

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

     var WasteIncinerationPlant = svg.selectAll(".plant")
        .data(point)
        .enter()
        .append("circle")
        .attr("cx", function(d,i) {
           return mapPath.projection()([d.Lng,d.Lat])[0];
        })
        .attr("cy", function(d,i) {
           return mapPath.projection()([d.Lng,d.Lat])[1];
        })
        .attr("r", 3)
        .style("fill", "#0000ff")
        .on('mouseover', function(d) {
           d3.select(this).style("fill","#ff0000");
           curText = d.PlaceName;
           SetTooltip(mapPath.projection()([d.Lng,d.Lat]),curText);
        })
        .on('mouseout', function() {
           d3.select(this).style("fill","#0000ff");
           var tooltip = svg.select(".tooltip")
           if(!tooltip.empty()) {
              tooltip.style("visibility", "hidden")
           }
        })

     setMapLegend();

  }

  function drawCanvas(initFlag) {

     var xRange = d3.extent([mapPath.projection()(mapBounds[0])[0],mapPath.projection()(mapBounds[1])[0]]);
     var yRange = d3.extent([mapPath.projection()(mapBounds[0])[1],mapPath.projection()(mapBounds[1])[1]]);
     var meshWidth = (xRange[1]-xRange[0])/(meshCountX- 1);
     var meshHeight = (yRange[1]-yRange[0])/(meshCountY- 1);

     var kde = GaussianKernelDensityEstimation2D(pointX, pointY, meshCountX, meshCountY, bandWidth, xRange, yRange);
     var drange = d3.extent(kde, function (d) { return d[2]; });

     var currentcolor;
     var canvas = document.getElementById('canvas')
     var ctx = canvas.getContext('2d');

     if (initFlag) {
        ctx.beginPath();
        var coordinate; 
        outline.features.forEach(function(d,i) {
           for (var i = 0; i < d["geometry"]["coordinates"][0].length; i++) {
              coordinate = mapPath.projection()([parseFloat(d["geometry"]["coordinates"][0][i][0]),parseFloat(d["geometry"]["coordinates"][0][i][1])])
              if (i==0) {
                 ctx.moveTo(coordinate[0], coordinate[1]);
              } else {
                 ctx.lineTo(coordinate[0], coordinate[1]);
              }
           }
        });
        ctx.closePath();
        ctx.clip();
     }

     ctx.clearRect(0,0,mapWidth,mapHeight); 
     ctx.globalAlpha = alpha;
     kde.forEach(function(d) {
        currentcolor = color((d[2]-drange[0])/(drange[1]-drange[0]));
        ctx.fillStyle = currentcolor;
        ctx.fillRect(d[0]-0.5*meshWidth,d[1]-0.5*meshHeight,meshWidth,meshHeight);
     });

  }

  function GaussianKernelDensityEstimation2D(x, y, nx, ny, h, rangex, rangey) {

     var hx,hy;
     var rx0,rx1,ry0,ry1;

     if(!h || isNaN(h) || ((h instanceof Array) && (h.length != 2))) {
        hx = SilvermanBandWidth(x);
        hy = SilvermanBandWidth(y);
     } else if (!(h instanceof Array)) {
        hx = h;
        hy = h;
     } else {
        hx = h[0];
        hy = h[1];
     }

     if((typeof rangex === 'undefined') || !(rangex instanceof Array) || (rangex.length != 2)) {
        rx0 = d3.min(x) - hx*1.5;
        rx1 = d3.max(x) + hx*1.5;
     } else {
        rx0 = rangex[0];
        rx1 = rangex[1];
     }
     if((typeof rangey === 'undefined') || !(rangey instanceof Array) || (rangey.length != 2)) {
        ry0 = d3.min(y) - hy*1.5;
        ry1 = d3.max(y) + hy*1.5;
     } else {
        ry0 = rangey[0];
        ry1 = rangey[1];
     }

     var gx = d3.range(rx0,rx1+Math.pow(10,-10),(rx1-rx0)/(nx-1));
     var gy = d3.range(ry0,ry1+Math.pow(10,-10),(ry1-ry0)/(ny-1));

     var kde = [];
     var dd;
     var i,j,ii;

     for (i in gx) {
        for (j in gy) {
           dd = 0;
           for (ii in x) {
               dd += (Math.exp(-0.5*Math.pow((gx[i]-x[ii])/hx,2))*Math.exp(-0.5*Math.pow((gy[j]-y[ii])/hy,2)))/(2 * Math.PI);
           }
           kde.push([gx[i],gy[j],(dd / (x.length * hx * hy))])
        }
     }
     return kde;
  }

  function SilvermanBandWidth(data,mode) {

     var n = 0;
     var Sxx = 0;
     var Ex,Sd;
     var d = [].concat(data),i;
     var IQR = d3.quantile(d.sort(d3.ascending), 0.75) - d3.quantile(d.sort(d3.ascending), 0.25);

     Ex = d3.mean(d);

     for (i in d) {
        Sxx += Math.pow(Ex - d[i],2);
     }
     Sd = Math.sqrt(Sxx/(d.length-1));

     if (mode == 0) {
        return 0.9 * d3.min([IQR/1.34,Sd]) * Math.pow(d.length,-0.2);  //Silverman's rule
     } else {
        return 1.06 * d3.min([IQR/1.34,Sd]) * Math.pow(d.length,-0.2); //Scott's rule
     }
  }

  function SetTooltip(tipXY,tipText) {

     var tooltip = svg.select(".tooltip")
     var fontSize = 12;
     var rectWidth = tipText.length * fontSize;
     if (tipText.length < 15) {
        rectWidth += 15 - tipText.length;
     }
     var rectHeight = 20;

     if(tooltip.empty()) {
        tooltip = svg
           .append("g")
           .attr("class","tooltip")

        tooltip
           .append("rect")
           .attr("height",rectHeight)
           .style("stroke","none")
           .style("fill","#ffffff")
           .style("opacity","0.5")

        tooltip
           .append("text")
           .attr("text-anchor","left")
           .attr("id","tooltiptext")
           .style("font-size",fontSize+"px")
           .style("font-family","sans-serif")
     }

     var tipleft = tipXY[0] + 2;
     if (tipleft + rectWidth > mapWidth) {
        tipleft = mapWidth - rectWidth;
     }
     var tiptop = tipXY[1] - rectHeight - 2;

     tooltip
        .style("visibility", "visible")
        .attr("transform", "translate("+tipleft+","+tiptop+")")

     tooltip.select("rect")
        .attr("width",rectWidth)

     tooltip.select("text")
        .text(tipText)
        .attr("transform", "translate(5,"+(fontSize+(rectHeight-fontSize)/2)+")")

  }

  function setMapLegend() { 

    var fontsize = 14;
     var legendLeft = 100;
     var legendTop = mapPath.projection()(mapBounds[0])[1] - 70;
     var legendWidth = 240;
     var legendHeight = 20;
     var rangeWidth = 4;
     var valuerange = d3.range(0,1.00001,1/(legendWidth/rangeWidth-1));

     var MapLegend = d3.select("svg")
        .append("g")
        .attr("id","maplegend")
        .attr("transform", "translate("+legendLeft+","+legendTop+")")

     MapLegend
        .append("circle")
        .attr("r", 3)
        .style("fill", "#0000ff")
        .on('mouseover', function(d) {
           d3.select(this).style("fill","#ff0000");
        })
        .on('mouseout', function() {
           d3.select(this).style("fill","#0000ff");
        })

     MapLegend
        .append("text")
        .attr("id","text1")
        .style("font-size",fontsize + "px")
        .style("font-family","sans-serif")
        .style("fill","#000000")
        .attr("transform", "translate("+10+","+ (fontsize * 0.5 - 2) +")")
      .text("清掃工場")

     MapLegend
        .append("text")
        .attr("id","text2")
        .style("font-size",(fontsize - 2) + "px")
        .style("font-family","sans-serif")
        .style("fill","#000000")
        .attr("transform", "translate("+20+","+ (fontsize * 1.5) +")")
      .text("※マウスオンすると名称が表示されます。")

     MapLegend.selectAll(".legend")
        .data(valuerange)
        .enter()
        .append("rect")
        .attr("class", "legend")
        .attr("transform", function(d,i) {
           return "translate("+(i*rangeWidth)+",32)";
        })
        .attr("width",rangeWidth)
        .attr("height",legendHeight)
        .style("fill", function(d){
           return color(d);
        })
        .style("opacity",alpha)

     MapLegend
        .append("rect")
        .attr("id","legendframe")
        .attr("transform", "translate(0,32)")
        .attr("width",legendWidth)
        .attr("height",legendHeight)
        .style("fill", "none")
        .style("stroke","#303030")
        .style("stroke-width",0.5);

     var MapLegendAxis = MapLegend
        .append("g")
        .attr("transform", "translate(0,"+ (32 + legendHeight + 5) +")")
        .attr("class","axis")
        .style("font-size",fontsize + "px")
        .style("font-family","sans-serif")
        .style("fill","#000000")

     MapLegendAxis
        .append("path")
        .attr("d", "M10,10L"+(legendWidth/2-(fontsize*3+10)/2)+",10M10,10L20,5M10,10L20,15"
                   +"M"+(legendWidth-10)+",10L"+(legendWidth/2+(fontsize*3+10)/2)+",10"
                   +"M"+(legendWidth-10)+",10L"+(legendWidth-20)+",5M"+(legendWidth-10)+",10L"+(legendWidth-20)+",15")
        .style("stroke","#808080")
        .style("stroke-width","0.5")

     MapLegendAxis
        .append("text")
        .attr("transform", "translate("+legendWidth/2+","+(10+fontsize*0.5)+")")
        .attr("text-anchor","middle")
        .attr("class","axistext")
      .text("密 度")

     MapLegendAxis
        .append("text")
        .attr("transform", "translate(0,"+(10+fontsize*0.5)+")")
        .attr("text-anchor","middle")
        .attr("class","axistext")
      .text("低")

     MapLegendAxis
        .append("text")
        .attr("transform", "translate("+legendWidth+","+(10+fontsize*0.5)+")")
        .attr("text-anchor","middle")
        .attr("class","axistext")
      .text("高")

  }

 Canvasへの描画はSVGよりもかなり速いです。これでメッシュを細かくしても(240×120)、さほどストレスなく表示されるようになりました。出力例を下図に示します。

図 東京都の清掃工場分布とカーネル密度推定(メッシュ数:240×120)

東京都の清掃工場分布とカーネル密度推定(バンド幅=10)

バンド幅=10

東京都の清掃工場分布とカーネル密度推定(バンド幅=20)

バンド幅=20

東京都の清掃工場分布とカーネル密度推定(バンド幅=30)

バンド幅=30

 実際には工場の煙の拡散には風向・風速等の気象条件が大きく影響するのですが、今回は考慮していません。前回作成したボロノイ図と比較しても、より「空間分析」っぽくなってきた気がします。
 清掃工場から排出される汚染物質がそれほど広い範囲に影響を及ぼすとは思えないので、上図の中では一番上のバンド幅=10が実態に合っているのではないかと思いますが、バンド幅=20もしくは30の方が見栄えがいいですね。見栄えで判断すべきではないのでしょうが、バンド幅=10だと単なる清掃工場分布図と大して変わらないですね。

コメント   

 空間分析といっても手法は様々であり、取っ掛かりがなかなかつかめません。
 「何から始めたらいいものか・・・」と思っていた所、d3.jsにはボロノイ図を作成するための関数が存在することを知りました。

 ボロノイ図とは「複数あるポイントのうち、最寄りとなるポイントは何処か?」という視点で、平面空間を領域に分けたものです(下図参照)。

ボロノイ図

ボロノイ図


 この例ですと、緑色の領域内ではFが最寄りのポイントとなります。
 新聞、雑誌、Web等の媒体でボロノイ図を見かけることはほとんどないので実際にどのような用途で利用されているのかは分かりませんが、イメージ的にはスーパーやコンビニの出店や携帯のアンテナ設置等の計画を立てる際にツールとして利用出来そうです。

 上図のようなシンプルなボロノイ図を出力するソースを以下に記載します。

<!DOCTYPE HTML>
<html>
<head>
<title>ボロノイ図サンプル</title>
<meta http-equiv="content-type" content="text/html; charset=UTF-8">
<meta http-equiv="Content-Style-Type" content="text/css">
<script type="text/javascript" src="js/d3.js"></script>
</head>

<body style="margin:0px;" onload="draw_voronoi()">
  <div style="padding:10px;">ボロノイ図サンプル</div>
</body>
</html>

<script type="text/javascript">

  var svg;
  var svgWidth =  400;
  var svgHeight = 300;
  var svgPadding = 5;

  function draw_voronoi() {

     svg = d3.select("body").append("svg")
        .attr("width", svgWidth)
        .attr("height", svgHeight)
        .attr("id", "svg")
        .style("margin-left","25px")

     var point = [[60,210],[120,90],[180,120],[240,240],[210,60],[300,150],[360,210]];

     var voronoiData = d3.geom.voronoi()
                           .clipExtent([[svgPadding, svgPadding], [svgWidth - svgPadding, svgHeight - svgPadding]]);

     var voronoiPath = svg.selectAll(".voronoi")
        .data(voronoiData(point))
        .enter()
        .append("path")
        .attr("class", "voronoi")
        .attr("d", function(d, i) { 
             return "M" + d.join("L") + "Z";
        })
        .style("stroke","#000099")
        .style("stroke-width","1")
        .style("fill", "none")

     var pointElm = svg.selectAll(".point")
        .data(point)
        .enter()
        .append("circle")
        .attr("cx", function(d,i) {
           return d[0];
        })
        .attr("cy", function(d,i) {
           return d[1];
        })
        .attr("r", 3)
        .style("fill", "#000000")
  }
</script>

 d3.geom.voronoi関数を使用すれば、ポイント情報をボロノイ図用に加工してくれます。

var voronoiData = d3.geom.voronoi()
                      .clipExtent([[svgPadding, svgPadding], [svgWidth - svgPadding, svgHeight - svgPadding]]);

 clipExtentは、ボロノイの作成範囲を指定するメソッドで必須ではありませんが、本家d3.jsのサイトでは使用を推奨しています。
 そして、

.data(voronoiData(point))

と、データを設定します。すると、各ポイントについてボロノイ領域を構成するポリゴンの座標が配列に収納されますので

.attr("d", function(d, i) {
   return "M" + d.join("L") + "Z";
})

 といった表記で、配列の各要素(座標)をLでつなげることにより、
  M140,5L165,75L245,115L355,5Z
 このような形式の、SVGのpath要素にポリゴン描画するための文字列を生成します。

 これらを応用して前回作成した東京都大気汚染常時監視測定局の地図をボロノイ図にしてみました。

東京都大気汚染常時監視測定局 ボロノイ図(1)

東京都大気汚染常時監視測定局 ボロノイ図(1)
このMAPを表示

 一応出来たのですが、見栄えが良くないです。ボロノイの線は東京都の輪郭内に限定出来ないものか・・・。
 いろいろ探して見つけたのが以下のサイト
  Areas of Clipped Voronoi Regions
 「Areas of Clipped~」とのタイトルの通り、ボロノイ図の描画エリアは地図の輪郭内に収まっています。当ページに公開されているソースを解読して、
  d3.geom.polygon().clip
 このメソッドで、クリップ(トリミングと言った方が分かりやすいでしょうか?)していることが分かりました。
 早速これを試してみたのですが・・・どうも当メソッドでクリップ出来る図形は1つだけのようでマルチポリゴンには対応していないようで(実際には対応しているのかもしれませんが実装方法が分かりません)、次のようなボロノイ図となりました。

東京都大気汚染常時監視測定局 ボロノイ図(2)

東京都大気汚染常時監視測定局 ボロノイ図(2)
このMAPを表示


 東京都のメインの輪郭は着色されていますが、中央防波堤付近の島は白のままとなっています(赤い丸の箇所)。この程度であれば無視してしまっても良いのですが、これでは佐渡ヶ島のある新潟県や、島だらけの長崎県などには適用出来そうにありません。また、上の図を良く見ると着色に関しては「クリップ」されているのですが、ボロノイの線に関しては所々でクリップエリアをはみ出してしまっています(青い丸の箇所)。

 他の方法を探った所、d3.js固有ではないのですがSVGにはクリップパス(clipPath)という切り抜き用の要素があることを知りました。
 まずはdefsタグ(表示対象ではない定義用の要素)を用意し、中にclipPathタグを入れてそこに切り抜きの「型」を記述します。そして、対象の図形にclip-path=”url(#クリップパスのID)”と付加することで切り抜きが実行されます。
 クリップパス実装のシンプルなソースを以下に記載します。

<!DOCTYPE HTML>
<html>
<head>
<title>SVG ClipPath Sample</title>
<meta http-equiv="content-type" content="text/html; charset=UTF-8">
<meta http-equiv="Content-Style-Type" content="text/css">
<script type="text/javascript" src="js/d3.js"></script>
</head>

<body onload="draw_clip()">
   <div style="padding:10px;">ClipPath Sample</div>
</body>
</html>

<script type="text/javascript">

  var svg;
  var svgWidth =  500;
  var svgHeight = 300;

  function draw_clip() {

     svg = d3.select("body").append("svg")
        .attr("width", svgWidth)
        .attr("height", svgHeight)
        .attr("id", "svg")
        .style("margin-left","10px")
        .style("background-color","#808080")

     var defs = svg.append("defs");
     var clip = defs.append("clipPath")
        .attr("id", "clip-two-circle")

     clip.append("circle")
        .attr("cx", 125)
        .attr("cy", 150)
        .attr("r", 100)
     clip.append("circle")
        .attr("cx", 375)
        .attr("cy", 150)
        .attr("r", 100)

     svg.append("rect")
        .attr("clip-path", "url(#clip-two-circle)")
        .attr("x", 0)
        .attr("y", 0)
        .attr("width", svgWidth)
        .attr("height", svgHeight)
        .style("fill", "#0000ff")
     svg.append("rect")
        .attr("clip-path", "url(#clip-two-circle)")
           .attr("x", 100)
           .attr("y", 120)
           .attr("width", 300)
           .attr("height", 60)
           .style("fill", "#ff0000")
  }
</script>

 まずは、

     var defs = svg.append("defs");
     var clip = defs.append("clipPath")
        .attr("id", "clip-two-circle")

 とsvgにdefs要素、さらにその中にclipPath要素を生成し、clipPathのIDをclip-two-circleとしました。
 このclipPath要素に

     clip.append("circle")
        .attr("cx", 125)
        .attr("cy", 150)
        .attr("r", 100)
     clip.append("circle")
        .attr("cx", 375)
        .attr("cy", 150)
        .attr("r", 100)

 と円を2つ横に並べていますが、これらは全てdegs要素内で定義しているにすぎないので、これだけでは何も表示されません。
 実際に表示させているのは、

     svg.append("rect")
        .attr("clip-path", "url(#clip-two-circle)")
        .attr("x", 0)
        .attr("y", 0)
        .attr("width", svgWidth)
        .attr("height", svgHeight)
        .style("fill", "#0000ff")
     svg.append("rect")
        .attr("clip-path", "url(#clip-two-circle)")
           .attr("x", 100)
           .attr("y", 120)
           .attr("width", 300)
           .attr("height", 60)
           .style("fill", "#ff0000")

 この部分。SVG全体と同サイズ(幅・高さとも)の大きな四角形(青)と300×60の小さな四角形(赤)を描画し、両図形とも
.attr(“clip-path”, “url(#clip-two-circle)”)
 とすることで、clipPathで定義した形に切り抜いています。結果として出力される図形は下画像の通りです。

SVG ClipPath利用の図形例

SVG ClipPath利用の図形例
この図を表示

 このclipPathを先の東京都大気汚染常時監視測定局のボロノイ図に実装してみました。

・ObservationPointsTokyoVoronoi.html

<!DOCTYPE html>
<html>
<head>
<title>東京都大気汚染常時監視測定局 ボロノイ図</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="js/d3.js"></script>
<script type="text/javascript" src="js/queue.v1.js"></script>
<script type="text/javascript" src="js/ObservationPointsTokyoVoronoi.js"></script>
</style>
</head>

<body style="margin:0px;" onload="initialize()">
  <div  style="background-color:#006666;padding: 10px;font-size: 16px;font-weight: bold;color: #FFFFFF;">東京都大気汚染常時監視測定局 ボロノイ図</div>
</body>
</html>

・ObservationPointsTokyoVoronoi.js

// *************************************************************** 
//  D3.jsで空間分析にチャレンジ
// *************************************************************** 

  var svg;
  var mapScale = 50000;  //マップスケール
  var mapWidth = 1000;  //マップ幅
  var mapHeight = 500; //マップ高さ
  var mapPath;
  var vcolor = d3.interpolateHsl("#66ff00", "#ff0000");
  var pMin = 0.045;  //濃度 Min
  var pMax = 0.080;  //濃度 Max
  var jsonCityFileName = "json/CityBorderTokyo.json?=" + new Date().getTime();
  var jsonOutlineFileName = "json/OutlineTokyo.json?=" + new Date().getTime();
  var pointFileName = "csv/T_Place_AirPollution.csv?" + new Date().getTime();
  var prefBorder = new Array();

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

  function ready(error, border, outline, observepoint) {

     svg = d3.select("body").append("svg")
        .attr("width", mapWidth)
        .attr("height", mapHeight)
        .attr("id", "svg")

     var mapbounds = d3.geo.bounds(border);
     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.forEach(function(d,i) {	
        prefBorder[i] = new Array();
        d["geometry"]["coordinates"][0].forEach(function(dd,j) {
           prefBorder[i].push(mapPath.projection()([parseFloat(dd[0]),parseFloat(dd[1])]));
        });
     });

     setVoronoi(observepoint);

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

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

     setMapLegend(mapbounds);

  }

  function setVoronoi(point) {

     var ObservePoint = svg.selectAll(".observepoint")
        .data(point)
        .enter()
        .append("g")
        .attr("class","observepoint")
   
     var strClip = "";
     for (var p = 0; p < prefBorder.length; p++) {
        strClip += mkPolygonPath(prefBorder[p]);
     }

     var defs = svg.append("defs");
     var clipBox = defs.append("clipPath")
        .attr("id", "clip")
        .append("path")
        .attr("d",function() { return strClip;})

     var voronoiData = d3.geom.voronoi(pickupCoordinates(point));

     var voronoiPath = ObservePoint.append("path")
        .attr("class", "voronoi")
        .attr("d", function(d, i) { 
           return(mkPolygonPath(voronoiData[i])); 
        })
        .attr("clip-path", "url(#clip)")
        .style("stroke","#000099")
        .style("stroke-width","0.5")
        .style("fill", function(d, i) { 
             var fcolor = vcolor((parseFloat(point[i].SPM2p) - pMin) / (pMax - pMin));
             return(fcolor)
        })

     var pointElm = ObservePoint.append("circle")
        .attr("cx", function(d,i) {
           return mapPath.projection()([d.Lng,d.Lat])[0];
        })
        .attr("cy", function(d,i) {
           return mapPath.projection()([d.Lng,d.Lat])[1];
        })
        .attr("r", 3)
        .style("fill", "#000099")
        .on('mouseover', function(d) {
           d3.select(this).style("fill","#ff0000");
           curText = d.PointName + " " + parseFloat(d.SPM2p).toFixed(3);
           SetTooltip(mapPath.projection()([d.Lng,d.Lat]),curText,"3");
        })
        .on('mouseout', function() {
           d3.select(this).style("fill","#000099");
           var tooltip = svg.select(".tooltip")
           if(!tooltip.empty()) {
              tooltip.style("visibility", "hidden")
           }
        })
  }

  function pickupCoordinates(point) {
     var coordinatesArray = new Array;
     for (var i = 0; i < point.length; i++) {
        coordinatesArray.push(mapPath.projection()([point[i].Lng,point[i].Lat]));
     }
     return coordinatesArray;
  }

  function mkPolygonPath(d) {
     return "M" + d.join("L") + "Z";
  }

  function SetTooltip(tipXY,tipText) {

     var tooltip = svg.select(".tooltip")
     var fontSize = 12;
     var rectWidth = tipText.length * fontSize + 5;
     var rectHeight = 20;

     if(tooltip.empty()) {
        tooltip = svg
           .append("g")
           .attr("class","tooltip")

        tooltip
           .append("rect")
           .attr("height",rectHeight)
           .style("stroke","none")
           .style("fill","#ffffff")
           .style("opacity","0.8")

        tooltip
           .append("text")
           .attr("text-anchor","left")
           .attr("id","tooltiptext")
           .style("font-size",fontSize+"px")
           .style("font-family","sans-serif")
     }

     var tipleft = tipXY[0] + 5;
     if (tipleft + rectWidth > mapWidth) {
        tipleft = mapWidth - rectWidth;
     }
     var tiptop = tipXY[1] - rectHeight -5;

     tooltip
        .style("visibility", "visible")
        .attr("transform", "translate("+tipleft+","+tiptop+")")

     tooltip.select("text")
        .text(tipText + "mg/m")
        .attr("transform", "translate(5,"+(fontSize+(rectHeight-fontSize)/2)+")")
        .append("tspan")
        .attr("id","super")
        .attr("dy","-2")
        .style("font-size",(fontSize-2)+"px")
        .text("3")

     tooltip.select("rect")
        .attr("width",rectWidth)

  }

  function setMapLegend(bounds) { 

    var fontsize = 14;
     var legendWidth = 200;
     var legendHeight = 20;
     var legendLeft = 30;
     var legendTop = mapPath.projection()(bounds[0])[1] - 60;
     var valuerange = d3.range(pMin,pMax+0.0001,(pMax-pMin)/(legendWidth-1));

     var MapLegend = d3.select("svg")
        .append("g")
        .attr("id","maplegend")
        .attr("transform", "translate("+legendLeft+","+legendTop+")")

     MapLegend
        .append("circle")
        .attr("r", 3)
        .style("fill", "#0000ff")
        .on('mouseover', function(d) {
           d3.select(this).style("fill","#ff0000");
        })
        .on('mouseout', function() {
           d3.select(this).style("fill","#0000ff");
        })

     MapLegend
        .append("text")
        .attr("id","text1")
        .style("font-size",fontsize + "px")
        .style("font-family","sans-serif")
        .style("fill","#000000")
        .attr("transform", "translate("+10+","+ (fontsize * 0.5 - 2) +")")
      .text("大気汚染常時監視測定局(一般局)")

     MapLegend
        .append("text")
        .attr("id","text2")
        .style("font-size",(fontsize - 2) + "px")
        .style("font-family","sans-serif")
        .style("fill","#000000")
        .attr("transform", "translate("+10+","+ (fontsize * 1.5) +")")
      .text("※マウスオンすると測定局の名称とSPM濃度(2013年度2%除外値)が表示されます。")

     MapLegend.selectAll(".legend")
        .data(valuerange)
        .enter()
        .append("rect")
        .attr("class", "legend")
        .attr("transform", function(d, i) {
           return "translate("+((fontsize*10)+i)+",30)";
        })
        .attr("width",1)
        .attr("height",legendHeight)
        .style("fill", function(d, i){
           return vcolor((d - pMin) / (pMax - pMin));
        })

     MapLegend
        .append("rect")
        .attr("id","legendframe")
        .attr("transform", "translate("+(fontsize*10)+",30)")
        .attr("width",legendWidth)
        .attr("height",legendHeight)
        .style("fill", "none")
        .style("stroke","#303030")
        .style("stroke-width",0.5);

     MapLegend
        .append("text")
        .attr("id","text3")
        .style("font-size",fontsize + "px")
        .style("font-family","sans-serif")
        .style("fill","#000000")
        .attr("transform", "translate(0,"+ (30 + (fontsize + legendHeight) * 0.5) +")")
      .text("SPM濃度(2%除外値)")

     var MapLegendAxis = MapLegend
        .append("g")
        .attr("transform", "translate("+(fontsize*10)+","+ (30 + legendHeight) +")")
        .attr("class","axis")
        .style("stroke","none")
        .style("fill","#303030")

     MapLegendAxis
        .append("rect")
        .attr("transform", "translate(0,10)")
        .attr("width",legendWidth)
        .attr("height",1)

     var MapLegendAxisX = MapLegendAxis.selectAll(".axis")
        .data([d3.min(valuerange),d3.max(valuerange)])
        .enter()
        .append("g")

     MapLegendAxisX
        .append("rect")
        .attr("transform", function(d,i) {
           return "translate(" + (i * legendWidth) + ",5)";
        })
        .attr("width",1)
        .attr("height",10)

     var pFormat = d3.format(".3f");

     MapLegendAxisX
        .append("text")
        .attr("transform", function(d,i) {
           return "translate("+(i * legendWidth)+","+(15+(fontsize+2))+")";
        })
        .attr("text-anchor","middle")
        .attr("class","axistext")
      .attr("id",function(d,i) {
           return "axistext"+i;
        })
        .style("font-size",(fontsize-2) + "px")
        .style("font-family","sans-serif")
      .text(function(d,i) {
           return pFormat(d) +"mg/m";
        })

     MapLegendAxisX.selectAll("text")
        .append("tspan")
        .attr("class","super")
        .attr("dy","-2")
        .style("font-size",(fontsize-4)+"px")
        .text("3")

  }

 クリップパス以外では、上で紹介したサイト「Areas of Clipped Voronoi Regions」でボロノイ領域の面積に応じて色を変化させているのを参考に、各測定局の濃度によってボロノイ領域の色を緑(低濃度)→赤(高濃度)へと変化させてみました。
 また、上記ソースではD3.js以外にライブラリ「Queue.js」を利用しています。これはcsvやjson等の外部ファイルの読み込みが完了するまで処理待ちをしてくれるものです。
  Queue.jsダウンロードページ(GitHub):https://github.com/mbostock/queue 

 出来上がったのが次の図です。

東京都大気汚染常時監視測定局 ボロノイ図(3)

東京都大気汚染常時監視測定局 ボロノイ図(3)
このMAPを表示


 ようやく「らしい」図になりました。
 
 しかし、ここまではSVGやd3.jsのメソッドを利用しているだけで、統計処理の計算は全く行っていません。次回はもう少し本格的な計算をしてみたいと思います。

コメント