ラベル GIS の投稿を表示しています。 すべての投稿を表示
ラベル GIS の投稿を表示しています。 すべての投稿を表示

2018/11/05

gdalでnetCDFをFITSに変換する

GEBCOを紹介したら,天文な方からFITSで無い?と言われたので,変換してみた.

地震津波の分野では,NetCDFがよく使われています.もちろん,リモセンの時にも見てたので色々な分野で使われていると思います.一方,MacPortsのデフォルトのgdalではサポートしていないので,variants で指定する必要があります.

1年前にgdalでnetCDFを使うために,以下のようにしていました.
port install gdal +expat +netcdf +geos +hdf4 +hdf5 +xerces

さらにFITS も扱えるようにするため,さらにproj4を忘れていたようなので,以下のように作り直し(正確には新しいMacBookで新規インストール)ました.
port install gdal +expat +netcdf +geos +hdf4 +hdf5 +xerces +cfitsio +proj4

% gdalinfo GEBCO_2014_2D.nc | head -1
Driver: netCDF/Network Common Data Format
こんな感じで読めたのを確認しました.

% gdal_translate GEBCO_2014_2D.nc -of FITS GEBCO_2014_2D.fits
で変換.どうやって確認しようかなと思ったら,MacPotsに ds9 が入ってました!これは驚き.

ds9 @7.6b8_1 (science, graphics)
    SAOImage DS9 astronomical imaging and visualization application

% foreach i ({0..23})
dd if=GEBCO_2014_2D.fits iseek=$i bs=80 count=1 2>/dev/null && echo
end
SIMPLE  =                    T / file does conform to FITS standard
BITPIX  =                   16 / number of bits per data pixel
NAXIS   =                    2 / number of data axes
NAXIS1  =                43200 / length of data axis 1
NAXIS2  =                21600 / length of data axis 2
EXTEND  =                    T / FITS dataset may contain extensions
COMMENT   FITS (Flexible Image Transport System) format is defined in 'Astronomy
COMMENT   and Astrophysics', volume 376, page 359; bibcode: 2001A&A...376..359H
LONGSTRN= 'OGIP 1.0'           / The HEASARC Long String Convention may be used.
COMMENT   This FITS file may contain long string keyword values that are
COMMENT   continued over multiple keywords.  The HEASARC convention uses the &
COMMENT   character at the end of each substring which is then continued
COMMENT   on the next keyword which has the name CONTINUE.
HIERARCH LAT#AXIS = 'Y       '
HIERARCH LON#AXIS = 'X       '
END

ヘッダも出来てた.

が,,,180度回転してるな...

ImageMagickのconvert で -rotate 180 で回転は出来たものの,出力されたFITSは正の値しか持ってなかった.一度値を変換して,ほげほげすれば出来そうな気もするけど,一旦おしまい.

(追記 11.8)
ファイルがでか過ぎてマカリで開けなかったという連絡が来たので,半分に解像度を落とす.
 % grdsample  GEBCO_2014_2D.nc -GGEBCO_2014_2D_1min.nc -I1m
 % gdal_translate GEBCO_2014_2D_1min.nc -of FITS GEBCO_2014_2D_1min.fits

これで,サイズは1/4と思ったら,Float32になっちゃってて半分にしかならなかった...


2014/08/29

Projで座標変換

ある1点の緯度経度からUTMに変換する方法.

proj.4 に含まれている cs2cs で出来る

% echo "140.0 36.0" | cs2cs +init=EPSG:4326 +to +init=EPSG:32654 -f %.4f
409870.9506     3984410.7881 0.0000

入力データの座標系(測地系+投影図法)は,EPSG:4326でWGS84の緯度経度座標.
出力データの座標系は,EPSG:32654でWGS84のUTM 54N.

EPSGコードは,proj.4 パッケージの epsg ファイルで確認できる

% grep -B 1 32654 /opt/local/share/proj/epsg
# WGS 84 / UTM zone 54N
<32654> +proj=utm +zone=54 +datum=WGS84 +units=m +no_defs  <>

ちなみに,ここ(Google Maps)

参考:
http://ryuiki.agbi.tsukuba.ac.jp/~nishida/lecture/GIS/proj.html

余談:
Macの場合,Ports に入っているのでついでによく使うgdalも入れちゃうのが良い.

% port install gdal

で,projも依存関係で入る.

2014/01/10

タイルサービス雑感

