神戸野外実習

昨年度は人文地理学野外実習の授業は中止になってしまいましたが、今年は際どいところで実施することができました。

場所は学生の希望で神戸市となりました。

2/11 13:00神戸市内のホテルに集合


行きの新幹線内では今昔マップで古い地図上に軌跡を表示して遊びました。

ホテル集合後、最初は芦屋市の有名な高級住宅地、六麓荘の見学に出かけました。
私は3回目ですが、住宅展示場みたいで楽しいところですね。
画像

これは有名な住宅案内図です。

画像

遠くに300mの高層ビル、阿倍野ハルカスが見えます。やはり高級住宅地は眺めがよくないといけません。

続いて、六甲山頂に六甲ケーブルで上りました。
画像

このケーブルカー、一両は窓があるのですが、もう一両は吹きさらしで、冬の吹きさらしは寒いかったです。

画像

ケーブルカーの上からバスで少し行った六甲ガーデンテラスです。
風が強く激寒だったのですぐ次のバスで戻りました。寒かったですが、初めて六甲山上に来れてよかったです。居留地の外国人にも気軽に行けるいい避暑地だったのではないでしょうか。

夜はハーバーランドのumieで晩御飯を食べました。umieってイオンが運営してるんですね。

2/12
2日目は、まず阪神・淡路大震災記念 人と防災未来センターに出かけました。
画像

小学生や中学生が見学に来てました。内容は今一つでした。前半に地図を使った説明が無いのが残念でした。

三宮で昼食をとって、北野異人館にでかけました。定番観光地ですが、ここも初めてでかけました。
画像

坂道が多くて疲れましたが、高級住宅地は一般に坂の上にあるので仕方ない。

続いてシティー・ループというバスに乗って神戸港のメリケンパークへ。
ここには阪神淡路大震災時に崩れた岸壁が保存されています。
画像


次に神戸ベイクルーズの遊覧船に乗って神戸港を一周しました。
画像

川崎重工・三菱重工の造船所があり、それぞれ自衛隊の潜水艦がドックに入っていました。写真では黒くてわかりにくいですが。

画像

ポートアイランドに作られた大学です。きれいな建物です。

画像

海上保安庁の巡視船が停泊していました。

埠頭に戻ってから、神戸海洋博物館・カワサキワールドを見学しました。ここは面白かったのでもう少し長時間いてもよかったかも。

夕方になったので、夕食をとるために三宮に向かいます。
画像

結局、南京町で中華料理を食べました。昨晩も中華料理だったけど....。


2/13
三日目は長田の土地区画整理事業地区を見てから西神ニュータウン方面へ。
新長田駅周辺は以前来たときよりも建物が増えていました。

画像

漫画家の故・横山光輝さんが神戸市出身だそうで、鉄人28号や三国志キャラで売り出しているそうです。

続いて西神ニュータウンへ。
画像

西神中央駅は郊外のターミナルとしては過剰なまでに立派です。
周辺のマンション・戸建て住宅を見学してから、バスで近くの江崎グリコの工場グリコピアに工場見学に。
西神ニュータウンは住宅地だけでなく工業団地も含まれています。
画像

係の人に写真を撮ってもらったのですが、その際お姉さんの丸い帽子が風で飛ばされ、どこまでも転がっていくので、ひやひやしました。

画像

ポッキーの製造工程を見学しましたが、ヨーロッパでは"MIKADO"という名前で売られてるんですね。

これで全行程終了、西神中央駅で解散となりました。

翌日は関東で大雪になったので、タイミングよく実施できてよかったです。

タイルマップで緯度経度からタイルを求める

Googleマップで採用されている地図画像の形式であるタイルマップは、地理院地図でも使われるなど地図をネットで公開する際に便利な形式です。

これをgoogle Maps APIなどから使う際には、自動的に表示範囲の地図画像が読み込まれて表示されるので、詳しい仕組みを知らなくてもいいのですが、スタンドアローンのソフトで表示する場合には、自作する必要があります。

そこで、必要そうな式を作ってみました。

