この時期恒例の府中市美術館の「春の江戸絵画まつり」。今年の切り口は「ファンタスティック」。

「江戸絵画の19世紀」ちらし


江戸絵画のテーマにカタカナ語、しかもファンタスティックって日常的にはそんなに使わない言葉だと思うのですが・・・。日本語にしてしまうと「素晴らしき江戸絵画」といった感じで、これじゃあテーマとしては弱すぎるので強引にカタカナにしたのでしょうかねえ。そろそろネタ切れでしょうか?

まあテーマは何であれ、個々の作品はファンタスティックなので楽しめました。
前期展示は4月10日(日)まで。後期は4月12日(火)から5月8日(日)まで。前後期で全作品の展示替えをするようです。

八尾狐図
≪狩野探幽≫


四条河原夕涼図屏風
≪山口素絢≫


蛙の大名行列図
≪川鍋暁斎≫

 ジョギング中、突如猛ダッシュしてきた大きな犬にタックルを食らい、脇腹を激しく地面に打ちつけて肋骨を5本骨折という重傷を負ってしまいました(犬の散歩中はリードをつけましょう!)。
 2週間経過し大分良くなってきたので、「はじまり、美の饗宴展 すばらしき大原美術館コレクション」を観覧に国立新美術館へ。当事故で腰も痛めたので、さすがに「ジョギングで」というわけにはいきませんが・・・。

はじまり、美の饗宴展 すばらしき大原美術館コレクション ちらし

はじまり、美の饗宴展 すばらしき大原美術館コレクション ちらし


 大原美術館は高校の修学旅行で観覧しましたが、当時は絵画には全く興味がなくピカソの絵を見たことくらいしか記憶に残っていません。ギュスターヴ・モローの作品を所蔵している日本では数少ない美術館の1つなので、いつかは再訪せねばと思っていたのですが、なかなか機会がなく・・・そんな中、当展が開催されることとなりモロー作品「雅歌」も出品されるとあれば、これは逃すわけにはいきません。
 目当ての「雅歌」は結構小さなサイズでしたが、モロー作品だとすぐに分かる繊細で神秘的な作品。私としてはこれだけで満足・・・なのですが、さすがに日本有数のコレクションを誇る美術館ですので他にも良い作品がたくさん。でも、以前に他の特別展で見たことのある作品が結構多かったです。つまり、大原美術館は所蔵作品を気前よく貸し出す余裕があるってことなんでしょうね。やはりいつかは再訪しなければ。
 4月4日(月)まで。

雅歌

雅歌
≪ギュスターヴ・モロー≫


りんご採り

りんご採り
≪カミーユ・ピサロ≫


風景

風景
≪ポール・セザンヌ≫

1月9日(土)
ユニークな企画展が度々開催される町田市立博物館。ただ、積極的にPRしないので知らないままで見逃してしまいがち・・・。
今回、町田に出向く予定が出来たため調べてみた所、ちょうどこの日から「懐中時計展」が始まるということで、観覧してきました。

懐中時計展 ちらし

懐中時計展 ちらし

特に時計に興味があるわけではない私でも、100点超の精巧で綺麗な時計の展示を楽しめました。
「アメリカでは鉄道衝突事故がきっかけとなり、鉄道時計として懐中時計が発展した」など、解説も興味深いものでした。

出展されている時計は、個人コレクターの方が亡くなられた後、ご家族より当博物館に寄贈されたものだそうです。どれも保存状態が良くピカピカに磨かれていて、元の持ち主の懐中時計への深い愛情が感じられました。

本展は写真撮影可。3月6日(日)まで。

懐中時計展

懐中時計展

懐中時計展

懐中時計展

懐中時計展

懐中時計展

懐中時計展

懐中時計展

懐中時計展

懐中時計展

懐中時計展

懐中時計展

 新年の恒例になりつつあるジョギングでの「七福神めぐり」。一昨年は日野、昨年は八王子でしたが、今年は調布にしました。

 7寺のうち5寺は比較的近い所に固まっているのですが残りの2寺(明照院、昌翁寺)が離れています。そのため、公式サイトでは、途中にバスや電車を利用するルートが紹介されています。
 全行程を歩けば(走れば)、どんなルートを選択しても10キロを超えるかと思われます。私は一番西の西光寺からスタートしました(というか、自宅からここまでもジョギングですが)。この場合はゴールは東の2寺(明照院、昌翁寺)のいずれかにするのが最短ですが、最後は深大寺で締めたかったので長めの距離になりました。

 調布の七福神めぐりは日野と似た感じです。新撰組と縁のある土地である点(日野は土方歳三、調布は近藤勇の出身地)や、都内有数の初詣参拝者数を誇る寺が含まれる点(高幡不動と深大寺)など・・・。
 でも調布の方は、あまりPRしてないようです。日野では随所に見られた「七福神めぐり」の幟(目印になります)の数も少なく、目的の寺が見つけられずに探し回ることも・・・。特に明照院が分かり辛かったです。
 日野や八王子では私同様の七福神めぐり客があちこちに見られましたが、こちら調布は少なめ。そのせいか朱印を押してもらうために、呼び鈴を押して寺の方に出てきてもらわなければならない所もいくつかあり、折角のくつろいだお正月の邪魔をしてしまっているようで申し訳ない感じでした。

 一番の難関は最後の深大寺でした。寺の場所はすぐ分かりますし、分からなくても大勢の参拝客の流れについていけば勝手に到着するのですが、朱印の場所が分からない・・・スタッフの方も殆どが臨時の雇いなので聞いても分からず、数少ない幟を頼りにたどって行っても全然見つからない。大混雑の中を右往左往の末、ようやく見つけた朱印所は幟とは全く逆の方向にありました。
 もう少し分かりやすく案内して欲しい所ですが、七福神めぐり客自体が少ないので仕方がないのか・・・。
 まあ、この程度の不便は「お正月だから」ってことで全然許せちゃうんですけどね。

 参考:色紙代が700円、朱印料は300円(×7)なので、費用は計2,800円です(2016年1月現在)。

