[QGIS] ボロノイポリゴンのバッファ領域について

QGISのメニューバーの「ベクタ」→「ジオメトリツール」→「ボロノイポリゴン」で、ポイントデータからボロノイ図が出力できるのですが、その入力パラメータとして「バッファ領域」というものがあります。
「これなんなん?」と問われたので実データでチェックしました。下図を参照ください。メニューバーの「ベクタ」→「調査ツール」→「レイヤ領域のポリゴン」(Layer extentを矩形で出力する)の出力結果が赤枠、バッファ領域を10%と設定して出力したボロノイ図が緑色のポリゴンになります。
(QGIS 2.14 x Windows 8.1)

qgis2webで出力したwebmapの背景図を変更する

先日宮城と岩手で行ったQGISハンズオンで、leaflet.jsベースのウェブマップを簡単に出力できる qgis2web プラグインを紹介しました。その際、背景図をOpenStreetMap以外に変更できないか?という質問がありました。帰りがけだったため、ざざっと「ここのURLを書き換えるだけでOKです」と説明したのみでしたので、再度その内容をここに書いておきます。

index.htmlをテキストエディタで開いて35-40行目付近(地図に載せたレイヤの数により前後します)にあるこの2行を書き換えれば地理院地図の背景図に変換できます。

<変更前:OSM>

<変更後:地理院タイル(標準)>

なお地理院地図で提供されているタイル地図であれば上記の
‘http://cyberjapandata.gsi.go.jp/xyz/std/{z}/{x}/{y}.png’
に相当する箇所のURLを変えるだけで背景を変更できます。URLのリストは以下のリンク先を参照ください。

地理院タイル一覧
https://maps.gsi.go.jp/development/ichiran.html

Pythonから使いやすくなったGDAL 2.1がステキだ

これは FOSS4G Advent Calendar 2016の13日目の記事です。