タイルサービスは似てるのが乱立してしまっていて以下の 3つある.
  • WMTS (Web Map Tile Service, OGC)
  • TMS (Tile Map Service, OSGeo)
  • WMS-C (Web Mapping Service - Cached, OSGeo)

で,そのタイルサービスを構築したいのだが,WMSもやりたいので手を抜きたい.というわけで,WMSを裏でオンデマンドで呼んでタイルサービスをやる.にあたってのソフトウェア比較.

  1. MapProxy
    • WMTS, TMS, WMS-C 全部に対応している
    • WMTS を使うか TMS を使うか選択する
      • 両方のサービスを提供するならサーバを2種類用意する必要がある
    • GeoWebCache より動作が軽い
    • WMS としても使える
    • (タイルを切り貼りして拡大縮小して返すこともできる)
  2. GeoWebCache
    • WMTS, TMS, WMS-C 全部に対応している
    • WMS としても使える
  3. TileCahce
    • 少し古い

MapProxyとWMTSを使うことにした.
タイルは EPSG:4326(緯度経度) と EPSG:900913(メルカトル図法, Google)
の両方を用意することにした.

GSJのサービス構築の時の検討内容でした.ここのWMSが基本的には全てWMTSで提供されています.


2013/11/12

Google Mapsのマイマップはどこへ?

久しぶりにGoogle Mapsを開いて,マイマップを編集しようと思ったら,新インタフェースにそれらしきリンクが見当たらない...と思ったら,ここにあるように歯車の中の「マイプレイス」だそうだ.これじゃ見つけられんよ....


2013/11/02

EasyWMSView version 1.3公開

FOSS4G 2013 Tokyo の発表に合わせて,WMTSとGetFeatureInfoに対応したEasyWMSView version 1.3.1を公開しました.次のバージョンでは位置情報付きPDF出力(印刷)サービスを公開する予定です.


2013/10/31

GeoSNSNotifier 1.0公開

産総研のオープンラボで講演してきました.その中でRSSから更新情報をTwitterやFacebook等に投稿するツールであるGeoSNSNotifierの紹介もしたので,合わせてVersion 1.0 を公開しました!


2013/08/05

google = 900913?

もはや正式扱いになりつつある EPSG:900913 ですが,projにはまだ入らないしどんなもんなんですかねぇ...

https://www.google.co.jp/search?q=google+epsg

proj のepsgリストが一応正式なのか?現実問題としては900913を扱えるようにするためにprojにあるepsgリストに追記したりするわけですが.

一つ上のディレクトリにはたくさんおまけが置いてある.

2013/05/23

新GSJ DBスタート?

昨年12月にFOSS4G Tokyo 2012でお話させて頂いた産総研地質調査総合センター(GSJ)のデータ配信が始まりました.せっかくなので,地質の日である5月10日に合わせてリリースしています.ただ,ライセンスについては検討中のため改定出来ていません.今後より使いやすいライセンスに書き換えていく予定です.

WMSはEasyWMSViewを使って見えるようにしてみました.

2012/12/26

WMS image/jpeg or image/png

JPEGの方がレンダリングが速いのは自明と言ってもよいのだが,透明(アルファチャンネル)が欲しい場合はPNGかGIFにせざるを得ない.一方で,basemapとして利用する場合は透明じゃなくてもよいのでJPEGを選択できる.が,OpenLayersで,

transparent: true, format: 'image/jpeg'

と設定してしまうと,transparentのtrueが優先されて,image/pngでリクエストが送られてた.

OpenLayers 2.12で確認

2012/07/06

GeoRSSやってみた

ひさしぶりの連合大会で発表したASTER高温領域検出システムで試してるのが,GeoRSSです.何が出来るか?というと,RSS(ホームページの更新情報フィードですよね)に位置情報をつけて,Geo化してしまおうというもの.つまり,「そこ」で新しい何かが起きたよということを通知できるのだろう.

んで,見ると緯度経度をともかく追加すれば完成できるみたい.
<geo:lat>緯度</geo:lat>
<geo:long>経度</geo:long>
とするか,
<georss:point>緯度 経度</georss:point>
でいいらしいので,

<geo:lat>-9.505607461</geo:lat>
<geo:long>-42.8624929092</geo:long>
と入れてみたのがこちら


さらに,Google MapsにURLそのものを渡してあげると勝手に描画までやってくれる!
という驚き.

Google MapsにGeoRSS URLを与える


2012/06/29

Google Earthを使ってWMSを見る

WMS (Web Map Service)を使えば色々なクライアントで勝手に閲覧ポータルを作ることも出来るからご自由にドウゾ的なところはあるわけです.

