ある点からの半径x 角度θの直線の書き方

[sakuya_izayoi]

こちらsakuya_izayoiのsakuyaです。お世話になります。
とある点からの半径x 角度θの直線の書き方をしりたいです。今まではcadで書いていたのですが、数が多かったのでフィールド計算機を活用して描けないかと思いましては。最終的には終点の緯度経度も求めたいのでラインの形で結果がほしいです。
すみませんがお願いします。

[sakuya_izayoi]

こちらizayoiです。

sakuyaへ。外からなので実際のフィールド計算内容を確かめられないので、申し訳ありませんが、考え方のみを記載します。

とりあえず、条件をまとめさせて貰うと終点の緯度経度が欲しいという事でいいのかな？

座標系が緯度経度だとちょっと面倒なので、

メルカトルか直法座標系に変換してもらった後に、

①真北をゼロ°とした場合で時計回りに角度が増えていくとする。その角度をθとする。

②線の開始座標が分かっている。開始座標を(Xs,Ys)とする。

③その線の長さは分かっている。その長さをRメートルとする。

という条件ならば、

開始座標を(Xs,Ys)とした場合に終点の座標を(Xe,Ye)とすると

終点のX座標、Xeは(R×sinθ)+Xs

終点のY座標、Yeは (R×cosθ)+Ys

となるはず。間違っていたら教えてください。

QGISでは角度はラジアンだから、θをラジアンに変換して計算させる。

角度とラジアンの変換式としては

ラジアン=θ×(π/180)=θ×(3.1415/180)=θ×0.01745

となる。

終点の座標を求めたら、ファイルから(shpファイルに保存したのちにdbfを使うのが楽かしら)XeとYeをカンマ区切りのファイルとして読み込ませれば終点座標は出来上がり。

ラインの書き方は分かりません。すみませんが始点と終点を力技で直接つないでください。

もっとスマートなやり方をご存知の方はすみませんがご投稿お願いいたします。

[adachi]

sakuya様、izayoi様

スマートな方法、なかなかなさそうですね…

pythonを使いこなせる方ならうまいコードを組めそうですが

基本的にはizayoi様ご提示のやり方しかなさそうな気がします。

ささやかな追加情報として

・角度→ラジアンの変換はradians() 関数でも可能

・始点群と終点群の両方のポイントレイヤに共通のIDを残しておけば、Connect pointsプラグインで自動的にラインは引ける

　（ただしこのプラグイン、できるラインレイヤのCRSが自動的にWGS84緯度経度（4326）になってしまうので、レイヤプロパティ画面でCRSを修正する必要あり）

ぐらいが省力化できそうなところでしょうか。

[sakuya_izayoi]

adashiさま

こちらsakuyaです。お世話になります。

izayoiのやり方も、最初は？が浮いていたのですが（sinとcosが逆じゃないかなと思っていたけど、先ほど解説受けて分かりました）どうにか理解できました。こういう時に数学って大事ですね。

・ラジアン変換の関数、ありがとうございました。izayoiの書き方だと四捨五入の時点で結構誤差が出そうなのでradians関数を使用させていただきます。

・プラグインの方なのですが、これは私はもちろん、izayoiも知らないものですね。ありがとうございます。

　そして、CRSが自動的にWGS84になってしまうというのがありますが、これはレイヤがWGS84で作られるという事？それとも、CRSが間違ってWGS84で定義されてしまうという事？

あと、izayoiの作成したラインの長さの表示についての投稿を受けて

https://groups.google.com/forum/#!topic/qgisshitumon01/TRRRgMKqK3o

なのですが、開始座標を(Xs,Ys)とした状態で、その点から500mの点を求めたとした場合、この500mというのは果たして何が500mになるのだろうという疑問が湧いてきました。

本日と明日外出なので帰宅してから検証してみたいと思います。izayoiで不明だった点があたしでできるかという問題も・・・。

[adachi]

sakuya様

正直私も最初見たとき一瞬、「逆じゃないか？」と思ったのですが

UTM、平面直角では一般にイメージするX軸、Y軸と逆になる、と言うことを改めて認識できました。

ご質問のプラグインの生成レイヤですが、

データ自体はインプットデータと同じCRSでできる（平面直角第XX系のデータを使うと、同じ座標系のデータが出力される）のですが

そのデータの表示上のCRSが強制的にデフォルトで4326になってしまう（=間違ったCRSが適用された状態）みたいです。

実際のデータのCRSと、適用されているCRSが違ってしまっている状態なので

レイヤプロパティ内でCRSを正しいものに変えると、ちゃんとしたデータになるかと思います。

[yoh_chan]

sakuyaさま：

こんにちは。

ご希望のラインの生成ですが、フィールド計算機を活用して、とのことですのでやってみました。

まず、フィールド計算機にはprojectという、あるポイントから距離と角度を指定して遷移する関数があります。

また、ラインの生成にはmake_lineという関数が利用できます。

このふたつを利用すれば、線を描くことができます。

ただし、IDなどのユニークな値を持つ共通したフィールドが必要になります。

以下、QGIS3.2でやってみたのでご参考いただければ。（平面直角12系でやったものです。）

QGISキャンバスに始点となるポイントレイヤ pt_2454 があるとします。（shpでもなんでも構いません。）

これをQGISレイヤメニューから[エクスポート]->[Save Features As]（名前をつけて保存）から、gpkgファイルに新たなレイヤとして保存します。（レイヤ名を ln_2454_gpkg とします。）

保存の際に、「ジオメトリタイプ」を「LineString」としてください。ただし、コピー元のレイヤが同じgpkgに入っているレイヤの場合は複製に失敗します。

