マーカ等のマウスホイールイベントが効かないのは最新版では解決されていました

openlayers

id:yellow_73:20120909 の続き。
当時はできなかったでいいと思うのですが、OpenLayers 2.13 からは、マーカなりの要素のCSSクラスに "olScrollable" を付けると、マウスホイールをパスしてくれるようになっています。…知らんかった。

CSSクラス名は HTML なら class属性値 で指定し、DOMなら className プロパティ名で指定しますね。

ogrでチャカチャカっとベクタデータを作成する

GIS

たとえば、ラスタデータの地理範囲が数字では分かってるけど、ベクタデータになっていなくて、GISのレイヤにできない、そういうことがたまにあるかも知れません。

…現に今直面してました。

WKT

数字として分かっている場合には WKT を使うとコピペとちょっと飾りをつけてやるだけで表現できます。

ポイント、ラインストリング、ポリゴンのWKTの例はそれぞれ次のようになります。

POINT(135 35)
LINESTRING(135 35, 135 36, 136 36, 136 35, 135 35)
POLYGON( (135 35, 135 36, 136 36, 136 35, 135 35) )

注意すべき点は次の通りです。

  • 点ごとの区切りはコンマで、座標値ごとの区切りは空白 (座標値もコンマ区切りをするとパーサがエラーと認識する)
  • 座標値は 東西、南北の順に固定されている (EPSGエントリに依存しない)

CSVデータの作成

CSVデータを作成します。属性値とジオメトリとをまぜて書いちゃって下さい。

  • 1行目はフィールド名
  • WKTはコンマを使うので、WKTの箇所はダブルクォートで括る

たとえば次のようになります。

code,geom
11,"LINESTRING(135 35, 136 35, 136 36, 135 36, 135 35)"
12,"LINESTRING(136 35, 132 35, 132 36, 136 36, 136 35)"
11,"LINESTRING(135 36, 136 36, 136 37, 135 37, 135 36)"
12,"LINESTRING(136 36, 132 36, 132 37, 136 37, 136 36)"

VRTファイルの作成

意味づけはVRTファイルで行います。基本的にはVRTがメインで、CSVはVRTからのリンクとなります。
詳細は、本家サイトの http://www.gdal.org/ogr/drv_vrt.html あたりを参照してください。

なお、さきほど作成したCSVファイルは "frame.csv" としています。


    
        frame.csv
        
        wkbLineString
  • GeometryTypeの値に "wkb" とついているのは、とりあえずそういうものと考えることにしています。
  • Field要素を入れているのは、ほうっておくとgeom要素をジオメトリフィールドと文字列フィールドとを兼用するからです (要調査)

LineStringでなくPointとする場合には、GeometryTypeの値を変えて下さい。

確認する

ogrinfo -ro -al frame.vrt

でファイルがどのようにogrに認識されているかが分かります。

INFO: Open of `frame.vrt'
      using driver `VRT' successful.