西光寺

西光寺


近藤勇坐像

近藤勇坐像(西光寺)


大正寺

大正寺


大正寺 門の前の行列

大正寺への参拝・・・ではなくて、奥にある布多天神社への参拝客の行列でした


祇園寺

祇園寺


常性寺

常性寺
深大寺以外で賑わっていたのはここだけでした。


明照院

明照院


昌翁寺

昌翁寺


深大寺

深大寺 当然の如く激混みでした


 
調布七福神色紙

調布七福神色紙


 

 12月26日(土)
 渋谷のBunkamuraザ・ミュージアムへ・・・今週から始まったばかりの「英国の夢 ラファエル前派展」を早速観覧してきました。

英国の夢 ラファエル前派展 ちらし

英国の夢 ラファエル前派展 ちらし


 ラファエル前派展は昨年の春(森アーツセンターギャラリー)に引き続いて2年連続ですが、前回がテート美術館からの出品だったのに対し、今回はリバプール国立美術館(リバプール市内・近郊の3美術館等の総称だそうです)所蔵の作品なので、重複出品はなく新鮮味が薄れるなんてことはありません。まあ、ロセッティの作品に関しては全部同じように見えてしまうのですが・・・。

 一番のお目当てはエヴァレット・ミレイの作品・・・全部で8点も出品されていました。さすがはミレイ、どれも素晴らしかったですが、中でも『ブラック・ブランズウィッカーズの兵士』は天才っぷりを一番発揮している作品だと思います。しばらく見とれてしまいました。
 本展はミレイ以外にも良作がずらり(ミレイ展ではないので当然なのかもしれませんが)。しかもサイズの大きな作品が多いので見応え充分。かなりお奨めです。
 2016年3月6日(日)まで。

ブラック・ブランズウィッカーズの兵士

ブラック・ブランズウィッカーズの兵士
≪ジョン・エヴァレット・ミレイ≫

女性のドレスの光沢・質感が素晴らしい!
また、この兵士の敵に該当するナポレオンの肖像(ダヴィッド作)が描かれているのもいいですね。

お気に入りの詩人

お気に入りの詩人
≪ローレンス・アルマ=タデマ≫

書見台での学習

書見台での学習
≪フレデリック・レイトン≫

テラスにて

テラスにて
≪エドワード・ジョン・ポインター≫

強く惹かれた作品の1つ。透明感がたまりません!

聖セシリア

聖セシリア
≪ジョン・メリッシュ・ストラヴィック≫

別の1点『おお燕よ、燕』を含め、ストラヴィックの作品にはギュスターヴ・モローと似た世界観を感じました。
(個人的な感想ですが・・・)

小さな召使い(乙女エレン)

小さな召使い(乙女エレン)
≪エレノア・フォーテスク=ブリックデール≫
ラストに展示されていた作品。最後まで見応えたっぷりの展覧会でした。

 貴金属や宝石の類には全く興味がないので、「黄金伝説展」というタイトルを聞いた時点では行く気ゼロだったのですが、ギュスターヴ・モローの作品が出品されていると知り、その観覧のみを目的に上野の国立西洋美術館へ。

ゴーギャンとポン=タヴァンの画家たち展 ちらし