が,単に見たいだけの人にとって(かつ,GISじゃない人)そんなこと言われても何のことだか分かりませんね.というわけで,まぁまぁ汎用なGoogle EarthでWMSを見るための手順を書いてみました.

今回対象にしているのは,産総研GEO Girdが出しているベースマップのWMSで,URLは,
です.まず,このリンクをクリックしたら,
"No query information to decode. QUERY_STRING is set, but empty."
とか言われておしまいですね.何もQueryが無いから何も答えられないと言っている(この場合は,MapServer)わけです.じゃ,WMSとして正しくやると, 今度は何か返って来ました.が,これを読んでどうしろ?と思いますね.じゃ,ここから先はGoogle Earthにお任せで.

まずは,イメージオーバーレイの追加を選択します.普通は分からないところで,【更新】メニューをみると,【WMSパラメータ】を設定できるようになっています.
イメージオーバーレイの追加の更新メニュー
【WMSパラメータ】ボタンを押すと,WMSのURLを選択できるようになっています.ここでは,自分が知っている他のWMSを貼り付けるので,【追加】ボタンをクリックします.
WMS URL追加メニュー
【追加】ボタンをクリックすると,URLを入力するウインドウが開くので,そこに,
http://ows.geogrid.org/basemap
を入力します.
WMS URL入力画面
【OK】を押すと,画面左側に何かがリストされます.実はこの時,上のGetCapabilitiesリクエストがサーバに送信されていて,よくわからないXMLが解釈されて左側のリストが出来上がっています.1行毎のエントリをレイヤと呼びます.このうちの1つ,Land Cover Map MOD12を選択して,【追加 ->】ボタンをクリックします.
レイヤーの追加画面
これで,WMS画像の貼り付けは完了なのですが,実際には見た目が汚いため,少し調整が必要です.まずは,画像のフォーマットがgifになっているので,jpegもしくはpngにするほうが良いです.特別な画面があるわけではないので,【リンク:】の中で直接修正する必要があります.
画像フォマット変更画面
もう一つは画像のサイズです.デフォルトでは,512 x 512ピクセルの画像がリクエストされています.これだと,小さい画像を拡大表示してGoogle Earthの上にレンダリングしてしまうので,非常にぼけた画像になってしまいます.ディスプレイのサイズに合わせて,
WIDTH=2048&HEIGHT=2048
などとして下さい.ここでの上限値(サーバの設定)は,2048です(いずれ大きくなるかもしれません).
リクエスト画像サイズ変更画面
最後に【OK】をクリックすると,表示されます.
Google EarthでMOD12プロダクト(土地被覆図)WMSを貼りつけた結果
(補足)
  • キャプチャの取得には,Mac OS X版Google Earthバージョン6.2.2.6613を使いました 
  • 画像サイズである,WIDTHとHEIGHTは,GoogleのWMSクライアントが画面のサイズじ応じてリクエストしてくれても良いはずと思ってるのですが... 
  • KMLの中でWMSを定義して提供することも可能です.地震動マップ即時推定システムで,各地震毎にKMLでくるんだWMSを提供しています.


2012/06/11

火山データベース

火山の「あたり」を定義する情報が欲しい.

スミソニアンが整備しているリストが一番大きいようで,情報がたくさんある.
一番単純なのはエクセルみたい.だけど,かっこよくない.が,KMLがある!しかもon-demand updateが出来るようになっていて,7日間で更新と定義されてる.やるなぁ〜.

参考


(思ったよりネタが少なかった)

on-demand updateなKMLを書くのはとても重要.普通にKMLに情報を書いて渡すのは良いけど,その後更新ができなくなってしまいます.これを回避するためには,小さいKMLを渡して,実体のデータはサーバ(CGIとか)で常時配信できるようにしておくのが良いです.

      <Link>
        <href>http://kml.example.com/volcano.pl</href>
        <refreshInterval>2</refreshInterval>
        <viewRefreshMode>onStop</viewRefreshMode>
        <viewRefreshTime>1</viewRefreshTime>
      </Link>

という感じでもろもろ設定できます.

2012/06/07

geometryカラムに対するPerl DBIからのinsert

geometryカラム(ここでは,the_geom)に対して,
     my $geom = "st_setsrid(st_makepoint($lon,$lat),4326)";
     my $sth = $dbh->prepare("insert into events (the_geom) values (?)");
     $sth->execute($geom);