まず、メルカトル図法に関して、その投影法の式は次のようになっています。

メルカトル図法の式
緯度lat 経度lon からメルカトル図法の座標に変換(Rは地球の半径)

X=π*R*(lon/180)
Y=-R*Log(Tan((45+lat/2)*π/180))

lat=90になるとタンジェントの値が無限大となるので、メルカトル図法の地図では極点を示せません。

タイルマップの表示範囲と分割
タイルマップは、このメルカトル図法を東西方向と南北方向を正方形で切り出して、内部を分割し(ズームレベル)、小縮尺から大縮尺までの地図を表示します。

地理院地図の説明 http://portal.cyberjapan.jp/help/development.html#siyou-zm

ズームレベルを1つあげると、ファイルの数は4倍になります。

タイルマップの南北端
メルカトル図法では、南北端は無限大なりますが、タイルマップでは南北に限界があります。
メルカトル投影後に南北方向の長さを東西方向に合わせ、正方形に切り出します。

東西方向の幅は、半球なら
X = π*R*(180/180) = π*R
となります。
これが赤道から南北端までの距離に等しくなる緯度を次式から求めていくと、
R*Log(Tan((45+lat/2)*π/180))=R*180*π/180
Log(Tan((45+lat/2)*π/180))=π
Tan((45+lat/2)*π/180)=e^π
Tan(π/4+π/360*lat)=e^π

lat=2 * atn(e^π) * 180 /π - 90
= 85.0511287798067

と求まりました。

緯度経度とズームレベルを指定して、当該タイルのXY値を取得
さて、いよいよ本題です。表示されている地図にタイルマップのタイルを表示するには、緯度経度とズームレベルを指定してタイルを特定する必要があります。

ズームレベル(z)ごとの縦横の分割数は、2^zとなります。
そのため、メルカトル投影後のタイルの幅・高さ(W(z))は、

W(z)=(2*π*R)/2^z

となります。

そこで、経度lonのX値は

X(lon,z)=int((π*R*(lon/180)+π*R))/W(z))
 =int(π*R*(lon/180+1))/W(z))
 =int(π*R*(lon/180+1))/(2*π*R)/2^z))
 =int((lon/180+1)*2^z/2)

緯度latのY値は、

