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

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

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


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

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


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

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

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

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


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

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

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

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

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

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


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

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

・MinamiKantoAirPollutionMakeCalpoint.html

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

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

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

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

・MinamiKantoAirPollutionMakeCalpoint.js

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

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

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

  function setCalPoint() {

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

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

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

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

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

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

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

     saveCalPoint(cpselect);  //CSV保存

  }

  function drawMap() {

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

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

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

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

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

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

  }

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

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

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

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

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

・savecalpoint.php

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

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

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

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

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

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

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

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

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

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

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

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

  

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

セミバリオグラム定義式


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

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


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

図3 セミバリオグラム


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

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


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

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

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

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

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

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

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

    各モデル式のパラメータ

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

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

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

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

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


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

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

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

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

 熊ノ平小屋は、最寄りの山である間ノ岳からも結構離れているので(山の地図によれば登りは3時間)、ご来光を見るにはロケーション的に難しいのですが、テラスからの朝風景はとても綺麗でした。夜の星も素晴らしく、本当に心地良い所でした。

 この日は熊ノ平小屋から北岳山頂を経由して広河原へ降りるルートです。前日よりも楽な工程ですし、広河原からのバスは伊那市バスとは違い夕方まで便数豊富なので特に急ぐ必要はないのですが、日曜日の夕方となると広河原は大混雑必至。小屋出発は早いに越したことはないので、夜明けとともに出て朝食は途中で、と当初は計画していたのですが、前日に朝食時間を伺うと4:30からと早め・・・それならばいただいてから出発しよう!と予定変更しました。朝食も夕食同様にとても美味しかったので正解でした。

 5時少し前に出発。間ノ岳まではずっと登りですが険しい部分はほとんどなく、あっさりと山頂に到着。昨年に標高が改定されて1m高くなり、日本3位タイに昇格した間ノ岳ですが、山標は変更されていないので以前の標高のまま。あくまでも個人的な印象ですが、相変わらずパッとしない山ですねえ(^^;)
 それでも北岳-間ノ岳ルートのルートに入ると、すれ違う登山者の数がずっと多くなりました。装備を見ると北岳山荘辺りに荷物を預けて間ノ岳までピストンしている人が結構多そうな感じです。間ノ岳ってそこまでして登る山ではない気がするんですが、やはり百名山になっている関係でしょうかねえ。なぜ深田久弥さんは間ノ岳を百名山に選んだのでしょう・・・選定ポイントは標高だけ?
 と、随分悪く書いてしまいましたが、間ノ岳は北岳を眺めるには最高のポイントだと思います。北岳はこちら側からのアングルが一番カッコよく見えます。
 
 今回のルートでは北岳山頂までの登りも結構楽でした。広河原から登る際の一番の難所は梯子が連続する八本歯のコル周辺なんでしょうが、難所といっても大したことはないので、北岳も塩見岳同様に難易度低めです(もちろん、北岳バットレスなどのバリエーションルートは別)。
 当然といえば当然ですが、北岳山頂は間ノ岳以上の賑わいでした。天気も良く眺望抜群。気持ち的なものですが、ここからの富士山の眺めは格別です! 富士山を見るには日本で一番高い場所なのですから・・・。

 広河原までの下りルートは、4年前に来た時には台風の影響で閉鎖されていた大樺沢コースを選択したのですが、当コースは渡渉がいくつもあるほか、道自体が沢になりかけている箇所も結構ありました。ゴアテックスなど防水性の高い登山靴でないと厳しいかも・・・特に登りで靴を濡らしてしまうと後が大変です。私は、今回は軽さ重視で普通のジョギングシューズ・・・奇跡的にほとんど濡らさずに走破出来ましたがコース選択を間違えました。

 毎日欠かさずジョギングをしていても、やはり登山では使う筋肉が違うので全身筋肉痛。いい運動したなあと実感します。

熊ノ平小屋の朝風景

熊ノ平小屋の朝風景

朝日を浴びる塩見岳

朝日を浴びる塩見岳

間ノ岳山頂付近

間ノ岳山頂付近