（shpだと、ポイントのジオメトリを持つデータをラインストリングタイプとして保存することができないのですが、gpkgであれば可能です。）

すると、コピー元のポイントレイヤと同じ属性値を持ったラインストリングのレイヤができます。

ラインストリングのはずのレイヤのジオメトリカラムにポイントのジオメトリが入っているので、QGISでは地物の描画はされませんが、属性テーブルは参照できると思います。

同じgpkgに、こんどは「ジオメトリタイプ」を「Point」としたレイヤもエクスポートします。（ pt_2454_gpkg とします。）

この操作で、ふたつのレイヤにはデフォルトで「fid」というユニークかつ共通のIDフィールドが付加されていると思います。

（もともと共通のIDフィールドがある場合は、ポイントのインポート操作は不要です。）

その上で、レイヤ ln_2454_gpkg の属性テーブルからフィールド計算機を開きます。

「既存のフィールドを更新する」として、更新するカラムは<geometry>を選択します。

入力する数式は以下です。

内容としては、

　・地物のfidから対応する点のジオメトリをレイヤ pt_2454_gpkg から参照し、

　・これをproject関数で指定した長さ、角度で遷移した点を生成し、

　・両者をmake_lineで線で結ぶ

というものです。

make_line(  -- 始点と終点から線を生成する

geometry(  -- 始点のジオメトリ

get_feature(  -- 属性から地物を検索する

'pt_2454_gpkg'  -- QGIS上のレイヤ名

, 'fid'  -- 共通ID名

, attribute( $currentfeature, 'fid' )

)

)

, project(  -- 終点のジオメトリ

geometry(

get_feature(

'pt_2454_gpkg'

, 'fid'

, attribute( $currentfeature, 'fid' )

)

)

, 100  -- 距離

, radians(45)  -- ラジアン角

)

)

試しに作成したポイントのshpと、そこから作成したgpkgを貼付しておきます。

ご参考になればさいわいです。

[ありた]

sakuya さん、こんにちは

yoh_chan さんと似たような感じですが、フィールド計算機の使用だと、

(1) ラインデータとしてレイヤを作成。参照系は平面直角座標系

(2) たとえば距離 dist 、角度 azimuth という名前のフィールドを用意

(3) 基準としたい点を始点に、適当な点を終点とした地物を作成。 dist / azimuth に数値を入力

(4) フィールド計算機で <geometry> の更新を指定し、

make_line(start_point($geometry),project(start_point($geometry),"dist",radians("azimuth")))

(5) ツールバーの再読込アイコンをクリックするなり、マップ表示位置を変更するなりする

（再描画されていないだけ）

(6) フィールド計算機で lat / lng といった（仮想）フィールドを計算。 CRS は適宜変更

x(transform(end_point($geometry),'EPSG:2447','EPSG:4612'))

y(transform(end_point($geometry),'EPSG:2447','EPSG:4612'))

これですと、一応、ひとつのレイヤで完結します。

あるいは基準点の座標値、距離、方位角がすでに与えられており、機械的に処理したい場合は

エクセルなどで WKT の項をもった CSV を作るのが楽だと思います。

たとえば1行目をヘッダ行として geom / x0 / y0 / dist / azimuth を入力

A 列 (geom) を下記のような WKT の形式となるようにします。

="LINESTRING("&TEXT(B2,"0.000 ")&TEXT(C2,"0.000 , ")&TEXT(B2+D2*SIN(RADIANS(E2)),"0.000 ")&TEXT(C2+D2*COS(RADIANS(E2)),"0.000")&")"

CSV で保存したのち QGIS で読み込みます。

このときジオメトリは WKT 形式で geom を選択。

とかでしょうか。

[sakuya_izayoi]

sakuyaです。お世話になります。

adachiさま

やっぱりsinとcosが逆だと思いますよね。そうなんですよ。どこをθと考えるかで変わってくるというのを今日初めて認識しました。

入試レベルまでの数学だと、cosの方がx座標で、sinの方がy座標になるという問題ばかり見ていたから新鮮だったです。

プラグインの解説ありがとうございます。気を付けて実施してみます。

yoh_chanさま 　ありた さま

良い方法をありがとうございます。

やはり距離という数字を直接入れる場合は平面直角座標系に変換してやらないとちょっとまずい感じでしょうか？

そうなるとizayoiがちょろっと述べているメルカトルも距離が長い場合は下記のリンクのように誤差が大きくなってしまう？

https://groups.google.com/forum/#!topic/qgisshitumon01/TRRRgMKqK3o

yoh_chanさま、projectという関数は初めて知りました。そしてgeometry関数の中にさらに入れ子上に計算をさせるというのはやったことなかったです。これは式を自分では作れませんね。見本感謝いたします。ここの距離や角度もフィールドに入力してあればその通り計算させれそうですね。

ありた さま

最初に書いてある式の方は、QGIS上で直接ポイントを落としたのちに、長さと角度を設定する場合に使うやり方ですか？

なるほど、アプローチの仕方でyoh_chanさまご提示の式ともまた違った感じになるのですね。御二方の式を見ていると、make_line関数とproject関数がキモということですね。ここに思い至ることができませんでした。

後半のWKTはごくまれにizayoiが作っているのを見たことがあるのですがあたしは経験なかったです。そうですね。もとの座標が分かっていた場合はフィールド計算機にこだわらなくても良かったのですね。視野狭窄していました。ありがとうございます。そしてizayoi書いてくれていた式と同じような内容があるので、あってたって安心しそうです。

今回は表計算レベルで材料がそろえれる状況であったため、WTKファイルの方を使わさせていただきました。

ありがとうございます。