黄金伝説展 ちらし

 出品されているモロー作品は『イアソン』『アルゴー船の乗組員』『ヘラクレスと青銅の蹄をもつ鹿』の3点。
 この中でも目玉の『イアソン』は昨年のオルセー美術館展でも見ているのですが、次はいつ見られるか分からないのでじっくりと・・・。やはり何度見ても感動モノです。鳥肌が立ちます。
 『アルゴー船の乗組員』と『ヘラクレスと青銅の蹄をもつ鹿』は初見。『イアソン』ほどの迫力はありませんが素晴らしいことには変わりなく、見られただけで幸せです。

 モロー以外の絵画も何点か・・・洋画で黄金といえば、やはりクリムト。そのクリムト作品『人生は戦いなり(黄金の騎士)』は2年半前に宇都宮美術館で見ているので、感動は半分でした。タイトルが『ダナエ』という作品もあったのですが、これがクリムトの『ダナエ』だったら最高なんですけどねえ・・・それはやっぱり無理ですかね。本来のメインである金の装飾品の数々は、ざーっと観覧しただけなので特に感想はありません。

 なお、この特別展のチケットで常設展も観覧することが出来ます。ここの美術館の常設展は、他の美術館の特別展並みもしくはそれ以上に充実した内容なのでお得です。金製品に全く興味がない私にも、観覧料1,600円の価値は十分にありました。
 1月11日(月・祝)まで。

イアソン

イアソン
≪ギュスターヴ・モロー≫


アルゴー船の乗組員

アルゴー船の乗組員
≪ギュスターヴ・モロー≫


ヘラクレスと青銅の蹄をもつ鹿

ヘラクレスと青銅の蹄をもつ鹿
≪ギュスターヴ・モロー≫

11月28日(土)
ここ最近は美術館に行く時間が全くとれず、禁断症状が出始めていたので、それを解消すべくパナソニック汐留ミュージアムへ。現在開催中の「ゴーギャンとポン=タヴァンの画家たち展」を観覧してきました。

ゴーギャンとポン=タヴァンの画家たち展 ちらし

ゴーギャンとポン=タヴァンの画家たち展 ちらし

個人的にはゴーギャンはそんなに好きではないのですが、ポール・セリュジエやモーリス・ドニなどゴーギャンの影響を受けたナビ派の画家は好きなので、それらの作品をまとめて見られる本展は見逃せません。

有名な作品は出品されていませんし、点数もさほど多くはありません(計74点)が、良作揃いですし解説も丁寧なので見応え充分でした。

ポン=タヴァンというのはブルターニュ地方にある町だそうで、ブルターニュと聞くと「ワインの美味しいところかな?」なんて思ってしまうのですが、それはブルゴーニュでした(間違えやすい?)。
この辺りは有名な観光地は少なそうですが(モン・サン=ミシェルは近いですが)、作品を見ているととても魅力的な場所に思えてきて訪れてたい気になります。これは、本展に限った話ではありませんが・・・。
12月20日(日)まで。

ブルターニュの眺め

ブルターニュの眺め
≪ポール・ゴーギャン≫


呪文或いは聖なる森

呪文或いは聖なる森
≪ポール・セリュジエ≫


ブレストの港

ブレストの港
≪モーリス・ドニ≫


礼拝行進

礼拝行進
≪アンリ=ガブリエル・イベルス≫

10月7日(水)
 「モネ展」を観覧しに東京都美術館へ・・・。

東京都美術館

東京都美術館


 「12月中旬まで開催されているので特に急いでいくこともないや」と思っていたのですが、本展の目玉作品『印象、日の出』の展示は10月18日までということを最近知り、これを見逃すわけにはいかないと、平日にも関わらず急遽向かいました。

 当作品は21年ぶりの来日とのことですが、21年前は絵画には全く興味を抱いていなかったので、私にとっては今回が初見です(パンの作品は2年半前に観覧しましたが)。21年前には、当作品のことを知ってすらいなかったような・・・ ^^;

 平日(水曜日)の昼にも関わらず大混雑。3年前にフェルメールの『真珠の耳飾りの少女』が来日した時ほどではないにしても、やはり『印象派』の由来となった有名作品だけあって凄い集客力です。まあ、それでも入場時の待ち時間は10分程度で済んだので、平日を選んで正解だったようです。ちなみに一番長く待ったのはグッズ売り場のレジです。
 12月13日(日)まで。

印象、日の出

印象、日の出
今回は今作を見られただけで満足です。


霧のヴェトゥイユ

霧のヴェトゥイユ
霧がかかったようなような作風が特徴のモネが「霧」を描いているという・・・That's Monet! といった感じ


ジヴェルニーの黄色いアイリス畑

ジヴェルニーの黄色いアイリス畑
『日傘の女』の背景と似た印象。カミーユの姿はありませんが・・・


睡蓮

睡蓮
これまで、モネの睡蓮は数十枚みてきましたが、ここまで抽象的な作品は初めてです。
この他、最晩年の作品には抽象画風のものが多く楽しめました。

2015年10月3日(土)
 今年のグローバルフェスタは昨年までの日比谷公園から会場を変えて、お台場にて開催。その関係なのか、恒例だったチャリティーランが行われないのは少々残念ですが、このイベント自体が好きなので見に行ってきました。
 
 普段はあまり馴染みのないアジアやアフリカ、中南米などの食べ物や文化に触れることが出来て楽しめます。また、出展されている各団体の皆さんの生き生きとした姿を見るのも楽しさの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という名称で保存・使用しています)