間ノ岳山頂

間ノ岳山頂。山標は旧標高のまま

間ノ岳より望む北岳

間ノ岳より望む北岳

北岳山頂

北岳山頂

賑わう北岳山頂

やはり北岳山頂は賑わっていました。

北岳山頂より富士山を望む

北岳山頂より富士山を望む

北岳山頂からの眺望 鳳凰三山

北岳山頂からの眺望
鳳凰三山

北岳山頂からの眺望 八ヶ岳

北岳山頂からの眺望
八ヶ岳

北岳山頂からの眺望 塩見岳

北岳山頂からの眺望
塩見岳

北岳山頂からの眺望 仙丈ヶ岳

北岳山頂からの眺望
仙丈ヶ岳

北岳山頂からの眺望 甲斐駒ヶ岳

北岳山頂からの眺望
甲斐駒ヶ岳

北岳バットレス

北岳バットレス
 

 週末(7月25・26日)を利用して南アルプスの塩見岳・北岳を縦走してきました。
 塩見岳はずっと登りたいと思っていたのですが、登山口までのバスの運行期間が短い(今年は7/18~8/30)ことなどから、なかなかタイミングが合わず、今回ようやく実現出来ました。北岳は4年前の秋以来2回目です。
 
 鳥倉登山口行きバスの発着点である伊那大島駅までは、夜行の毎日あるぺん号を利用しました。毎日あるぺん号は各路線とも竹橋始発ですが、駒ヶ根・伊那大島・戸台口行きは八王子から乗車出来るので府中在住の私には便利です。中央道のルートでは駒ヶ根よりも伊那大島(松川インター)の方が遠いのに、当バスは最初に伊那大島駅に着きます(戸台口が中央道から離れている関係ですかね)。八王子発が24:15なのに3時過ぎには到着しました。結構近いですね。登山バスの出発時刻は6:45と時間が空いてしまうので、「同じ塩見岳行きの人が何人かいればタクシーに相乗りさせてもらおうかな」なんて淡い期待を持っていたのですが、伊那大島で降りたのは私一人・・・塩見岳人気ないんですかね? さすがに一人でタクシーに乗るのは気が引けるので(料金12,000円ですし、そもそも予約してません)、朝まで駅で仮眠・・・伊那大島駅は拍子抜けする程小さな駅で駅前には店の1件もありませんが、夜中でも開いていたので助かりました。

 鳥倉登山口までのバスの料金は2,490円。「確かネットで調べた時には1,600円くらいだったはず・・・記憶違いかな?」と思ったのですが、この日山小屋で一緒になった方から聞いた話では荷物代が50%加算されるとのこと。なるほど、大切な荷物なので子供料金相当は仕方なしということなのでしょう。まあ、伊那市も運営が大変なことでしょうしバス路線を維持してくれるだけでも感謝しないと・・・。
 バスは「よくこんな急坂を登れるなあ」なんて思うような道を通り、登山口到着は8時半くらい。塩見岳・北岳を2日間で縦走するとなると、宿泊地は熊ノ平小屋のほぼ一択となりますが、鳥倉登山口から熊ノ平小屋までの所要時間は地図によれば13時間弱。私の場合は機動力重視なので、そこまではかからないでしょうが9時間くらいは見込んでおいた方が安全。その場合、8時半スタートだと17:30着、暗くなる前には着けそうですが、これだと夕食を用意してもらえない可能性があるので、予め食料は多めに準備しておきました。
 熊ノ平小屋ではちょうど逆ルートで縦走されている方とお会いしましたが、伊那大島駅・松川インター行きのバスは鳥倉発14:25なので、そちらのルートの方が時間的にタイトですね(その方は広河原から熊ノ平小屋まで7時間で走破されたそうなので、きっと間に合ったことでしょう)。

 塩見岳山頂までの登山道は整備されていてとても歩きやすかったです。三伏峠小屋までは10分の1毎に案内板があり参考になります。山頂近くは岩がゴツゴツしていますが、鎖やロープが必要となるほどではなく、子供でも登れそうな感じです。ただ、今年は塩見小屋が建替工事で休業中なので、ゆっくり登山される方には工程的に少しキツイですかねえ。
 一方、山頂から熊ノ平小屋までのルートは登山者も多くないのか、倒木や生い茂った草で道が分かり辛くなっている箇所もありました。といっても、コースロスしてしまう程には荒れていませんから、悪天候でなければ危険度は低いです。

 そんな歩き易い道にも助けられ所要時間は6時間半少々。計画よりも随分早い15時過ぎに到着出来ました。熊ノ平小屋はトイレが離れた所にあるなど設備が最小限でこじんまりとした造りですが、とても居心地の良い所です。美味しいと評判の食事も噂に違わず素晴らしい! 過去に泊まった山小屋の中で最高かも。
 でも、この時期の土曜日にも関わらず宿泊客は15人程度とガラガラ。元々、メジャーなルート上に位置していないので、というのもあるのでしょうが、今年の場合は塩見小屋が利用出来ないことが影響しているかもしれません。三伏峠小屋まで(から)では遠すぎるとの理由で、このルートを避けてしまう人も結構いるのではないかと思います。
 お奨め出来る素晴らしい小屋なので、利用客が少なくて閉鎖なんてことにならないようもっと賑わって欲しいのですが、静岡市営なのでそんな心配は不要かな?
 しかし、静岡県の上に突き出た部分の先っぽに位置するこんな場所(下図参照)まで静岡市とは驚きです! さすがは静岡市、一時期は面積日本一でしたからねえ。
 ※平成の大合併の関係で、現在は5位まで下がってしまいました。

