他にもやり方はありますが QGISですと processing にある「Create points along lines(ラインに沿ってポイントを作る)」というgdalのツールを使うと楽です。フィールド計算機でもポリラインと始点と終点の座標値は出力できますが、さくっとポイントにしたいよ、という需要もあるでしょう。たとえば、そのポイントデータを使って point sampling tool で各種レイヤの値を抽出する、などが想定されます。
さてこのツールですが、各ラインフィーチャの特定の地点にポイントを発生させるツールです。その地点はフィーチャの長さに対する割合で指定することができます。始点が0で終点が1として指定するので中間点は0.5ってことになりますね。なので始点と終点にポイントを発生させたい場合は以下の通りになります。なおポリラインフィーチャの属性はすべて出力されるポイントデータに受け継がれます。
仕事/実習/講義でネットワーク分析を取り上げる際によく使わせてもらっているGRASS GIS の v.net.iso コマンドですが、これの移動コストの設定方法について質問が来たので備忘録として書いておきます。
しかもQGISしばりでとのことだったのですが、3.4xでは processing の処理がコケるので2.18xでやることにします。
そんなのGRASSでやればいいんですよ！
v.net.iso のシンプルな使い方はこちらを見てみてください(資料が雑でごめんなさい)
v.net.iso の処理はデフォルトではネットワークデータであるポリラインの長さがコストとして取り扱われます。
たとえば逃げ地図のように「その距離を何分で歩けるか」ということを地図化したい場合など、距離ではない単位でコストを表現したい場合がありますよね(とはいえこの場合も距離は重要なコスト因子ですが)。
たとえばこんな単純なケースを考えてみます…とりあえずv.net.isoのコストとは、がわかればいいので地形要因とかはいったん置いときましょう。
課題：自宅から0.5時間(30分)/1時間(60分)/1.5時間(90分)で到達できる範囲をそれぞれ出力したい
設定：徒歩(時速4km) と ジョギング(時速8km)の2パターンで算出したい
これはポリラインデータに距離がm単位で格納されているlength列があるとして walkcost と jogcost 列を作ったうえで
walk_cost = length/4000
jog_cost = length/8000
を算出してこれをコストとして用いればよいです。もちろんジオメトリから直接算出してもよいです。
これはフィールド計算機で処理できます。
つまり下図のように各ポリラインフィーチャに「通過所要時間(単位:時間)」が属性値として格納されることになります。
そのうえで v.net.iso処理 を行う際には
processing で v.net.iso を選択し
ネットワークデータと自宅のポイントデータを入力したうえで
先ほどの walk_cost または jog_cost 列を
Arc forward/both direction(s) cost column (number)
に指定してやり
cost for isoline を
0.5,1.0,1.5
と求めたい時間で入力してやればよいです。
こんなかんじ。
すると徒歩の到達圏とジョギングの到達圏が求められますね。
もちろん各フィーチャ毎に制限速度や勾配などの移動コストに影響する属性があれば
それらをコストに変換する計算式を作成して新しいコスト列を作ってやればよいわけですね。
….またまだるっこしい説明をしてしまいましたが、まあよいです。
ではでは。
QGISのメニューバーの「ベクタ」→「ジオメトリツール」→「ボロノイポリゴン」で、ポイントデータからボロノイ図が出力できるのですが、その入力パラメータとして「バッファ領域」というものがあります。
「これなんなん？」と問われたので実データでチェックしました。下図を参照ください。メニューバーの「ベクタ」→「調査ツール」→「レイヤ領域のポリゴン」(Layer extentを矩形で出力する)の出力結果が赤枠、バッファ領域を10%と設定して出力したボロノイ図が緑色のポリゴンになります。
(QGIS 2.14 x Windows 8.1)
先日宮城と岩手で行ったQGISハンズオンで、leaflet.jsベースのウェブマップを簡単に出力できる qgis2web プラグインを紹介しました。その際、背景図をOpenStreetMap以外に変更できないか？という質問がありました。帰りがけだったため、ざざっと「ここのURLを書き換えるだけでOKです」と説明したのみでしたので、再度その内容をここに書いておきます。
index.htmlをテキストエディタで開いて35-40行目付近(地図に載せたレイヤの数により前後します)にあるこの2行を書き換えれば地理院地図の背景図に変換できます。
<変更前：OSM>
|
1 2 |
var basemap0 = L.tileLayer('http://{s}.tile.openstreetmap.org/{z}/{x}/{y}.png', { attribution: '© <a href="http://openstreetmap.org">OpenStreetMap</a> contributors,<a href="http://creativecommons.org/licenses/by-sa/2.0/">CC-BY-SA</a>', |
<変更後：地理院タイル(標準)>
|
1 2 |
var basemap0 = L.tileLayer('http://cyberjapandata.gsi.go.jp/xyz/std/{z}/{x}/{y}.png', { attribution: '<a href="http://www.gsi.go.jp/kikakuchousei/kikakuchousei40182.html" target="_blank" rel="noopener noreferrer">国土地理院</a>', |
なお地理院地図で提供されているタイル地図であれば上記の
‘http://cyberjapandata.gsi.go.jp/xyz/std/{z}/{x}/{y}.png’
に相当する箇所のURLを変えるだけで背景を変更できます。URLのリストは以下のリンク先を参照ください。
これは 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でこれらのコマンドと同じ処理を行いたいなーというときには
という手続きを踏んでいました(私は)。
もちろんこれでも十分楽に扱えるのですが、
大変ものぐさなわたしは
もすこしシンプルに同等の処理を書けたらいいなー
とか思っておりました。
これが誰かに伝わったのか(そんなわけないですが
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.
とあります。
まえよりわかりやすくオプションを指定できる便利なラッパーができたよー とのことですね。
それではみてみましょう。
いちばん単純なフォーマットの変換処理を例に取りますと、
|
1 |
gdal_translate -of HFA C:/temp/in.tif C:/temp/out.img |
というコマンドと同じ処理を、いままでPythonでこんな風に書いていました(私は)。
|
1 2 3 4 |
from osgeo import gdal srcrs = gdal.Open("C:/temp/in.tif") driver = gdal.GetDriverByName("HFA") dstrss = driver.CreateCopy("C:/temp/out.img", srcrs) |
んで、GDAL 2.1 x Python ではどーなるかといいますと..
|
1 2 |
from osgeo import gdal gdal.Translate("C:/temp/out.img", "C:/temp/in.tif", format = 'HFA') |
と、とてもシンプルに書ける様になりました。
もひとつ例を挙げるとラスタデータの投影変換はこんなコマンドですが
|
1 |
gdalwarp -of GTiff -s_srs EPSG:4326 -t_srs EPSG:32654 -tr 90 90 -r bilinear C:/temp/in_gcs.tif C:/temp/out_u54.tif |
GDAL 2.1 x Pythonでは…
|
1 2 3 4 |
from osgeo import gdal gdal.Warp("C:/temp/out_u54.tif","C:/temp/in_gcs.tif", srcSRS="EPSG:4326",dstSRS="EPSG:32654",xRes=90,yRes=90, resampleAlg="bilinear") |
という感じです。オプションの指定がとてもわかりやすくてヨイです。
もちろんラスタデータに対して複雑な処理をするときはarray(配列)にしたほうが便利なのですが、単純なデータ処理(フォーマットの変換 投影変換 クリップなどなど)をしたいときにすっきり書けてうれしいです。
その他の処理やオプションについての詳細はこちらを参照ください。
http://gdal.org/python/
すこし古い話ではありますが私にとってありがたいアップデートだったので書いてみました。
ではでは。よいお年を（まだ早い
とあるところでタイトルのとおりの処理のやり方を提案してみたところ
ぽつぽつとリアクション(そーやってもできるのか！的な)があったりしたので書いてみます。
コネタでも書いておくのだいじなのね..と思いました。
さて、とある湖岸線(琵琶湖,支笏湖などなど)を示すラスタデータがあるとします。
湖岸線を示すピクセル値が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と定義されるかはラスタのデータタイプ次第です)
として定義されていることとおもいます。
もじばっかりな記事でアレですが、お試しください。
久しぶりに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
A: MapWindow をつかったらいいのでは（Windowsユーザなら)
MapWindow http://www.mapwindow.org/
FOSS4G Advent Calendar 2015 二個目の記事として書いています。
「DEMから流域界ポリゴンを作成したいのだけど、QGISのプロセッシングツールボックス以外にもっと簡単な方法ない？」
ときかれたので、別件で使っていたMapWindowでためしにやってみたところさっくりできました。ツールを使うための設定等が一切いらないので、悩みたくない人にはおすすめです。(Windowsしばりなことが残念な点です)
DEMと流域界を定義したい場所のESRI Shapefile 形式のポイントデータがあればOKです。ここでDEMおよびポイントデータの座標系はUTM等のMetricのものである必要があることに注意してください。
DEMとポイントデータは、空のディレクトリを作成した後、その中に格納するとよいです。その場所が作業ディレクトリとなり以降出力される中間ファイルが格納されます。
インストールのためのバイナリはここから入手してください。
http://mapwindow4.codeplex.com/releases/view/110244
インストール自体はさっくりおわります。起動した後、watershed delineation プラグインをアクティブにします。
メニューバーにWatershed Delineation がでてきます。
処理プロセスを選択しながら実行するやり方もありますが、ここではいちばんかんたんなAutomaticを選択します。
データ入力のためのウィンドウがでてくるので、
“Setup and Preprocessing” でDEMを選択します。ここではDEMの指定のみでほかの場所は触りません。
“Network Delineation by threshold method”では、処理プロセス中に地形(DEM)から流路(谷)をポリラインデータとして出力する際に、どのくらいの集水面積を持つ場所を流路と定義するかを決めます。ここではDEMのセル数(DEMの解像度は90mです)で定義しています(5000)。
“Custom Outlet/Inlet Layer”でポイントデータを指定します。(snap距離はデフォルトのままでだいたいOKと思われます)
ファイルの指定が終わると下図のように表示されます。DEMはSRTM(Shuttle Radar Topography Mission)の北海道日高山脈付近のデータを利用しています。
あとは「Run All」 をクリックするだけです。
しばらく待つと、さっくり流域界が出力されました。
なおより細かく流域界を抽出したい場合は”Network Delineation by threshold method”の値を調整するとよいでしょう。以下の出力結果はセル数を500にしたものです。
なお中間ファイルと出力結果のファイルは入力DEMと同じディレクトリに格納されています。
おためしください。ではでは。
選挙も終わって長らく気絶していた官公庁がもぞもぞと動き出すころですね。
各政党の公約は今後どのように実行されていくのでしょうか。
公約違反だ！…というセリフを聞くたびに
そういえば私も頼まれものをまだやっていないよなー、と思ってしまいます。
そこで、自分のところで滞留している頼まれ物リストのうち
FOSS4G(というかQGIS)関連で3つを選んで挙げて、自分のダメっぷりを晒していこうとおもいます。
++++++++++++++++++++++++++++++++++++++++
1) XXXXをQGISで見られるようにYYYYで配信する
のっけから伏せ字で申し訳ありません。
これ、ちまちまと作り出していたのですが
先週、あるかたがすでに同じテーマに取り組んでいることが判明し
しかも目指すクオリティがはるかに高い…
…これは素直にお任せすることにしました。
近々素敵なものが公開されるはずです。
楽しみにしてます！(他力本願)
2) QGIS上でポイントフィーチャを最近隣のポリゴン/ポリラインのエッジにフィットさせる(くっつける)
これをQGISでやりたいんだけど…という問い合わせを複数のかたから頂きました。
河川関係の仕事をやっているからでしょうかね。
たぶんこれとこれを組み合わせればできそうです。
NNjoin Plugin
最近隣エッジ上のlocationとそこまでの距離を返す
後者のほうを実行すると以下のような答えがかえってきます。
|
1 |
Closest point found at: POINT (0.5 1), with a distance of 0.50 units. |
この結果に、フィットさせた先の
ポリゴンまたはポリラインのフィーチャIDを付与できるような処理が行えればOKです。
3) QGIS上でNETWORK解析の到達圏解析をする
魚類の仕事をしているかたから
QGIS上で、一定の条件を与えて(ex.距離,標高 etc..)
遡河回遊を行う魚類の遡上可能範囲(河川内)を抽出したい
とのご相談を受けました。
GRASSかpgRoutingつかって到達圏解析をしたらどうですか？と依頼主に提案したら
“GRASSとか面倒くさくてイヤ！”
と若干かぶせぎみにいわれたので
GRASSのヘビーユーザとしてすこし（かなり）傷つきました。
1)といい2)といい私の周りの人々のQGIS好きは異常です
まあ最短経路解析(Shortest Path)ができるんだから
到達圏解析もできるんだろうな、と思い調べてみたら、
Network Analysis:Areas of the availability
…できそうですね。
サンプルコードを試しに実行してみた結果は以下のとおりです。
(なお下図の河川は適当にえらんだ場所であり、ここで解析を行うというわけではありません)
Input と Output をQGISのレイヤに変更してやればよさそうで
たとえば input のあたりは以下のように変えて
複数の地点を一気に処理できればよいかな
と思っています。
|
1 2 3 4 5 6 7 |
vlayer = QgsVectorLayer("river_mouth.shp", "river_mouth", "ogr") iter = vlayer.getFeatures() for feature in iter: geom = feature.geometry() x1 = geom.asPoint().x() y1 = geom.asPoint().y() pStart = QgsPoint( x1, y1 ) |
++++++++++++++++++++++++++++++++++++++++
結局 “QGISで” というしばりのために放っておいたものばかりですね。
1)は来週には公開されるとして
2) 3) はクリスマスまでに使えるようにしたいです。
….が、もしかしたらお年玉の時期になるかもしれません。
ぼちぼちとPyQGISをお勉強しながら取り組みたいと思います。
なお、解決してくれる素敵なツールやプラグインがあればぜひご一報ください。
（ものぐさなので大変助かります）
ではでは皆様よいお年を。