すると,


     DBD::Pg::st execute failed: ERROR:  parse error - invalid geometry
     HINT:  You must specify a valid OGC WKT geometry type such as POINT, LINESTRING or POLYGON
とか言われてしまう.


が,
     my $sth = $dbh->prepare("insert into events (the_geom) values ($geom)");

     $sth->execute();
として,executeの引数じゃなくしちゃえば動く...

そんな話なのか....


あ,でもこれじゃ引数変えて実行出来ないじゃん.


2012/06/05

Geospatial PDF (GeoPDF)

GeoPDFって,OGCで話題にはなっていたので知ってはいたけど全然状況を把握してなかった....

日立ソリューションズのTerraGo Product Suitesで作ったり見たり出来るのね.で,作るソフトは有料だけど,見るのは無償ですよと.
(だけど,そのリンクが, .exe 直リンだったりするので窓な匂いしかしないからやめる)

が,同じページにサンプルが幾つか転がってたので見てみた.数値地図ってのがあるので,筑波山を見ようと思ったけど,水戸で開いたPDFがクリッカブルなようでクリックできない.ん?と思ったらChromeのPDFエンジンがダメなだけのようだ.Safariで開き直したら,筑波山が表示できた!

Safariで筑波山を拡大

さらに,よく見れば左側のメニューでレイヤー選択が出来るようになっている.そんなことが出来るのね.もう一つ東京駅周辺GeoPDFというのを見てみた.が,今度はSafariで見ても普通のPDFを見ているのと何も変わらない.それはおかしいと思ってPDFをダウンロードして,Adobe Readerで開いてみると,右上の「拡張」の中に「地図位置ツール」ってのが隠れていた!
Adobe Readerの地図位置ツールを使って座標を表示

まぁ,だからなんだというのはおいといて,座標はとれた.ものさしツールを使えば,直線距離,周囲距離,面積を測ることも出来た.なるほど.

最後にオバマ大統領就任パレードの警備配置に実際に使われたらしいGeoPDFを見る.これもダウンロードしてAdobe Readerで見た.たぶんたくさん色々な機能があるんだろうけど,ピンをクリックしたら,写真が出てきた.

ピンをダブルクリックしたらプレビューで写真が開いた
ただ,この写真をオンラインで持ってきて表示したのか,PDFに含まれていたのかはよくわからない.Tmpディレクトリにjpegファイルは出来ていた.あと,こんな警告は出るのでどうしたものか.....



そもそも,GeoPDFという言葉はOGC(Carlが喋ってた)では聞いたけど,標準規格リストに無い.つうことは,OGCの標準ではないはず(?)だ.ということは何か.
でGoogle様に聞いてみると,08-139でBest Practiceっぽいということが分かる.

ここまでくれば簡単で,OGCのサイトでBest Practiceを見ると確かに08-139r3が出ていた.タイトルは,
になってました.とすると,仕様は?と調べを進めると,WikipediaにGeoPDFがあった.GeoPDFは規格じゃなくてTerraGoのフォーマットか.さらに,Adobe's proposed geospatial extensions to ISO 32000ということは,OGCじゃなくてISOか?もう一つ,Geospatial PDFというリンクがあるのでたどってみた.で,こっちの方に明確にThe georeferencing metadata for geospatial PDF is most commonly encoded in one of two ways: the OGC best practice and as Adobe's proposed geospatial extensions to ISO 32000となっていた.

つわけで,タイトルも変える...


つか,ISO 32000は,PDFそのものだったのか.proposed geospatial extensions to ISO 32000なわけだから,PDFの規格そのものとして georeferencing extension を追加しようという試みだけど,まだ決まりきってはいないのかな?

やっぱり,GeoPDFはTerraGO Technologies社の登録商標って書いてあった(Iページ).


ブログで追記って,反則な気もしないでもないけど....Adobe Readerのバージョンについて補足.

上記のスナップショットは,バージョン 10.1.3 です(見て分かると思いますが,Mac OSです.10.6).まだ9.4.6のままになっていたMac Book Airがあったのでそっちで試してみましたが,ツール・注釈・拡張という右上のメニューが出ませんでしたが(レイヤーのON/OFFは表示されて操作も出来る),「ツール」→「分析」に隠れています.「分析ツールバーの表示」でも出てきます.また,数値地図GeoPDF利用ガイドによると,9.5.0で動作確認済とありますが,9.4.6でも問題はなさそう.

Adobe Reader 9.4.6 (Mac)での地図ツールの場所

さらに,Adobeのバージョン9.1のリリースノートにGeospatial Locationの話が出てるからこのタイミングくらいから使えるようになってたのか?