今年はどんなことを書こうかな、とか思ってたら..ここまでみなさん トバしてますね…
あがりっぱなしなハードルをここでぐっと下げるために
今年(2016年)にあった各種のアップデートの中で、あ、これ便利だなー と思ったことをそのまま書きます。 以降、超地味な内容です(派手なのむり

結構前になりますが2016年の5月にGDAL 2.1.0 がリリースされました。
http://trac.osgeo.org/gdal/wiki/Release/2.1.0-News
いろいろあるアップデートのなかで、Pythonでの取り扱いが便利になったところが私にはぐっときました。

GDALに用意されているコマンドラインユーティリティ(gdal_translateとかgdalwarpなど)はとても便利です。

Pythonでこれらのコマンドと同じ処理を行いたいなーというときには

  • subprocessとしてGDALのコマンドを呼び出す
  • array(配列)として吐き出して各種の処理をする

という手続きを踏んでいました(私は)。
もちろんこれでも十分楽に扱えるのですが、
大変ものぐさなわたしは
もすこしシンプルに同等の処理を書けたらいいなー
とか思っておりました。
これが誰かに伝わったのか(そんなわけないですが
GDAL/OGR utilities as a libraryというRFCがGDAL2.1で採用されておりました。
Python bindingの項をみますと

For Python bindings, convenience wrappers are created to allow specifying options in a more user friendly way.

とあります。
まえよりわかりやすくオプションを指定できる便利なラッパーができたよー とのことですね。
それではみてみましょう。

いちばん単純なフォーマットの変換処理を例に取りますと、


というコマンドと同じ処理を、いままでPythonでこんな風に書いていました(私は)。


んで、GDAL 2.1 x Python ではどーなるかといいますと..


と、とてもシンプルに書ける様になりました。
もひとつ例を挙げるとラスタデータの投影変換はこんなコマンドですが


GDAL 2.1 x Pythonでは…


という感じです。オプションの指定がとてもわかりやすくてヨイです。

もちろんラスタデータに対して複雑な処理をするときはarray(配列)にしたほうが便利なのですが、単純なデータ処理(フォーマットの変換 投影変換 クリップなどなど)をしたいときにすっきり書けてうれしいです。

その他の処理やオプションについての詳細はこちらを参照ください。
http://gdal.org/python/

すこし古い話ではありますが私にとってありがたいアップデートだったので書いてみました。
ではでは。よいお年を(まだ早い

ラスタの領域を拡張する(QGIS)

とあるところでタイトルのとおりの処理のやり方を提案してみたところ
ぽつぽつとリアクション(そーやってもできるのか!的な)があったりしたので書いてみます。
コネタでも書いておくのだいじなのね..と思いました。

さて、とある湖岸線(琵琶湖,支笏湖などなど)を示すラスタデータがあるとします。
湖岸線を示すピクセル値が1でありその他の領域がNODATAとされている
と仮定します。

その周辺域のにおいて湖岸線からの距離を面的に求めたい場合は
QGISのラスタ→解析→プロキシミティ(ラスタ距離)というツールをつかうと便利です。
プロキシミティってなんじゃらほいって解説は以下のサイトを参照ください。

gdal_proximity.py
(QGISのプロキシミティはこのコマンドを生成して実行するためのインターフェイスです)
http://www.gdal.org/gdal_proximity.html
QGIS | 道からの距離を視覚的に表現
http://y4su.seesaa.net/article/430793300.html

さてこの処理ですが、出力ラスタの領域がもとのラスタデータと同一領域となるため
湖岸線ラスタが湖岸線が丁度収まる領域のデータであった場合
知りたい場所からの湖岸線までの距離はもとめられないことが多いのですよね。
(ex.琵琶湖の範囲ぴったりの湖岸線ラスタだと京都駅からびわ湖畔までの距離はもとめられない…)
なので湖岸線ラスタの領域を拡張したくなることがあります(拡張領域はNODATAでよい)。
その手順を以下より書きます。

  • 元のラスタの領域情報を抽出する

QGIS(今回はvesion 2.16を利用しています)のメニューバーのラスタを開き
→その他→情報を選択してください。
入力ファイルとして湖岸線ラスタデータを指定してOKしてください。
対象ラスタの情報が出力されますが
その中に以下のようなラスタデータの領域を示す
四隅の座標値が含まれています。

Upper Left ( 491925.257, 4781017.207) (140d54′ 2.30″E, 43d10’54.76″N)
Lower Left ( 491925.257, 4723777.207) (140d54′ 5.27″E, 42d39’59.12″N)
Upper Right ( 534495.257, 4781017.207) (141d25’28.06″E, 43d10’52.08″N)
Lower Right ( 534495.257, 4723777.207) (141d25’15.37″E, 42d39’56.48″N)

これをテキストエディタ等にコピーしておいてください。

  • クリッパーツールを使ってラスタの領域を拡張する

メニューバーのラスタ→抽出→クリッパー
を開き、湖岸線ラスタデータを入力ラスタとし
任意の出力ラスタデータ名を指定してください。

クリッピングモードを「領域」とします。

ここででてくる「1」と「2」のxy座標は
出力ラスタデータの
「1」はUpperLeft
「2」はLowerRight
の座標値にあたります。

この領域指定はラスタの切り取りだけでなく
領域拡張にも利用することができます。
(拡張範囲はNODATAとなります)

今回領域を拡大したい範囲を
X方向(東)に100セル分
Y方向(北)に100セル分
とし
ラスタの解像度を90mとしますと
90×100セル = 9000m
分XYの正の方向へ領域を拡大することになります。

この値を先ほどの「情報」ツールで得られた
UpperLeftのY座標値とLowerRightのX座標値に
加算して領域を指定すると以下のようになります。

「1」(Upperleft)
x:491925.257
y:4790017.207(北に+9000m)

「2」(LowerRight)
x:543495.257(東に+9000m),
y:4723777.207

上記の値を入力し処理を実行すると
湖岸線ピクセルの位置がずれることなく
領域のみが拡張された新たな湖岸線ラスタデータが得られるかと思います。
拡張領域のセル値はNODATA
(どの値がNODATAと定義されるかはラスタのデータタイプ次第です)
として定義されていることとおもいます。

もじばっかりな記事でアレですが、お試しください。

Sentinel-2のサンプル画像

久しぶりにESA(欧州宇宙機関)のSentinels Scientific Data Hubにアクセスしてみたら、Sentinel-2 Pre-operational Hub : pre-operational access point for all users to Sentinel-2 data.
というものができていました。いまのところユーザ登録がいらずにログインできるようになっているようですね。

Sentinel-2ってなんじゃらほいというひとはこちらをざっとみてみてください。Overviewをみてるだけで結構ぐっときます(私は)。https://en.wikipedia.org/wiki/Sentinel-2

Multi-spectral data with 13 bands もよいですが  free and opendata policy というのが素敵です。

アクセスしてみると日本を撮影したデータも公開されていたので、やんごとない人にデータを見てもらうためにサンプルサイトをさきほど作成しました。(忙しい人は、1clickでぱっと表示されるように用意しておかないと見てくれないし、見てくれないとその後の話も通じないのです)

福井県の三方五湖付近

2015/12/8に撮影されたデータの band 2,3,4 (空間解像度10m) を合成して作成してあります。表示可能ズームレベルは 5-15 です。

Sentinel-2についての日本語の解説(緒元概要)はこちら
https://www.restec.or.jp/satellite/sentinel-2-a-2-b

Copernicus Sentinel Data の Terms and Conditionsはこちら
https://sentinel.esa.int/documents/247904/690755/Sentinel_Data_Terms_and_Conditions

Q:簡単に流域界データをつくりたい

A: MapWindow をつかったらいいのでは(Windowsユーザなら)

MapWindow   http://www.mapwindow.org/

FOSS4G Advent Calendar 2015 二個目の記事として書いています。

「DEMから流域界ポリゴンを作成したいのだけど、QGISのプロセッシングツールボックス以外にもっと簡単な方法ない?」

ときかれたので、別件で使っていたMapWindowでためしにやってみたところさっくりできました。ツールを使うための設定等が一切いらないので、悩みたくない人にはおすすめです。(Windowsしばりなことが残念な点です)

DEMと流域界を定義したい場所のESRI Shapefile 形式のポイントデータがあればOKです。ここでDEMおよびポイントデータの座標系はUTM等のMetricのものである必要があることに注意してください。

DEMとポイントデータは、空のディレクトリを作成した後、その中に格納するとよいです。その場所が作業ディレクトリとなり以降出力される中間ファイルが格納されます。

mapwindow03

インストールのためのバイナリはここから入手してください。

http://mapwindow4.codeplex.com/releases/view/110244

インストール自体はさっくりおわります。起動した後、watershed delineation プラグインをアクティブにします。

mapwindow01

メニューバーにWatershed Delineation がでてきます。

処理プロセスを選択しながら実行するやり方もありますが、ここではいちばんかんたんなAutomaticを選択します。

mapwindow02

データ入力のためのウィンドウがでてくるので、

“Setup and Preprocessing” でDEMを選択します。ここではDEMの指定のみでほかの場所は触りません。

“Network Delineation by threshold method”では、処理プロセス中に地形(DEM)から流路(谷)をポリラインデータとして出力する際に、どのくらいの集水面積を持つ場所を流路と定義するかを決めます。ここではDEMのセル数(DEMの解像度は90mです)で定義しています(5000)。

“Custom Outlet/Inlet Layer”でポイントデータを指定します。(snap距離はデフォルトのままでだいたいOKと思われます)

mapwindow04
ファイルの指定が終わると下図のように表示されます。DEMはSRTM(Shuttle Radar Topography Mission)の北海道日高山脈付近のデータを利用しています。

mapwindow05

あとは「Run All」 をクリックするだけです。

しばらく待つと、さっくり流域界が出力されました。

mapwindow06
なおより細かく流域界を抽出したい場合は”Network Delineation by threshold method”の値を調整するとよいでしょう。以下の出力結果はセル数を500にしたものです。

mapwindow07

なお中間ファイルと出力結果のファイルは入力DEMと同じディレクトリに格納されています。

mapwindow08

おためしください。ではでは。

年内に終わらせたい宿題 (QGIS編)

選挙も終わって長らく気絶していた官公庁がもぞもぞと動き出すころですね。
各政党の公約は今後どのように実行されていくのでしょうか。

公約違反だ!…というセリフを聞くたびに
そういえば私も頼まれものをまだやっていないよなー、と思ってしまいます。

そこで、自分のところで滞留している頼まれ物リストのうち
FOSS4G(というかQGIS)関連で3つを選んで挙げて、自分のダメっぷりを晒していこうとおもいます。

++++++++++++++++++++++++++++++++++++++++

1) XXXXをQGISで見られるようにYYYYで配信する

のっけから伏せ字で申し訳ありません。

これ、ちまちまと作り出していたのですが
先週、あるかたがすでに同じテーマに取り組んでいることが判明し
しかも目指すクオリティがはるかに高い…
…これは素直にお任せすることにしました。

近々素敵なものが公開されるはずです。
楽しみにしてます!(他力本願)

2) QGIS上でポイントフィーチャを最近隣のポリゴン/ポリラインのエッジにフィットさせる(くっつける)

これをQGISでやりたいんだけど…という問い合わせを複数のかたから頂きました。
河川関係の仕事をやっているからでしょうかね。

たぶんこれとこれを組み合わせればできそうです。

NNjoin Plugin
最近隣エッジ上のlocationとそこまでの距離を返す

後者のほうを実行すると以下のような答えがかえってきます。


この結果に、フィットさせた先の
ポリゴンまたはポリラインのフィーチャIDを付与できるような処理が行えればOKです。

3) QGIS上でNETWORK解析の到達圏解析をする

