リアルGeoHexにいってきました！

Togetter - 「リアルGeoHexを探せ！」に登場する
毛馬（けま）にあるリアルGeoHex、淀川大堰（よどがわおおぜき）を見てきました。

Wikipediaの淀川大堰の項目には

淀川大堰へは阪急千里線・柴島駅（くにじまえき）で下車するか、自動車では阪神高速12号守口線・長柄出入口下車が最寄となる。

と書かれているのですが、
地図でみて1kmぐらいありそうだし、ヘタレなのでもっと近いのを探したところ
大阪市バスの37番がすぐ近くを通る事を発見しました。
大阪市交通局のバスのページは天才的に分かりにくいデザインになっていて
ホームページにある路線図ではその路線が確認できず、
PDFの路線図と乗り場案内で
ようやくたどり着けたのでした。

大阪駅前のバス乗り場

これに写っている緑っぽいバスに乗っていきます。

大阪駅前のバス乗り場は、乗り場案内ではオープンスペースっぽいですが
上に建物が建っていて、二階にちまちまお店があります。
わりと寂れてて、コインロッカー探してるときはいくといいと思うよ。

バス停。37番井高野車庫前行き

途中全部省いて、15分ほどで長柄橋北詰バス停に到着。
右も左も何にも無いとこで、ちょっとビビる。

左（北東）側

右（西南）側

本当に何も無い上に、この高くなってる道路が邪魔で周りが全然見えない。

そこでこんな風に歩くのですが...

トンネル

歩いていいか迷ってたらおっちゃんがスタスタ入っていったので
ついていきました。

明かりの向こうにはHex!とおもいきや、まだ土手。

道もないのでのぼることに。

上に橋っぽいのがみえる！


登りきると

淀川大堤！
最初に思ったのはただただ「でかい」そんだけです。

写真ではスケール感がすごく分かりにくいのですが、
堤の高さは4階か5階建てのビルくらいあります。

さらにこのすぐ下流側には橋となんかのパイプが合計3本通っていて
それぞれ鉄橋なもんでとても巨大感にあふれる場所です。
橋マニアは好きかもしれません。

写真下部にHexが！黒いのは石っぽい。
もうこの目下にはHexが広がっています。

あと、Hexの上でバーベキューしてる団体が見えます。

この時点で「でかい」と何回思ったか...。
せいぜい数人乗れるくらいだろうと思っていたHexひとつに
余裕で15〜16人がバーベキューしてスカスカというサイズ。でかい。。


そしてついにHexの大地に立つ！

でかすぎてよくわかりませんがな。。



気を取り直して、観察してみましょう。
Hexはコンクリの淵で6角形がつくってあって
その中に入っているか、かぶせてある物が3種類ぐらいあるようです。

黒い石のHex、茶色っぽい石のHex、草のHex

黒っぽい石は形も大きさもバラバラ。

固定もされてなくて、礫石っぽいので、大きい写真だと河原で撮ったもの
といっても通じるでしょう。

茶色っぽい石は整然と並んでいます。

サイズはだいたい30センチ四方になっているようです。

これも固定はされていないようですが
さすがに掘り返すのは怒られそうなのでがまんしました。


そして4番目のタイプ。土に石が埋まっています。

これらの差が何なのか私には分かりません。

Hexのサイズは一片5m20cmほど。

Hexの構造のコンクリはちょうど1m。
大きいわけです。


ふらふらする間にきになるアレも発見。

きになるアレ

ふ、古い。
もう壊れそうなくらいぼろぼろ。

ふたのマークは大阪市のもの。

もういっこきになるアレ

影が。。


これまた3階建てのビルくらいあって、最初に大堤が見えた写真にも写っています。

写真ではスケールがよく分からないのですが
この塔の手前側のすこし色がついている部分で
大人の頭より上ぐらいの大きさです。
（よくみると塔の左側の奥に自転車に乗っている人が見えます。

とにかく、でかい。


Hexは地上で見るとこのへんで唐突に終わっているように見えます。

このへん。

地上で見てよくわかったのは、
6角形はYの字型のモジュール（？）を3つ組み合わせて作ってあることです。
6辺の各真ん中で分かれていて、継ぎ目があります。
GoogleMapでも拡大すると少し見える箇所があります。


しかしこの端はその継ぎ目も無視して作られていて、
よほどこのサイズにしなければならない何かがあった。。のかもしれません。



ここまでに文字が出てきたのは、大阪市のマークの気になるアレの
横にくっついていた「止水栓」のみ。
まだなんだか分からない私。


さて、気になるものその2に似た塔が上流側にもう一つあったので
そちらに向かってみました。

まわりは草ぼうぼうというほか表現しようがない状態。

しかしその途中のここで、何となくヒントを発見。

監視カメラ...じゃなくてこれは取水口ですね。