ERROR 1: Failed to open datasource `frame.csv'.

Layer name: frame
Geometry: Line String
Feature Count: 0
Layer SRS WKT:
(unknown)

>ogrinfo -ro -al frame.vrt
INFO: Open of `frame.vrt'
      using driver `VRT' successful.

Layer name: frame
Geometry: Line String
Feature Count: 4
Extent: (132.000000, 35.000000) - (136.000000, 37.000000)
Layer SRS WKT:
(unknown)
code: String (0.0)
OGRFeature(frame):1
  code (String) = 11
  LINESTRING (135 35,136 35,136 36,135 36,135 35)

OGRFeature(frame):2
  code (String) = 12
  LINESTRING (136 35,132 35,132 36,136 36,136 35)

OGRFeature(frame):3
  code (String) = 11
  LINESTRING (135 36,136 36,136 37,135 37,135 36)

OGRFeature(frame):4
  code (String) = 12
  LINESTRING (136 36,132 36,132 37,136 37,136 36)

なお、ジオメトリフィールドが表示されないまたは文字列フィールドと認識されている( "geom(String) =" が付く)場合には、パースエラーが起きています。
パースエラーのエラーメッセージは出ませんので、おかしいと思ったら WKT が正しいかを確認して下さい。

シェープファイルにする

md frame
ogr2ogr -f "ESRI Shapefile" frame frame.vrt

ここでは、frameフォルダをつくり、そこにシェープファイルを入れています。

gdalをWindowsでビルド

GIS Windows

標題通り。計算用の機械が他の仕事で動きが怪しくなっているため、Windows機でも仕事をさせたい。Windowsならosgeo4wを入れると早いのですが、osgeo4wのGDALは32ビットバイナリなので GDAL_CACHEMAX に2Gのカベがあります。そこで、64ビットバイナリを作りたい、と。

VC持ってたら、特に問題なくコンパイルできました。
http://trac.osgeo.org/gdal/wiki/BuildingOnWindows を参照。
ソースを展開したら、そのツリーのルートに移動。上述WikiのBasic Options節を参照しつつ、nmake.opt をゴニョゴニョする。

それから次を実行(64ビットバイナリの場合)。

call "C:\Program Files (x86)\Microsoft Visual Studio 10.0\VC\vcvarsall.bat" amd64

これによって、ローカルのコンパイル環境を反映してくれる…でよかったのかな。少なくともこれを実行しないと、パスを通しててもコンパイルできなかったりします。

インストールまで終えると、インストール先ディレクトリに bin, data, html ができます。

直感的に (gdal)\bin に Path を通す必要があるのは分かると思いますが、もうひとつ。

GDAL_DATA を (gdal)\data に設定します。そうしないと、EPSGコードリスト等のデータファイルが読めずに、GDAL_DATAを設定しろエラーが発生します。

FGD対応ogrを更新してくださいました

GIS 雑談

http://www.osgeo.jp/foss4g-mext/ から id:yellow_73:20130527 で言及した問題をクリアして下さったものが出ています。バイナリ版については確認しました。

1.9の問題箇所だけならすぐ話が終わったかと思いますが、バージョンを上げたものでの対応だったので、これは対応早かったと感じました。

今となっては、実はここまで書いてる mext_gdal を実際に使ってるのは DEM だけ(今回の問題とは直接関係無い)で、道路とか建物のベクタデータをmext_gdalで扱ってはいない、というのは言えない…。

float型とdouble型

雑談

C言語では、float型はよく考えたうえで使うべきで、基本的にはdouble型にしたほうが良いと思いますという話。

数値演算ユニットが入ってないとか、メモリや外部記憶が少ない機械で使うとかいった場合には、float型を用いる方が良い場合がかなりあると思います。

しかしながら、C言語の場合は、float型とdouble型が混ざる演算ではfloat型の値をいったんdouble型の値にしてから演算するので、逆に余計なステップを要します。
たとえば http://www.pro.or.jp/~fuji/mybooks/cdiag/cdiag.4.4.html 等参照。もちろん"0.1"を"0.1F"にして単精度定数にすると、全てfloat型の演算になるので、printfに渡すところを除いては、余計なステップは不要です。

また、引数型が定まっていない引数にfloat型の値を「値渡し」で渡す場合にはdouble型に変換したうえで渡します。fprintf系等を呼ぶ場合には余計なステップが発生します。

なので、float型でも十分な精度であって、アーキテクチャ上float型の方が高速である場合であって、かつfloat型とdouble型とが混じらないように慎重にコーディングするなら、float型を使うべきですが、そうでないならdouble型を使うべきです。

つまり、今のパソコンなら、ほとんどdouble型にした方が良いと思います。

余談 #1

double型にした方が良いと*思います*、と自信が無い口ぶりなのは、GPUとかは触ってないためです。もしかしたら何か間抜けなところがあるかも知れないと戦々恐々としている模様。

余談 #2

fprintf系ではfloatがdoubleに変換されるため、かつては"%lf"が無かったんです。現在は"%lf"と"%f"とは同等です*1。なお、fscanf系は「ポインタ渡し」なので、"%f"と"%lf"の区別がしっかりあります。

余談 #3

Z80は浮動小数点数演算どころか整数の乗除算さえなかった。ていうか今気づいたけど単精度の仮数部がレジスタペアに入らんのでないか…。

*1:JIS X3010:2003によると(p.199)「l(エル) … a, A, e, E, f, F, g又はG変換指定子の場合，何の効果も無い。」

基盤地図情報対応ogr2ogrの出力が多少残念になる問題

GIS

2013/5/31 この問題について対応して下さっています。

id:yellow_73:20130411 で出てきた基盤地図情報対応GDAL/OGRの続き。

問題無いという問題な判断

金曜日にtwitterで流れてきた話。基盤地図情報2500の建物とか道路縁とかのベクタデータを ogr2ogr で他の書式に変換しようとした際に、地理院さんが公開しているコンバータFGDVと結構ちがってくる、というのが証拠付きで出しておられる方がいらっしゃいました。ogr2ogrの出力の方が問題がありそうな雰囲気で、ogr2ogrでの出力は「ガクガクしている」「箱型になるべきが箱型になってない」ようなかんじでした。

当初私はogr2ogrの結果と基盤地図コンバータの結果とを、重ね合わせずに見比べて、「こちらではそういう症状は全く出ない」と判断しました。FreeBSDをソースからコンパイルしたのを使っていたので、Windowsバイナリ特有の問題ではないか、と。これで無事解決。


ところが、月曜日になって、金曜日に作成したogr2ogrの出力とFGDVの出力とを自分の手元で実際に合わせてみる。


(JGD2000経度緯度 青がogr2ogr, 赤がFGDV)

なんということか、Windowsバイナリ特有の問題ではなかったのです。完全に間違ってたorz

いや意外と気付きにくかったんです、本当です信じて下さい、かなり拡大しないと分かりにくかったんです、それにウィンドウを並べて見比べるだけでは分かりにくいんです…重ねて表示させろってことですねorz

floatをdoubleにすればOK

で、ソースを見ると、経度緯度文字列をsscanfで数値型に変える際に、floatを使っていました。これをdoubleに変えれば、精度の問題なら解決するはず。

変更に先立ち、他のところでもfloatを使ってないかチェックしないといけない。出力までの間にfloatに落とされるなら、結局はdoubleでやっても意味がないです。

文字列からsscanfで読み取った値を納めている変数 fLon, fLat は、文脈に応じて m_poPoint, m_poLineString, m_poPolygon にポイントを上書きしたり追加したりします。
いずれも src/ogr/ogr_geometry.h にある OGRRawPoint を使っています。これのメンバ変数(x,y)はdouble型変数でした。
よって、floatをdoubleに変えると何かいいことがあるかもしれない、と判断。

src/ogr/ogrsf_frmts/fgd/fgdreader.cpp で
727行の

float fLon, fLat;

を

double fLon, fLat;

にします。

読み込みには sscanf() を使っていて、%f だと float型を読むので、double型だと %lf にしないといけません。
ということで同じく src/ogr/ogrsf_frmts/fgd/fgdreader.cpp で
730行の

sscanf( m_szTempChar, "%f %f", &fLat, &fLon );

を

sscanf( m_szTempChar, "%lf%lf", &fLat, &fLon );

にします。
(5/28 1:55 "%lf %lf"の真ん中の空白文字を消しました。仕様上区切りの文字列が"SPCひとつと0文字以上の空白文字"と定まってないのではないかと思った(確認してない)ため)

結果はこのようになります。

(JGD2000経度緯度 緑が訂正ogr2ogr、赤がFGDV 重なっている)

他への影響は検討していませんが、とりあえずうまくいきます。

出典

上記の画像は基盤地図情報(縮尺レベル2500)を使用しています。

gdal2tilesでEPSG:900913がズレる問題は900913を使わなければ済む

GIS

id:yellow_73:20110408#p1 (2年前)で、こっそり gdal2tiles ... -s epsg:900913 と、EPSG:900913 を指定しています。
本来はGeoTIFFにEPSG:900913を埋めているので問題ないはずですが、なぜかgdal2tilesとかMapTilerとかではうまくいかない。ちょっと北にずれる。

これ、根本的な問題解決をするには GDAL に手を加えないといけないのではないかと思われます。gdalinfoで見たらEPSG:900913だとWGS84メルカトルと解釈されるような空間参照系WKTだった、と記憶しています。

解決策が分かりました。

EPSG:900913 ではなく EPSG:3857 を使えばOK。

ああ、石投げないで。