魚類の仕事をしているかたから
QGIS上で、一定の条件を与えて(ex.距離,標高 etc..)
遡河回遊を行う魚類の遡上可能範囲(河川内)を抽出したい
とのご相談を受けました。

GRASSpgRoutingつかって到達圏解析をしたらどうですか?と依頼主に提案したら

“GRASSとか面倒くさくてイヤ!”

と若干かぶせぎみにいわれたので
GRASSのヘビーユーザとしてすこし(かなり)傷つきました。

1)といい2)といい私の周りの人々のQGIS好きは異常です

まあ最短経路解析(Shortest Path)ができるんだから
到達圏解析もできるんだろうな、と思い調べてみたら、
Network Analysis:Areas of the availability

…できそうですね。

サンプルコードを試しに実行してみた結果は以下のとおりです。
(なお下図の河川は適当にえらんだ場所であり、ここで解析を行うというわけではありません)
fig1
Input と Output をQGISのレイヤに変更してやればよさそうで
たとえば input のあたりは以下のように変えて
複数の地点を一気に処理できればよいかな
と思っています。

++++++++++++++++++++++++++++++++++++++++

結局 “QGISで” というしばりのために放っておいたものばかりですね。

1)は来週には公開されるとして

2) 3) はクリスマスまでに使えるようにしたいです。
….が、もしかしたらお年玉の時期になるかもしれません。