静岡県

熊ノ平小屋の位置。こんな所でも静岡市です!


 
 

伊那大島駅

伊那大島駅。小さな駅ですが夜中でも開いてました

鳥倉登山口

鳥倉登山口

塩見岳登山道

登山道は整備されていて歩きやすいです

三伏峠小屋までの案内板

三伏峠小屋までは10分の1毎に案内板があります

塩見小屋は建替工事中

塩見小屋は建替工事中。2016夏に再開予定
 

塩見岳

塩見岳。深田久弥が表現した「鉄兜」の言葉がまさにピッタリ。勇ましい山容です

塩見岳山頂付近

塩見岳山頂付近。ゴツゴツしてますが難易度は低め

塩見岳山頂(西峰)

塩見岳山頂(西峰)

塩見岳山頂(東峰)

塩見岳山頂(東峰)
こちらの方が5m程高いです

塩見岳山頂からの眺望(南東方面)富士山

塩見岳山頂からの眺望(南東方面)
富士山

塩見岳山頂からの眺望(南方面)

塩見岳山頂からの眺望(南方面)
荒川岳など

塩見岳山頂からの眺望(西方面)中央アルプス

塩見岳山頂からの眺望(西方面)
中央アルプス

塩見岳山頂からの眺望(北方面)仙丈ヶ岳・甲斐駒ヶ岳・北岳・間ノ岳など

塩見岳山頂からの眺望(北方面)
仙丈ヶ岳・甲斐駒ヶ岳・北岳・間ノ岳など

塩見岳山頂より甲斐駒ヶ岳を望む

塩見岳山頂より甲斐駒ヶ岳を望む
甲斐駒は夏でも白いので遠くからでもすぐに分かります

ちょっとしたお花畑

ちょっとしたお花畑

塩見岳(竜尾見晴付近より)

塩見岳(竜尾見晴付近より)

熊ノ平小屋

熊ノ平小屋

熊ノ平小屋の目の前にそびえる農鳥岳

熊ノ平小屋の目の前にそびえる農鳥岳

 これまで、「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圧縮)はこちら

府中市の郷土の森では、現在あじさいまつりが開催されています。
今がちょうど見頃。期間は7月5日(日)までですが、陽当たりの良い所の紫陽花の中には枯れかかっているものもあったので、7月を迎える頃には既に終わってしまっているのではないかと思います。

あじさいまつり 2015

あじさいまつり 2015

あじさいまつり 2015

あじさいまつり 2015

あじさいまつり 2015