監視カメラは人の監視をしてるのかと思いきや
水量の監視が主なようです。

このわきに1979年竣工のプレートがありました。
やっと文字にであえた！


きになる塔その2は何も無く空振り。

どうしようもないくらい、塔。
あいかわらず文字なし。ちょっとボロい。

監視カメラ。
あー、川の様子でも見てるんだなー
って思ってました。

この周辺で私の行く手を阻むやたら育ちのいい草。
一番背が高いのは3mぐらいありました。

ふと川をみるとこんな管が。
このときはなんだか分からなかったけど、これもヒントになりました。

この辺です。

さて、このHexはいったい何者なのか。。次に続く...

GeoHex拡張のまえにおさらい

GeoHexのどこかを拡張する話があるのですが、
正直中身や周辺の理屈を理解しているわけではないので
現状を考えてみました。

GeoHexのコードと呼んでいる文字列、"quhO0cp" などは、
レベル（＝Hexのサイズ）文字列に続けて、
赤道を軸とした斜交座標系 [wikipedia] の交点を示すXY座標を
60進数に変換した値を表しています。

六角形はこの交点を中心に平面充填していきます。
（とりあえず地球の球面上での境界条件を無視して進めます）

図はフリーハンドです。すみません。。。

一つ大きなレベルでは六角形は一つ小さくなり、
XYの座標軸はもとのレベルのXY軸の間に入ります。

このとき、同じXY座標をもつHexはレベルが一つあがると
XとYの値に2をかけた物になります。
例えば、
level 16で 37502/-13386 のHexと同じXY座標の
level 17のHexは75004/-26772となります。
75004 = 37502*2、
-13386 = -26772
です。


さて、斜交座標と、レベルとは別に、
隣接Hexについて考えてみます。

通常の直交座標（普通のXY軸）では

こうなるもので、考えるまでもないのですが
斜交座標のHexでは次のようになります。

この、隣接するHexの数値は相対的に固定であり、
どこのHexを取り出しても同じことが言えます。

つまり、XY座標がわかると、
* 同じXYの交点を中心に持つそれ以上のLevelのHexのXY座標
が分かり、さらに
* 自分の周辺のHexのXY座標
* 同じ交点を中心に持つHexのXY座標
も簡単な加減算で求めることができます。
さらに、
* 小さなレベルと同じ交点を共有する場合そのXY座標
もわかります。

個数を数えてみると、
レベル差1では中心+周辺6つのHexで合計7つ、
レベル差2ではそのさらに周辺12で合計19のHexについて
簡単な計算でXY座標を求めることができます。
これはつまり、小さなレベル（大きなHex）とXY座標を共有しているHexでは
大きなHexはそのまま
それよりも大きなレベル（小さなHex）の複数のHexのインデックスになっており、
1つのインデックスで19個のHexを表すとは、
あるレベルでのXYに対して、
それと同じ交点を座標に持つ二つ小さなレベル（大きなHex）が
GeoHex的に存在していなくても
ある物とするコード、ではないかと思ったのでした。

緑が現状のGeoHexでは存在していないHex（インデックス？）

地理情報学会第19回研究発表大会にお邪魔してきました。

残念ながら写真無しです。


GISA（地理情報学会、GIS Assosiation of Japan）の研究発表に縁あってお邪魔させていただきました。
案内のページに情報があります。
しかしいつまで残ってくれるか分からないURLです。

経緯

私はGISAの会員ではなく、学生でも研究者でもありません。
学会の開催を知ってはいましたが、行く予定ではありませんでした。
一般参加費4000円という金額におののいたものの、
@tossetoさん、@ogugeoさんからやや強引（笑）にお誘いを受けました。
この時点ではGISAがどんな学会かも殆ど分かっていなかったので
公開されている発表タイトル一覧や
twitterでお知り合いになった方達の研究室いくつかから
論文を斜め読みさせていただいて、
ガチガチな学会でもなさそうなので行ってみることにしました。

会場

会場は立命館大学衣笠キャンパスでした。
当日まで衣笠が読めなかったレベルに全く分かっていませんでしたが
事前に @tosseto さんのお手伝いで衣笠キャンパスのタグづけをしていたので
なんとなく周辺把握ができてたどり着けました。

各セッション

いろいろなセッションを聴かせていただきました。
全くの素人でおそらく内容も殆ど分かっていないので、かなりおとなしめに過ごしました。
しかし英語のセッションとか、まだ聴ける耳が残ってたんだなぁとしみじみ。

こゆい

初日、 @ogugeoさん、 @niyalisut さんとお食事させていただきました。
なんともこゆいメンバーで、
とてもとても楽しくお話できました。
その後もお二方は常に気にかけてくださり、ほぼ全員初対面の私には
とてもとてもありがたく、助かりました。
この場を借りてお礼申し上げます。ありがとうございました。