Y(lat,z)=int((-R*Log(Tan((45+lat/2)*π/180))+π*R)/W(z))
 =int((R*(-Log(Tan((45+lat/2)*π/180))+π)/W(z))
 =int((R*(-Log(Tan((45+lat/2)*π/180))+π)/(2*π*R)/2^z)
 =int(((-Log(Tan((45+lat/2)*π/180))+π)*2^z/(2*π))

となります。地球の半径Rは計算では使わなくてもすみます。


指定したタイル(x,y,z)の北端緯度、西端経度を求める

指定した位置のタイルが特定されたら、当該タイルの四隅の緯度経度を求めることで、地図上に表示できます。

西端経度は、

lon(x,z)=((x*W(z)/(2*π*R))*360-180
 =((x*(2*π*R)/2^z/(2*π*R))*360-180
 =(x/2^z)*360-180

北端緯度は、まず投影後の座標mapyを求め、

mapy=(y / 2 ^ z) * 2 * π- π

メルカトル投影の逆変換を行って求めます。

lat(y,z)=2 * atn(e^(-mapy)) * 180 / π - 90

これをプログラムに組み込めば、自作ソフトにタイルマップを表示できるようになります。

「今昔マップ on the web」のWebサイトへの埋め込み

時系列地形図閲覧サイト「今昔マップ on the web」でインラインフレームを使って今昔マップの地図をWebサイトやブログに埋め込めるようになりました。



「共有」ボタンをクリックすると、リンクURLと埋め込みコードが表示されます。

リンクをクリックして表示させたい場合は、「リンクURL」を貼りつけます。

Webサイトやブログに埋め込みたい場合は、「埋め込みコード」をコピーして貼りつけます。
フレームの高さと幅は、iframeタグの後ろの、widthとheightのパラメータで設定してください。

幅が狭い場合は、2画面表示でない状態で表示すると見やすくなります。その場合、dual=のパラメータをfalseにします。

左端の操作パネルは、埋め込んだ場合も最初は表示されています。

国土地理院の標高API

国土地理院のWebサービスである標高APIは、緯度経度を送るとその地点の標高を返してくれる便利なサービス。

http://portal.cyberjapan.jp/portalsite/info/dem.html

Google Maps APIでも標高を取得できますが、日本国内なら国土地理院のデータがより確実でしょう。

JavaScriptで取得する場合は次のようなコードで取得できます。
関数に緯度と経度を渡すと、標高とデータソースが取得できます。

function getAltitude(lon, lat){
var xhr = new XMLHttpRequest();
xhr.onreadystatechange = function(){
if (xhr.readyState === 4 && xhr.status === 200){
var json_data = eval( '('+xhr.responseText +')');
var txt = xhr.responseText + '\n\n';
txt += '標高:' + json_data.elevation + 'm\n';
txt += 'データソース:' + json_data.hsrc + '\n';
txt += '経度:' + lon + '\n';
txt += '緯度:' + lat + '\n';
var result = document.getElementById('data_result');
result.value = txt;
}
};
var url = "http://cyberjapandata2.gsi.go.jp/general/dem/scripts/getelevation.php?lon=" + lon + "&lat=" + lat +"&outtype=JSON";
xhr.open('GET', url);
xhr.send(null);
};


実際に試したい方はこちら

※現在このサンプルは動きません。CORSによる制限と思われます。
ただし、IEでローカルで実行した場合は、なぜか取得できます。

Google Maps APIのエンコード化ポリラインアルゴリズム

Google Maps APIでは線や多角形を描けますが、緯度経度の座標が多くなると、ファイルが座標で埋め尽くされてしまいます。そこで、座標を圧縮する方式として、エンコード化ポリラインアルゴリズムが用意されています。この方法では、座標列は文字列に置き換えられます。

Google Maps API V3では、ジオメトリライブラリを組み込んで、 google.maps.geometry.encoding.encodePath(path)でエンコーディングされた文字列を取得することができ、decodePath() により元の座標列に戻すことができます。

この機能を使って、エンコード化するWebサイトを作成しました。
http://ktgis.net/lab/software/encoding/index.html

しかし、JavaScriptを使用せず、他言語からエンコードするには、アルゴリズムを理解しておく必要があります。その点については、Google Developersに親切な説明があります。
https://developers.google.com/maps/documentation/utilities/polylinealgorithm?hl=ja

これを順番に見ていきます。※は私のコメントです。

1.最初の符合付き値を取得します:
-179.9832104

2.10進値を取得して 1e5 で乗算すると、結果は丸められます:
-17998321
※ここで注意すべき点は、精度が小数点以下5桁という点です。

3.この10進値を2進値に変換します。負の値は2の補数を使用して計算する必要があります。それには2進数を反転し、結果に1を加えます:
00000001 00010010 10100001 11110001
11111110 11101101 01011110 00001110
11111110 11101101 01011110 00001111

4.2進値を左に1ビットシフトします:
11111101 11011010 10111100 00011110
※3と4は、2進数に直して演算するように書いてあるりますが、実際は普通に元数値を2倍するだけであす。-17998321*2=-35996642

5.元の 10 進値が負の場合は、このエンコーディング結果を反転します:
00000010 00100101 01000011 11100001
※これもビット演算子で反転させることができます。
JavaScriptなら~-35996642=359966411

6.2進値を5ビットの集合に分割します(右側から):
00001 00010 01010 10000 11111 00001
※ここから2進数文字列として操作する必要があります。

7.5ビットの集合を逆の順序に並べ替えます:
00001 11111 10000 01010 00010 00001
※文字列を分割して個別の6つのバイト変数に入れ、並べ替えます。

8.後続のビット集合が続く場合は、各値に 0x20 の論理和演算を行います:
100001 111111 110000 101010 100010 000001
※「後続のビット集合が続く場合」とは、6バイト目だけでなく、後ろから00000が続く場合はその直前も含まれます。そうしたケース以外は、32を加算します。

9.各値を10進値に変換します:
33 63 48 42 34 1

10.各値に 63 を加算します:
96 126 111 105 97 64

11.各値を対応する ASCII 文字に変換します:
`~oia@
※8.で書かれているように、後ろに0があるケースでは、6文字よりも少なくなります。
※元数値が0の場合は、説明には書いてありませんが、アスキーコード63の?となります。


こんな感じで変換し、緯度の場合も経度の場合も同じです。そして、緯度経度の順に文字列を追加していきます。

ただし、座標列の2地点目からは、前の地点からのオフセット値を使用します。

また、注意すべき点として、変換後の座標の文字列に\(円マーク)が含まれることがあります。JavaScriptでは、\が特殊文字として扱われるため、JavaScriptの文字列変数に代入するには、\\としておかなくてはいけません。


最後に、精度の問題があります。2.のように、小数点以下5桁までの精度となっており、その上2つめの座標からオフセット値となるため、誤差が発生していきます。単独のポリゴンであれば、目立ちませんが、隣接するポリゴンの場合、拡大して表示すると、重複したり、隙間ができたりします。

この問題を避けるためには、実際の座標値と、小数点以下精度5桁でのオフセット値の累積を比較し、誤差が大きくなったらオフセット値を調整してやる必要があります。

ktgis.netのWebサイトへの接続方法

※9/26現在、おおむね復旧して問題は解消したようです

昨日24日から、私の諸々のWebサイト ktgis.net に接続できなくなっています。

これは、ドメインの更新を忘れていたためで、手続きをしたのでそのうち回復すると思いますが、それまでは

http://sv53.wadax.ne.jp/~ktgis-net/

から接続して下さい。
ただし、一部のページについてはうまく表示できないことがあります。

ジオコーディングは使えません。

また、「今昔マップ2」の地図画像のダウンロードは出来なくなっています。

ご不便をおかけしますが、復旧までしばらくお待ち下さい。

人口ピラミッド作成サイトを公開ました

人口ピラミッドは社会科や地理の資料でよく使われますが、Excelで作るのは結構面倒です。
参考:Excel2007で人口ピラミッド作成 http://standardization.at.webry.info/200812/article_6.html

そこで、男女・年齢別人口を貼りつけるだけで人口ピラミッドを作成するWebサイトを作成しました。
http://ktgis.net/lab/etc/pop_pyramid/index.html

この作成にはHTML5のグラフィック表示機能を利用しているので、古いブラウザ(Internet Explorer8以前など)では表示できないことがあります。

使い方は、Excelから年齢別人口をコピーしてテキストボックスに貼りつけるだけで簡単です。色の設定なども可能です。
取りあえず「サンプル表示」を試して頂ければと思いますが、参考までに2010年国勢調査による各歳人口ピラミッドを表示してみました。

人口ピラミッドでは5歳階級なども使われますが、各歳で見るとコーホートごとの戦争や丙午による出生数の減少、ベビーブームの影響が読み取れて興味深いものがあります。

全国:二つのベビーブーム世代が特徴です。第三次ベビーブームはなさそうです。
画像


東京都:20~30歳代が多くのですが、合計出生率が低く、子どもが少ないという特徴があります。
画像


埼玉県:郊外型です。35歳前後の第二次ベビーブーム世代が多くなっています。
画像


群馬県:平均型です。全国の人口ピラミッドとだいたい同じ形状をしています。
画像


青森県:地方型です。20歳代が少なく、35歳前後の第二次ベビーブーム世代の規模が小さいことが特徴です。
画像


沖縄県:合計出生率が高く、子どもが多くなっており、日本国内では特異な形状を示します。
画像


類型についてはこちらを参考にして下さい。
谷 謙二(2008)1920年から2005年にかけての都道府県ごとの年齢構造の変化とその類型化-コーホートごとの人口分布変動-.埼玉大学教育学部地理学研究報告,28,1-24.