あじさいまつり 2015

あじさいまつり 2015

あじさいまつり 2015

あじさいまつり 2015

あじさいまつり 2015

あじさいまつり 2015

あじさいまつり 2015

春と秋の2回、限られた期間のみ入場出来る生田緑地ばら苑。今春の公開期間は5月14日~5月31日です。

生田緑地ばら苑

生田緑地ばら苑

今月は夏日の暑い日が続いているせいか、もう終わりかけている品種も結構ありました。香りはあまり・・・午後に行ったせいかもしれません。香りを楽しむには早朝が良いそうですが、ここは9時開場(平日は10時)のようです。
まあ、姿を見るだけでも十分楽しめました。

インカ

インカ

エスメラルダ

エスメラルダ

キャラメルアンティーク

キャラメルアンティーク

ブランデンブルグ

ブランデンブルグ

アメリカンヘリテイジ

アメリカンヘリテイジ

エリナ

エリナ

トラディション

トラディション

ストロベリー・ダイキリ

ストロベリー・ダイキリ

グラハム・トーマス

グラハム・トーマス

マチルダ

マチルダ

ブルー・フォー・ユー

ブルー・フォー・ユー

ガーデンオブローゼズ

ガーデンオブローゼズ

サマースノー

サマースノー

クンバヤ

クンバヤ

宮城野

宮城野

ユメオトメ

ユメオトメ

しかし、ここが向ケ丘遊園の跡地だったとは・・・初めて知りました。

5月17日(土)
世田谷美術館にて開催中の企画展「速水御舟とその周辺 大正期日本画の俊英たち」を観覧してきました。

「速水御舟とその周辺 大正期日本画の俊英たち」リーフレット


御舟作品だけでも企画が成り立つくらいの点数なのに、本展では今村紫紅や松本瓢湖など御舟と縁が深い画家の作品も出品されているので、なんとも贅沢な構成になっています。
中でも、御舟とそのライバル小茂田青樹両名の「猫」「朝顔」「月」など同じ題材の作品を、まるで競作であるかのように並べて展示する試みが一番楽しめました。

本展は前期と後期とで大幅な展示替えが行われます。前期の目玉は「洛北修学院村」でしょうか・・・素晴らしかったです。リーフレットの表紙に採用されている「菊火図」は後期展示なので、もう一度行かなければ。

前期は5月31日(日)まで。後期は6月2日(火)~7月5日(日)

洛北修学院村

洛北修学院村
≪速水御舟≫


四季草花図(夏季)

四季草花図(夏季)
≪小茂田青樹≫

5月5日(火・祝)
練馬区立美術館へ初訪問。「没後100年 小林清親展 文明開化の光と影をみつめて」を観覧してきました。

練馬区立美術館

練馬区立美術館

小林清親は「最後の浮世絵師」などと呼ばれているようですが、その肩書きは確か川瀬巴水にも使われているような・・・。東京の風景を描いた作品が多いので、描く対象も日本橋、増上寺、亀戸の藤など川瀬巴水と共通するものが結構あり、巴水好きの私としては「ああ、これはやっぱり巴水の方がいいかなあ・・・」なんてついつい思ってしまうのですが、もちろん清親の作品も素晴らしいです。
風景画の他、動物を描いた作品も秀逸です(点数は多くないですが)。風刺画や戦争画もありましたが、個人的にはあまり惹かれませんでした。また、本展には肉筆画も数多く出品されています。こちらも悪くはないのですが、同時期に活躍した素晴らしい画家はたくさんいるのでねえ・・・。やはり小林清親の真骨頂は版画ではないかと思います。

日本橋夜

日本橋夜


亀井戸藤

亀井戸藤


東京両国百本杭暁之図

東京両国百本杭暁之図


猫と提灯

猫と提灯


鉄砲打猟師

鉄砲打猟師


練馬区立美術館は、駅から近いのに長閑で雰囲気のとても良い所ですね。今年は開館30周年に当たるそうで、本展はその記念企画の一環のようです。秋にはシスレー展が開催されるので再訪したいと思います。