そしてこゆい方たち

懇親会と、それに続くGeo飲みなどで多種多様な多くの方達にお会いすることができました。
ありがとうございました。
ここ数年で外に出て人にあうのは
二度目ぐらい（一度目は夏に @tosseto さんにお会いしたとき）
でしたが、とても良い経験をさせていただきました。

お会いできた方達: @DIG11さん、@geo80kさん、@heromiyaさん、@heshikikさん、@mapconciergeさん、@megumeruさん、@niyalistさん、@nissyyu さん、@nro2daiさん、@ogugeoさん、@picaosgeoさん、@pokopenさん、@sada1110さん、@Say_noさん、@tagchanさん、@tossetoさん、@wing83さん、@yan_yakkuさん、@yaskondoさん、@usuyuさん、ラガワン先生

一日中走り回っておられた@tossetoさん、本当におつかれさまでした。
今回参加できたのはほとんどが @tossetoさん、@ogugeoさんのおかげです。
ありがとうございました。

最後に一言つぶやいて
お礼に変えさせていただきます。

GeoHex v2 PHP

GeoHexのPHPポーティング (http://github.com/geohex/geohex-php) の
世界対応版(v2.01)で作成しました。
元になるjsがアップデートするとのことですのでまだお試し版です。
問題等ありましたら教えていただけると助かります。

2010/09/11追記
特定の値でloc2zoneがちょっとずれるバグがあるようです。
geohex2_coreのリファクタリングを取り込みつつ調査しています。

＜?php
/* GeoHex module for php */
/*  DISCLAIMER OF WARRANTY

BECAUSE THIS SOFTWARE IS LICENSED FREE OF CHARGE, THERE IS NO WARRANTY
FOR THE SOFTWARE, TO THE EXTENT PERMITTED BY APPLICABLE LAW. EXCEPT WHEN
OTHERWISE STATED IN WRITING THE COPYRIGHT HOLDERS AND/OR OTHER PARTIES
PROVIDE THE SOFTWARE "AS IS" WITHOUT WARRANTY OF ANY KIND, EITHER
EXPRESSED OR IMPLIED, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. THE
ENTIRE RISK AS TO THE QUALITY AND PERFORMANCE OF THE SOFTWARE IS WITH
YOU. SHOULD THE SOFTWARE PROVE DEFECTIVE, YOU ASSUME THE COST OF ALL
NECESSARY SERVICING, REPAIR, OR CORRECTION.

IN NO EVENT UNLESS REQUIRED BY APPLICABLE LAW OR AGREED TO IN WRITING
WILL ANY COPYRIGHT HOLDER, OR ANY OTHER PARTY WHO MAY MODIFY AND/OR
REDISTRIBUTE THE SOFTWARE AS PERMITTED BY THE ABOVE LICENCE, BE
LIABLE TO YOU FOR DAMAGES, INCLUDING ANY GENERAL, SPECIAL, INCIDENTAL,
OR CONSEQUENTIAL DAMAGES ARISING OUT OF THE USE OR INABILITY TO USE
THE SOFTWARE (INCLUDING BUT NOT LIMITED TO LOSS OF DATA OR DATA BEING
RENDERED INACCURATE OR LOSSES SUSTAINED BY YOU OR THIRD PARTIES OR A
FAILURE OF THE SOFTWARE TO OPERATE WITH ANY OTHER SOFTWARE), EVEN IF
SUCH HOLDER OR OTHER PARTY HAS BEEN ADVISED OF THE POSSIBILITY OF
SUCH DAMAGES.

*/
/*

GeoHex oliginaly written by @sa2da (http://geogames.net/, http://geohex.net/) in Javascript

Ported to PHP by Mage Whopper  at 2010.9.8

Copyright (C) 2009 sa2da (http://twitter.com/sa2da)
License : CC-BY-SA-2.1-JP
  You can find the full legal code at
  http://creativecommons.org/licenses/by/2.1/jp/legalcode (Japanese)
  or in the local file cc-by-sa-2.1-legalcode.html.
  such like http://creativecommons.org/licenses/by-sa/3.0/legalcode (English)
  Here is only an abstract :
 
   You are free:
        to Share – to copy, distribute and transmit the work
        to Remix – to adapt the work

  Under the following conditions:
        Attribution – You must attribute the work in the manner specified by the author or licensor 
             (but not in any way that suggests that they endorse you or your use of the work).
        Share Alike —  If you alter, transform, or build upon this work, you may distribute the
              resulting work only under the same or similar license to this one. 

   For any reuse or distribution, you must make clear to others
        the license terms of this work. The best way to do this
        is with a link to this
        (http://creativecommons.org/licenses/by-sa/3.0/) web page.

  Any of the above conditions can be waived if you get
        permission from the copyright holder.

  Nothing in this license impairs or restricts the author's moral rights.

  SPECIAL THANKS TO @_hfu_ for projection transform mathematical expression
  (http://twitter.com/_hfu_/status/23144088252)
 
*/

class GeoHex{
 var $h_key = 'abcdefghijklmnopqrstuvwxyz0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZ';
 var $h_base = 20037508.34;
 var $h_deg;
 var $h_k;
 var $h_size;
 var $h_x;
 var $h_y;
 var $h_lon;
 var $h_lat;

 function __construct(){
  $this->h_deg = pi()*(30/180);
  $this->h_k = tan($this->h_deg);
 }

 // this function only  for downward compatibility
 function geohex2latlng( $code ){
   return $ths->getZoneByCode($code);
 }

 // this function only  for downward compatibility
 function latlng2geohex($lat, $lon, $level){
   return $ths->getZoneByLocation($lat, $lon, $level);
 }

 function setHexSize($_level){
  if($_level < 1 || $_level > 24){
   return -1;
  }
  return $this->h_base/pow(2,$_level)/3;
 }

 function getZoneByLocation( $_lat, $_lon, $_level){
  $this->h_size = $this->setHexSize($_level);
  if($this->h_size < 0){
   die("invalid level error!");
  }
  $level = $_level;

  list($lon_grid, $lat_grid) = $this->loc2xy($_lon,$_lat);
  $unit_x = 6* $this->h_size;
  $unit_y = $unit_x*$this->h_k;
  $h_pos_x = ($lon_grid + $lat_grid/$this->h_k)/$unit_x;
  $h_pos_y = ($lat_grid - $this->h_k*$lon_grid)/$unit_y;
  $h_x_0 = floor($h_pos_x);
  $h_y_0 = floor($h_pos_y);
  $h_x_q = floor(($h_pos_x - $h_x_0)*100)/100;
  $h_y_q = floor(($h_pos_y - $h_y_0)*100)/100;
  $this->h_x = round($h_pos_x);
  $this->h_y = round($h_pos_y);

  $h_max=round($this->h_base/$unit_x + $this->h_base/$unit_y);

  if($h_y_q>-$h_x_q+1){
   if(($h_y_q<2*$h_x_q)&&($h_y_q>0.5*$h_x_q)){
    $this->h_x = $h_x_0 + 1;
    $this->h_y = $h_y_0 + 1;
   }
  }else if($h_y_q<-$h_x_q+1){
   if(($h_y_q>(2*h_x_q)-1)&&($h_y_q<(0.5*$h_x_q)+0.5)){
    $this->h_x = $h_x_0;
    $this->h_y = $h_y_0;
   }
  }
  $this->h_lat = ($this->h_k*$this->h_x*$unit_x + $this->h_y*$unit_y)/2;
  $this->h_lon = ($this->h_lat - $this->h_y*$unit_y)/$this->h_k;

  list($z_loc_x, $z_loc_y) = $this->xy2loc($this->h_lon,$this->h_lat);
  if(($this->h_base - $this->h_lon) < $this->h_size){
   $z_loc_x = 180;
   $h_xy = $this->h_x;
   $this->h_x = $this->h_y;
   $this->h_y = $h_xy;
  }

  $h_x_p =0;
  $h_y_p =0;
  if($this->h_x<0){ $h_x_p = 1;}
  if($this->h_y<0){ $h_y_p = 1;}
  $h_x_abs = abs($this->h_x)*2+$h_x_p;
  $h_y_abs = abs($this->h_y)*2+$h_y_p;
  $h_x_10000 = floor(($h_x_abs%77600000)/1296000);
  $h_x_1000 = floor(($h_x_abs%1296000)/216000);
  $h_x_100 = floor(($h_x_abs%216000)/3600);
  $h_x_10 = floor(($h_x_abs%3600)/60);
  $h_x_1 = floor(($h_x_abs%3600)%60);
  $h_y_10000 = floor(($h_y_abs%77600000)/1296000);
  $h_y_1000 = floor(($h_y_abs%1296000)/216000);
  $h_y_100 = floor(($h_y_abs%216000)/3600);
  $h_y_10 = floor(($h_y_abs%3600)/60);
  $h_y_1 = floor(($h_y_abs%3600)%60);

  $hl_key = $this->h_key;
  $h_code = "" . substr($hl_key, ($level%60), 1);

  if($h_max >=1296000/2) $h_code = $h_code . $hl_key[$h_x_10000] . $hl_key[$h_y_10000];
  if($h_max >=216000/2) $h_code = $h_code . $hl_key[$h_x_1000] . $hl_key[$h_y_1000];
  if($h_max >=3600/2) $h_code = $h_code . $hl_key[$h_x_100] . $hl_key[$h_y_100];
  if($h_max >=60/2) $h_code = $h_code . $hl_key[$h_x_10] . $hl_key[$h_y_10];
  $h_code = $h_code . $hl_key[$h_x_1] . $hl_key[$h_y_1];

  return array(
   "lat"     => $z_loc_y,
   "lon"    => $z_loc_x,
   "x"       => $this->h_x,
   "y"       => $this->h_y,
   "code" => $h_code
  );
 }

 function getZoneByCode($_code){
  $c_length = strlen($_code);
  $zone = array();
  $hl_key=$this->h_key;
  $level = strpos( $hl_key, $_code[0]);
  $scl = $level;
  $this->h_size =  $this->h_base/pow(2,$level)/3;
  $unit_x = 6* $this->h_size;
  $unit_y = 6* $this->h_size* $this->h_k;
  $h_max=round($this->h_base/$unit_x + $this->h_base/$unit_y);
  $this->h_x=0;
  $this->h_y=0;

  if( $h_max >= (1296000/2)) {
   $this->h_x = strpos( $hl_key, $_code[1]) *1296000+strpos( $hl_key, $_code[3] )*216000+strpos( $hl_key,  $_code[5] )*3600+strpos( $hl_key,  $_code[7] )*60+strpos( $hl_key,  $_code[9] );
   $this->h_y = strpos( $hl_key,  $_code[2] )*1296000+strpos( $hl_key,  $_code[4] )*216000+strpos( $hl_key,  $_code[6] )*3600+strpos( $hl_key,  $_code[8] )*60+strpos( $hl_key,  $_code[10] );
  }else if($h_max >=(216000/2) ) {
   $this->h_x = strpos( $hl_key,  $_code[1] )*216000+strpos( $hl_key,  $_code[3] )*3600+strpos( $hl_key,  $_code[5] )*60+strpos( $hl_key,  $_code[7] );
   $this->h_y = strpos( $hl_key,  $_code[2] )*216000+strpos( $hl_key,  $_code[4] )*3600+strpos( $hl_key,  $_code[6] )*60+strpos( $hl_key,  $_code[8] );
  }else if($h_max >=3600/2){
   $this->h_x = strpos( $hl_key,  $_code[1] )*3600+strpos( $hl_key,  $_code[3] )*60+strpos( $hl_key,  $_code[5] );
   $this->h_y = strpos( $hl_key,  $_code[2] )*3600+strpos( $hl_key,  $_code[4] )*60+strpos( $hl_key,  $_code[6] );
  }else if(h_max >=60/2){
   $this->h_x = strpos( $hl_key,  $_code[1] )*60+strpos( $hl_key,  $_code[3] );
   $this->h_y = strpos( $hl_key,  $_code[2] )*60+strpos( $hl_key,  $_code[4] );
  }else{
   $this->h_x = strpos( $hl_key,  $_code[1] );
   $this->h_y = strpos( $hl_key,  $_code[2] );
 }

  $this->h_x = ($this->h_x%2)?-($this->h_x-1)/2:$this->h_x/2;
  $this->h_y = ($this->h_y%2)?-($this->h_y-1)/2:$this->h_y/2;
  $h_lat_y = ($this->h_k*$this->h_x*$unit_x + $this->h_y*$unit_y)/2;
  $h_lon_x = ($h_lat_y - $this->h_y*$unit_y)/$this->h_k;

  list($this->h_lon, $this->h_lat ) = $this->xy2loc($h_lon_x,$h_lat_y);

  return array(
    "code" => $_code,
    "lat"     => $this->h_lat,
    "lon"    => $this->h_lon,
    "x"       => $this->h_x,
    "y"       => $this->h_y
  );
 }

 function loc2xy($_lon, $_lat) {
  $x = $_lon*$this->h_base/180;
  $y = log(tan((90+$_lat)*pi()/360))/(pi()/180);
  $y = $y*$this->h_base/180;
  return array($x, $y);
 }

 function xy2loc($_x, $_y) {
  $lon=($_x/$this->h_base)*180;
  $lat=($_y/$this->h_base)*180;
  $lat=180/pi()*(2*atan(exp($lat*pi()/180))-pi()/2);
  return array($lon, $lat);
 }

 function _gh_strpos( $num ){
  $hl_key = $this->h_key;
  return strpos($hl_key, $_code[$num]);
 }
 
}/* class geohex */

使い方

//モジュールの読み込み
require_once('geohex2.01.php');
 
// GeoHexオブジェクト生成
$ghx = &new GeoHex();
 
// 各計算
// 緯度経度からGeoHexコードを求める
print $ghx->getZoneByLocation($latitude, $longitude, $level) ."\n";
 
// GeoHexコードから緯度経度に変換
print $ghx->getZoneByCode($geohex)."\n";

Spatialiteのデータ保存形式

Spatialiteという、SQLiteをベースに
拡張したOpenGISベースの便利なDBシステムがあります。

SQLiteなのでお手軽！とおもいきや
SQLを直たたきで位置情報を取ったりすることはできなくて
PostGIS同様関数でgeometry情報を格納しています。

spatialiteの持つテーブルは

CREATE TABLE spatial_ref_sys (
 srid INTEGER NOT NULL PRIMARY KEY,
 auth_name VARCHAR(256) NOT NULL,
 auth_srid INTEGER NOT NULL,
 ref_sys_name VARCHAR(256),
 proj4text VARCHAR(2048) NOT NULL
)
CREATE TABLE geometry_columns (
 f_table_name VARCHAR(256) NOT NULL,
 f_geometry_column VARCHAR(256) NOT NULL,
 type VARCHAR(30) NOT NULL,
 coord_dimension INTEGER NOT NULL,
 srid INTEGER,
 spatial_index_enabled INTEGER NOT NULL
)

の二つで、地物をを地図に追加すると
地図のtableにalter tableで追加、、、ってそういう使い方ですか。

試しにspatialiteのページから
テストDBをDLして中をのぞいてみる。

sqlite> .tables
HighWays           Towns              geometry_columns 
Regions            geom_cols_ref_sys  spatial_ref_sys  

geometry_columns、spatial_ref_sysは
spatialiteがinitするときに作成しています。

geom_cols_ref_sysはVIEWなので割愛します。

これにHighWaysとTownsとRegionsの三つのデータが入っている感じなのですが
この中身は

INSERT INTO "Towns" VALUES(4847,'Druento',8235,1,0,0,X'0001787F000048E17A14BCAF17413D0AD733BE11534148E17A14BCAF17413D0AD733BE1153417C0100000048E17A14BCAF17413D0AD733BE115341FE');

こんな感じになっていて読めません。
スキーマではTownsは

CREATE TABLE Towns (
 PK_UID INTEGER PRIMARY KEY AUTOINCREMENT,
 Name TEXT,
 Peoples INTEGER,
 LocalCounc INTEGER,
 County INTEGER,
 Region INTEGER,
 "Geometry" POINT
);

という構造です。
POINTはspatilaiteの持つ形式で

0 - 7   X coordinate   a double value corresponding to the X coordinate for this POINT
8 - 15  Y coordinate  a double value corresponding to the Y coordinate for this POINT

という構造のBLOBをALTERTABLEで追加しています。

もうすこし言うと、
TownsというのはPOINT構造体の配列をBLOBに格納した物で、
DBの機能や最適化はガン無視しています。

私もどちらかというと力技でプログラムを書いてきた方なので
言いたいことは痛いほどわかります。
たとえばLineのデータを持つときに
一点一点ポイントを

CREATE TABLE point(
 point_id integer NOT NULL,
 line_id integer,
 x    integer,
 y    integer
);

CREATE TABLE linestring(
 line_id integer NOT NULL,
 number_of_points integer
);

などとするよりも、
（簡便のためlinestring以外のケースなどは省略しています）
データをバイナリで持って、処理する側で解析も変更もやったほうが
便利なこともあるでしょう。
まとまったデータでのパフォーマンスも制御しやすくなります。

そしてこの、バイナリを数値やテキストとして扱えるようにした
OpenGISのインターフェースがついたのがspatialiteです。

地図DBの内容を使うだけの立場を取る限りは
内部構造がどうであろうと自分が望む形式で取り出せればそれで問題ないでしょう。
しかしさしあたって、
DBはDBでものすごいがんばって最適化やら効率化をいろいろしてるのに
BLOBだと殆どのDBではそういう恩恵が受けられないこと、
データ部分だけ交換ができないこと
がもったいないと思う点でしょうか。
何より、一般的なDBでの開発経験はほぼまったく生かせません。
もったいない。。

せっかく叩く側のI/Fがきちんと決まってるので
データベースとモデル化をちゃんと分離したら
システムの可搬性もスケーラビリティも選択の余地も
ぐっとあがると思うのですが、何故この状態なのか、実はそうでもないのか
もうすこし探ってみる必要がありそうです。

GoogleMapのデータモデルの妄想

GoogleMapsは皆さんご存知の事と思います。
タイルの格納の仕方だとか、スケール問題がどうとか、
そういう話は見た事があるのですが、
このタイルをどうやって生成しているのか、
地図をどうモデル化しているのかはあまり公開されていません。

解っているのは、PostGIS的何かをHadoop的何かで
やたらスケールできるように作り直したDBのようなデータを、
多言語対応したMapServer的何かを使って
庶民には想像もつかないレベルのCPUパワーで
タイル生成を行っている、という事ぐらいです。

そこで生成されたタイルから、少し考えてみようというのが今日のお題です。

GoogleMapsと一言で言っても、少なくともMapとSatelliteの
2種類のデータを使っているのはご存知の通りです。
しかし簡単な操作で、それだけではない事が解ります。

本好きなら行かずにおれない古本屋街小川町駅付近のGoogleMapsです。



これをMapモードにするとこうなります。


もうひとつ、Satelliteのまま、ラベル表示を消すとこうなります。


この、Satelliteの状態+ラベル表示は、OpenLayersを上げるまでもなく
地図表示のレイヤーにラベルのレイヤーを重ねて表示が実現されています。
重ねているレイヤーだけ取り出すとこうなります。



見やすいように背景を透明からグレーにしてみました。
実際にはこれが透過PNGになっていて、衛星画像のjpgに重ねて表示する事で
衛星+地図となっています。
当たり前すぎることですがMapで表示しているものとは明らかに異なります。


衛星画像に重ねている地図画像を、Mapsの地図と比べてみると、
地図画像は歩道など一つ一つのディティールが角度まで書かれているのに対して、
重ねている道路レイヤーは道路が始点と終点、太さのみで書かれていて
おそらくは道路地図のベクタから生成されている物と考えられます。

これより、GoogleMapsでは
* 衛星画像のラスタレイヤ
* 地図画像のラスタレイヤ
* 道路と、鉄道と、地名表示用のベクタレイヤ
をそれぞれ持っている事が何となくわかります。

次に、地図を縮小してみます。地図表示にして、

zoom=19では


*細かい建物名の表示 + 地図記号（銀行など）

zoom=18では
*バス停表示
*大きな建物名表示
*小さな建物の輪郭

zoom=17
*中くらいの建物の輪郭

zoom=16



* 大きな建物表示

などなど、段階があるのはわざわざ見ないでもご存知と思いますが
大切なのは、ズームレベルごとに見える物が変わっていくという事は
それぞれの地物は別々の地物としてshapeを持っていて
決して一枚のラスター画像を縮小コピーしているのではない、というところです。
それぞれの地物は表示すべきズームレベルをあらかじめ持たせていて、
タイル生成の際に
zoom=16なら中くらいの建物もレンダリング、
zoom=17ならバス停もレンダリング、
zoom=18なら小さな建物名もレンダリング、

というように、タイルに書くべき物を選択していると考えるのが自然です。


同じ建物でも、zoomレベルによって表示するもの、しない物があったり、
道路もzoom5〜zoom17まで
各レベルで書くもの書かない物を細かく調整しており、
少なくとも12段階の表示レベル的指標を各道路が持っているようです。
ラスタも、ラスタではなくベクタかラスタ集合で持っていて
差分更新的何かで更新していっている物と考えられます。


...すいません、元気なくなってきました。
早足でまとめると、

現状、どのズームレベルでレンダリングすべきか、は
Mapserverやosmarender、mapnickが決定していて
タグなど(highway=primary)で決定されますが
そういった論理的なタグと、レンダリングすべきレベルは別に持っているらしい点が一つ。

もうひとつ、論理的な道路の開始/終了と、その経路のベクタと、
各交差点の形状などのラスタは別のレイヤにあるということ。
は見習うべきだと思いました。
どちらも現状ではOSMのデータも、論理構造に終止するPostGISにもなさげです。

きれいにまとまってませんが、今日はこの辺で。

タイルとlat lonの関係

ちゃんとスクリプトそのうち作りたいけど
今やってる方法

1. タイルの位置を知りたい大体の場所を http://osm.org/ で広めに開く。
2. エクポートをクリック
3. 左側に出るカラムにある、別の領域を指定するをクリック
4. 右側の地図上で範囲を選択する
5. 左の枠にlat/lonが表示される
上 maxlat 右maxlon 左minlon 下minlon
6. http://home.provide.net/~bratliff/tilecalc/を開いて、lat/lon/zoomを入力。
ここで入力できるのはbboxではなくて一枚のみです。
osmのサイトで言うと、lat=maxlat（上）、lon=minlon（左） がタイルの左上頂点、
lat=maxlat、lon=minlon がタイルの左上頂点、(wms2tile開始点)
lat=maxlon、lon=minlon がタイルの右下頂点(wms2tile終了点)で二回やる必要があります。
タイルでは、北極点が0の番号の振り方になっています。
タイルとlatlonはmax/minが逆になりますので注意してください。
TMSのzoom17はここでは4と入力
7. calcをクリックして
googleの欄に表示されるタイルのXとY、
こんなURL http://mt0.google.com/vt/lyrs=m@150&x=112882&y=53366&z=17
のxとyがタイルになります。

wms2tile

wms2tileの使い方が全然なくてソース読まないと動きも意味不明だったので書いてみる。
多分にOSMでのオルソキャッシュの情報に偏っています。

wms2tileはwmsで配信されているデータをラスタ、ベクタ関係なく
タイルを作成する物です。
作成したタイルはタイルキャッシュやTMSのサーバで配信して
Openlayersやその他TMSを使えるもので利用可能です。

基本的には
1. config.phpを編集 (config.php.sampleをrename)
2. tile.json を編集 (tile.json.sampleをrename)
3. コマンドラインPHPで動作(php range.php xxxxxxxx -)

で動きます。
2010年3月の時点ではDLしたtile.jsonの設定は
オルソ化空中写真ダウンロードシステムWMSのものになっていて、
変更のは必要ありません。コピーすればそのまま動きます。

config.phpは結果を集計するサーバの設定が必要ですが
わからないときにはとりあえずlocalhostにしておけばキャッシュはできます。

'tracker_domain'=>'localhost',
'tracker_upload'=>'',
'tracker_password'=>'password',
'user_id'=>3,

こんなかんじで。
pngcrushはパスを指定してください。

連続してタイルをダウンロードするのはrange.phpの動作になります。
PHPが上手に設定できていないとconfig.phpの
'tile_local_base'=>
にフルパス指定が必要な場合があります。

また、PHP5以降ではfunction.phpの最後の方にある
check_dir中の

  if(!file_exists($config['tile_local_base']."$layerName/$z/$x")){
-    @mkdir($config['tile_local_base']."$layerName/$z");
-    @chmod($config['tile_local_base']."$layerName/$z",0777);
-    @mkdir($config['tile_local_base']."$layerName/$z/$x");
+    @mkdir($config['tile_local_base']."$layerName/$z/$x",0777,true);
    @chmod($config['tile_local_base']."$layerName/$z/$x",0777);
  }
}

