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

クリギング計算式

クリギング計算式


 上式の重み係数(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という名称で保存・使用しています) 
 

コメント