2012/05/27

まっぱーもどき

某所で工事中(新規)だった道路が開通していたので,GPSロガー (myTracks)をOnにして車で走ってきた.これ使えば,OSMで道路繋げられるの?どんくらい簡単にできるものかやってみようかと思いつつ,時間が無い!


新しい道路があったからGPS logger onにして走ってきた。 #OSMjp... on Twitpic

いずれ....やろうかな....(と言ってるうちに繋がってそう).関連

2012/05/23

GPX

myTracksというiPhoneで動くGPS ロガーを使って少し歩いてみた.

実は前にも試したんだけど,Macに取り込むには,有償のマック版を入れる必要があるということで,放置してたのだが,メールでデータを送信というのがあったので改めて試してみた.

送られてきたデータは,GPX形式のファイル.ともかく見てみたかったのでで,TrailRunnerで表示.もうちょっと目的があったらまじめに考えるかな.


なんと....見るだけなら,Google Earthで読み込んでしまえば良いだけだった...


久しぶりの連合大会

連合大会に久しぶりに行って来ました.正式名称は日本地球惑星科学連合2012年度連合大会.って長すぎ...

国際標準規格を用いたASTER高温領域検出システムの開発
View more presentations from Naotaka Yamamoto.

もうちょっとで公開できるはず...

iPad Google Earthで見るとこんな感じ
iPad Google Earthでみるとこんな感じ
http://ekbo.blogspot.jp/2012/05/bl... on Twitpic

2012/05/17

MapServer 2048

WebGISサーバとしては,一番お手軽でラスタに対しては最速(と思う)MapServerですが,余り大きい画像が作れません.

msWMSLoadGetMapParams(): WMS server error. Image size out of range, WIDTH and HEIGHT must be between 1 and 2048 pixels.

というエラーが表示されます.

これを回避する場合は,マップファイルのMAP
MAXSIZE 4096
とかすると,大きい画像も返せるようになります.

当然,処理に時間がかかるようになるのでサーバには負担がかかるので注意っす.


2011/12/14

ASTER TIR輝度温度

ASTERのTIRの輝度温度を求める方法を整理.

TIR (Band 10-14) で,DN値からとりあえず大体の温度を知りたいので輝度温度Tを求めるのだが,まずは,放射輝度Rを

R = (DN - 1) * UCC [W/m^2/str/um]
で求める.UCC (Unit Conversion Coefficients)は,バンド毎に求められている単位変換係数.実際の値はプロダクトに含まれるが,
Band 10 ... 0.006882
Band 11 ... 0.006780
Band 12 ... 0.006590
Band 13 ... 0.005693
Band 14 ... 0.005225

となっている.なお,DN=0は,データ欠損のため,DN-1をする時か,その前にDN=0を何らかの形で除外する必要がある.

放射輝度Rを用いて,輝度温度Tは,
T=K2/ln(K1/R+1)
で求めることが出来る.この時,K1とK2は,
K1 = 2hc^2/λ^5 [W/m^2/str/um]
K2 = hc/kλ [K]
であり,実際の値ははバンド毎に異なり(波長が異なるため),
Band 10 (λ=8.300) ... K1=3.00e3, K2=1.73e3
Band 11 (λ=8.650) ... K1=2.44e3, K2=1.66e3
Band 12 (λ=9.100) ... K1=1.90e3, K2=1.58e3
Band 13 (λ=10.60) ... K1=8.84e2, K2=1.35e3
Band 14 (λ=11.65) ... K1=6.42e2, K2=1.27e3
となっている.

例えば,Band 13で,DN=2640の場合,
R = ( 2640 - 1 ) * 0.005693 = 15.023827 [W/m^2/str/um]
T = 1350 / ln( 884 / 15.023827 + 1) = 329.938459616334 [K]
= 329.938459616334 - 273.15 = 56.7884596163339 [C]
となる.有効数字無視しちゃってるけどご容赦.

出展)




2011/11/07

FOSS4G 2011 Tokyo

GIS学会からやっと日本のFOSS4G / OSGeo / OSMの世界に飛び込んだ感じ.これまでFOSS4G本体には2回ほど参加してたけど,国内は初めてでした.なにせ,土日というのはなかなか参加の踏ん切りがつかない....

GEO Gridの紹介ということで少しまとめたスライド.

もうちょっと面白おかしく&未来話しをするつもりだったんだけど明らかに準備不足.

おまけ:
FOSS4G Tシャツなのだ