しないとディレクトリの生成に失敗してしまうことがあります。

動作には、wms2tile直下の Tilesのアクセス権を666など全ユーザ書き込み可能な状態にする必要があります。

range.phpの引数は
php range.php レイヤ名 zoom minX minY maxX maxY

オルソ化空中写真ダウンロードシステムWMSではレイヤ名を
MLITJにしているそうです。
zoomは17固定
範囲は割り当て範囲となります。

範囲指定について

基本的にはbboxと同じ考え方でタイルの頂点座標を指定します。
なぜこれが必要かというと、帯域が非常に狭いからで
DLする環境によって4〜16、場合によってはもっと多くのDLを並行した方が
効率よく回収できるためです。

割り当て範囲が(114240,52064) - (115040,52704)
の場合
php range.php MLITJ 17 114240 52064 115040 52704
で全領域となります。

これを複数に分けるにはbbox的に
php range.php MLITJ 17 114240 52064 114450 52704
php range.php MLITJ 17 114450 52064 115040 52704
とすることで分割ができます。
タイル一枚だけキャッシュするには縦横+1の数値指定で行えます。
php range.php MLITJ 17 114240 52064 114241 52065

環境によって効率のいい分割数は大きく異なります。
負荷的にmaxの位置を見極めて、
その8割〜9割程度のDL速度を保てる位置が効率よいようです。
しかし一概には言えませんので
狭い範囲で少し試してから8とか10など
適当な数でやってみるとよいと思います。