ぼちぼちとPyQGISをお勉強しながら取り組みたいと思います。

なお、解決してくれる素敵なツールやプラグインがあればぜひご一報ください。
(ものぐさなので大変助かります)

ではでは皆様よいお年を。

QGIS:OpenLayers Pluginの背景画像をコンポーザで印刷するとずれる

タイトルの質問がとどいたのでちょっとしらべてみました。
2年前から問題が指摘されていたみたい。
http://hub.qgis.org/issues/5827
#13のコメントの内容を再現できたのでメモしておきます。

コンポーザでたとえばA4サイズで以下のような地図を印刷するとします。
以下の画像のように余白を設定することはわりとよくありますね。
p001
そうすると重ねてあるベクタレイヤがずれて印刷されてしまいます。
質問の状況が再現できました。
p002
これはアイテムプロパティの”位置とサイズ”が以下のように
印刷サイズであるA4のサイズと食い違っているために発生するようです。
p003
印刷地図の余白がなくなってしまいますがとりあえず
“位置とサイズ”をA4サイズと同じ値に修正し印刷してみます。
p004
ずれは修正されているようです。
p003
OSMの地図を背景にした印刷需要は高そうなので
ちょっと不便だなーと思っているヒトもかなりいるだろうな。

航空レーザ測量データをひたすら ESRI shapefile に変換するだけ

ディレクトリに格納されている複数の航空レーザ測量のデータ(xxxx.txt)を
さっくりESRI shapefile形式(ポイント)に変換したいという相談があったので
pythonで書いてみました。

カンマで区切られたテキストの2列目がx 3列目がy 2列目がz
となっていたのでそれにあわせて書いてあります。

ogr2ogr使えばいいんだけど
筋トレ的に書いたので気にしない気にしない。

pyshpというライブラリが必要です。
コマンドプロンプトから
pip install pyshp
または
easy_install pyshp
でインストールできます。

pip も easy_install もインストールされていない場合は
このあたりからバイナリ(pyshp-1.2.0.win-amd64-py2.7.exe)をダウンロードして
インストールしてください。

その後以下をテキストエディタ等で適宜名前(ex. lp2shp.py)をつけて
Python27\Lib
直下に保存してください。

# 10/9追記:余計な処理をしてメモリを無駄に消費していたので修正



その後、python(command line) 等で



とテキストデータの入っているディレクトリを指定して実行してください。
同じディレクトリ内に xxxx_p.shp というファイル名でshapefileが出力されます。