わたしははじめはざっくり大きめに実行し、
php range.php MLITJ 17 114240 52064 114460 52704
php range.php MLITJ 17 114250 52064 114470 52704
php range.php MLITJ 17 114260 52064 114480 52704
php range.php MLITJ 17 114270 52064 114490 52704
と10単位で大体作成していくつか実行、
5、2、1単位を状況に応じて追加していっています。

range.phpは数が少ない方から順にDLしていきます。
同じ箇所を追っかけで動作させるとおかしなことになるので
前からくるスクリプトに追い抜かれないよう調整する必要があります。
それか、そういう小細工をしないかどちらかです。

キャッシュファイル

wms2tileはTiles以下に

Tiles/レイヤ名/x/y
という階層を切ってファイルを保存します。
できるファイルはpngですが拡張子はつきません。

そのため、上記のようにx方向での切断の方がフォルダ単位作業的には楽です。
y方向で切断すると確認がとても面倒なことになります。

エラー

いろいろエラーが発生しますが、
取得できているところの再取得はされないので、
エラーが発生したときには全く同じ領域で再度実行すればOKです。
キャッシュできていないところだけ再度取りにいってくれます。

エラーのケース
* WMSからタイル取得失敗
* pngcrushやOSでのファイルオープン失敗
* png変換失敗
* ファイル名変更失敗
再度実行後、
エラー無くdoneの表示がでたらその領域は終了です。

報告するサーバを設定しないとreportエラーは永遠に出続けます（が、キャッシュはできています）