postGISで接触する線を抽出する方法を探す

2023年5月5日

 

プレゼンテーション1

http://cse.naro.affrc.go.jp/yellow/pgisman/2.2.0/ST_ClosestPoint.html」

http://postgres.cn/docs/postgis-2.3/ST_Buffer.html

このジオメトリ/ジオグラフィからの距離が指定された距離以下となる点全てを表現するジオメトリ/ジオグラフィを返します。

■線の接触を検知するテスト用の線

ちなみに、現実世界の座標では、0,150や100,100などの座標はないのだけど、

create extension postgis
create extension pgrouting

を仕込んだpostGISのDBであれば、こんな訳の分からない座標でもちゃんと計算してくれます。

■線に幅10を与える方法

yoko_db=#
yoko_db=# SELECT ST_Buffer(
yoko_db(# ST_GeomFromText(
yoko_db(# 'LINESTRING(50 50,150 150,150 50)'
yoko_db(# ), 10, 'endcap=round join=round');
st_buffer
-----------------------------------------------------------------------------
01030000000100000031000000AA07FFCFB9DD614056F8003046A26340E82FBFAF370E624050929AFD11CA6340B196BE928A45624085F91736A4E76340E1F76833928162409FA91FEED9F963400000000000C0624000000000000064401F0897CC6DFE62409FA91FEED9F963404F69416D753A634085F91736A4E7634018D04050C871634050929AFD11CA634056F8003046A2634056F8003046A2634050929AFD11CA634018D04050C871634085F91736A4E763404F69416D753A63409FA91FEED9F963401F0897CC6DFE624000000000000064400000000000C06240000000000000644000000000000049409FA91FEED9F9634082DFA3CD4806484085F91736A4E76340C35AFA4A2A16474050929AFD11CA6340A0BFFCBEDE38464056F8003046A26340A71EFC3FE776454018D04050C8716340BEB69509B8D744404F69416D753A6340ED19A0276F6144401F0897CC6DFE624084598147981844400000000000C062400000000000004440E1F76833928162408459814798184440B196BE928A456240EC19A0276F614440E92FBFAF370E6240BEB69509B8D74440AA07FFCFB9DD6140A71EFC3FE7764540B06D6502EEB561409FBFFCBEDE3846407B06E8C95B986140C25AFA4A2A1647406156E0112686614081DFA3CD48064840000000000080614000000000000049400000000000806140A81EFC3FE7765F4058E103C018894C40A81EFC3FE77645405E40034121C74B40BEB69509B8D744403AA505B5D5E94A40EC19A0276F6144407B205C32B7F949408459814798184440FFFFFFFFFFFF4840000000000000444083DFA3CD480648408459814798184440C45AFA4A2A164740ED19A0276F614440A1BFFCBEDE384640BEB69509B8D74440A81EFC3FE7764540A81EFC3FE7764540BFB69509B8D74440A0BFFCBEDE384640EE19A0276F614440C45AFA4A2A164740855981479818444082DFA3CD480648400000000000004440FFFFFFFFFFFF484084598147981844407A205C32B7F94940EC19A0276F6144403AA505B5D5E94A40BDB69509B8D744405D40034121C74B40A81EFC3FE776454058E103C018894C40AA07FFCFB9DD614056F8003046A26340
(1 row)

■幅10の線を交差させる

SELECT ST_AsEWKT(ST_Intersection(ST_Buffer(ST_GeomFromText('LINESTRING(50 50,150 150,150 50)'), 10, 'endcap=round join=round'), ST_Buffer(ST_GeomFromText('LINESTRING(0 100,100 100,150 150,0 150)'), 10, 'endcap=round join=round')));

で、その出力値をグラフにするとこうなっていました。

重ね併わせてみます。

線の幅を1にしてみました。

SELECT ST_AsEWKT(ST_Intersection(ST_Buffer(ST_GeomFromText('LINESTRING(50 50,150 150,150 50)'), 1, 'endcap=round join=round'), ST_Buffer(ST_GeomFromText('LINESTRING(0 100,100 100,150 150,0 150)'), 1, 'endcap=round join=round')));

併わせてみます。

■領域と線の点(ブルーの線)

# SELECT ST_AsEWKT(ST_Intersection(ST_Intersection(ST_Buffer(ST_GeomFromText('LINESTRING(50 50,150 150,150 50)'), 10, 'endcap=round join=round'), ST_Buffer(ST_GeomFromText('LINESTRING(0 100,100 100,150 150,0 150)'), 10, 'endcap=round join=round')),ST_GeomFromText('LINESTRING(0 100,100 100,150 150, 0 150)')));

st_asewkt

----------------------------------------------------------------------

LINESTRING(85.857864376269 100,100 100,150 150,135.857864376269 150)

(1 row)

■領域と線の点(オレンジの線)

yoko_db=# SELECT ST_AsEWKT(ST_Intersection(ST_Intersection(ST_Buffer(ST_GeomFromText('LINESTRING(50 50,150 150,150 50)'), 10, 'endcap=round join=round'), ST_Buffer(ST_GeomFromText('LINESTRING(0 100,100 100,150 150,0 150)'), 10, 'endcap=round join=round')),ST_GeomFromText('LINESTRING(50 50,150 150, 150 50)')));

st_asewkt

------------------------------------------------

LINESTRING(90 90,150 150,150 135.857864376269)

(1 row)

■実世界におけるST_Bufferが作る"幅"は?

先程、『現実世界の座標では、0,150や100,100などの座標はないのだけど、create extension postgis
、create extension pgroutingを仕込んだpostGISのDBであれば、こんな訳の分からない座標でもちゃんと計算してます』と記載しましたが、現実世界の座標では、どうやってこのサイズを決めればいいのかが分かりません。

という訳で実際の座標を使って試してみます。

139.622931 35.489653 と、 139.6401889 35.502831をつかって試してみます。

yoko_db=# SELECT ST_Distance('SRID=4326;POINT(139.622931 35.489653)'::GEOGRAPHY,'SRID=4326;POINT(139.6401889 35.502831)'::GEOGRAPHY);
st_distance
---------------
2142.35290029
(1 row)

ざっくり2km強ですね。

で、ここで、「描画しないと確認が面倒だなぁ」ということに気がつきました。どうやらpgAdminを使えばこれが可能なようです(先程、最新版をインストールして、古いバージョンをアンインストール(デジレクトリレベルで手動で消去しないと、最新版の起動に失敗します)

上記は古いpgadmin4のビューです。新しいpgAdmin4では、これで地図が表示されるようです

 

 

これで、SQL文を入力すると、その絵が出てくる環境が整いました。

で、まずは、

select ST_Buffer(ST_GeomFromText('LINESTRING(139.622931 35.489653, 35.502831 139.6401889)'),10, 'endcap=round join=round')

をやってみました。

定規で図ってみたところ、1/10くらいの幅になっているので、実際この場合、この線の幅は200メートルくらいでしょう。

誤差として50メートルとして、"2"くらいで試してみます。

select ST_Buffer(ST_GeomFromText('LINESTRING(139.622931 35.489653, 35.502831 139.6401889)'),2, 'endcap=round join=round')

よく分からなけど、GPSの誤差を鑑みて、このくらいの幅を取っておくべきかな。とりあえず、これをデフォルトとして、まずはDBを作ってみます。

テストとして、バスルートを記載してみます。

select ST_GeomFromText('LINESTRING(139.595873 35.378154,139.595779 35.37812,139.595766 35.378045,139.595922 35.377752,139.596212 35.377193,139.596704 35.376296,139.596804 35.376139,139.597005 35.375862,139.597443 35.375322,139.597832 35.374847,139.598029 35.374589,139.598145 35.374415,139.598302 35.374138,139.598369 35.374024,139.598422 35.373919,139.598491 35.373738,139.599019 35.372353,139.599495 35.371156,139.599468 35.371081,139.599338 35.371056,139.59899 35.370963,139.597059 35.370473,139.595393 35.370028,139.594419 35.369764,139.59407 35.369659,139.593798 35.369572,139.593476 35.369457,139.592966 35.369259,139.592779 35.369176,139.592465 35.369025,139.591683 35.368627,139.590731 35.368134,139.590073 35.367794,139.589176 35.367331,139.588262 35.366863,139.587735 35.366595,139.587508 35.366492,139.587332 35.366423,139.587176 35.366366,139.586864 35.36627,139.586549 35.36619,139.586207 35.366124,139.585855 35.366079,139.58554 35.366056,139.585204 35.366045,139.584925 35.366056,139.584641 35.366077,139.584445 35.366098,139.584247 35.366129,139.583902 35.366197,139.583561 35.366286,139.583199 35.366401,139.582906 35.366519,139.582565 35.366683,139.582309 35.366816,139.582119 35.366939,139.581933 35.36708,139.581569 35.367314,139.580995 35.367738,139.58035 35.368227,139.57968 35.368715,139.579054 35.36918,139.57882 35.369336,139.578604 35.369493,139.578352 35.369657,139.57811 35.369804,139.577878 35.369934,139.577738 35.37001,139.577589 35.370082,139.577418 35.37016,139.57732 35.37024,139.577346 35.370322,139.577492 35.37064,139.577597 35.370922,139.577664 35.371296,139.577672 35.371576,139.577637 35.371872,139.577594 35.372092,139.577509 35.37234,139.577434 35.372512,139.577012 35.373505,139.576949 35.373583,139.576554 35.374506,139.576609 35.374598,139.576838 35.374662,139.577089 35.374742,139.577114 35.374811,139.577011 35.375071,139.576926 35.375085,139.576557 35.374996,139.576471 35.374952,139.576471 35.374876,139.576599 35.374564,139.576916 35.37381,139.577 35.373789,139.577313 35.373871,139.577775 35.373993,139.57802 35.374046,139.578323 35.374087,139.578673 35.37412,139.578925 35.374131,139.579275 35.374129,139.579628 35.374106,139.580017 35.374068,139.581646 35.373846,139.582008 35.373789,139.582074 35.373858,139.581963 35.374391,139.581837 35.374905,139.581765 35.375135,139.581633 35.375534,139.581499 35.375863,139.58136 35.376167,139.581301 35.376281,139.580964 35.376851,139.580326 35.377869,139.578402 35.380965,139.577251 35.382825,139.577101 35.383021,139.577006 35.383152,139.576828 35.383389,139.576656 35.383641,139.576268 35.38417,139.575843 35.384739,139.575346 35.385415,139.57515 35.385681,139.574998 35.385894,139.574891 35.386055,139.574838 35.386181,139.574775 35.386347,139.574692 35.386563,139.574638 35.386697,139.574542 35.386724,139.574012 35.386604,139.572348 35.386248,139.571268 35.386011,139.570946 35.385959,139.570646 35.385931,139.570348 35.385922,139.570073 35.385934,139.569805 35.385966,139.569487 35.38603,139.568931 35.386182,139.568571 35.386286,139.568349 35.386354,139.568135 35.386442,139.56792 35.386576,139.567769 35.386716,139.56766 35.386853,139.56657 35.388671,139.566268 35.389177,139.566173 35.389362,139.566128 35.389529,139.566101 35.389744,139.566109 35.389929,139.566138 35.390063,139.566177 35.390179,139.566254 35.390338,139.567047 35.391736,139.567153 35.391917,139.567294 35.39212,139.567409 35.392281,139.56767 35.392568,139.568151 35.393053,139.568317 35.393216,139.568846 35.393747,139.569264 35.394174,139.569335 35.394232,139.569419 35.394281,139.569508 35.394292,139.570315 35.394297,139.571173 35.3943,139.571345 35.394307,139.571454 35.394351,139.57165 35.394627,139.572146 35.395394,139.573316 35.397176,139.573368 35.39726,139.573414 35.397345,139.573465 35.397477,139.573485 35.397559,139.573544 35.398079,139.573553 35.398228,139.573516 35.398354,139.573465 35.398447,139.57337 35.398541,139.57296 35.398854,139.57273 35.399024,139.570967 35.400511,139.570926 35.400555,139.570918 35.400594,139.571046 35.400985,139.57111 35.401124,139.571253 35.401304,139.571308 35.401365,139.571382 35.401397,139.571453 35.401399,139.57153 35.401349,139.571716 35.401204,139.571767 35.401149,139.571763 35.401079,139.57172 35.401031,139.571503 35.400935,139.571434 35.40092,139.57138 35.400931,139.57135 35.400995)')

では、次に、このルートに幅を与えてみます。
select ST_Buffer(ST_GeomFromText('LINESTRING(139.595873 35.378154,139.595779 35.37812,139.595766 35.378045,139.595922 35.377752,139.596212 35.377193,139.596704 35.376296,139.596804 35.376139,139.597005 35.375862,139.597443 35.375322,139.597832 35.374847,139.598029 35.374589,139.598145 35.374415,139.598302 35.374138,139.598369 35.374024,139.598422 35.373919,139.598491 35.373738,139.599019 35.372353,139.599495 35.371156,139.599468 35.371081,139.599338 35.371056,139.59899 35.370963,139.597059 35.370473,139.595393 35.370028,139.594419 35.369764,139.59407 35.369659,139.593798 35.369572,139.593476 35.369457,139.592966 35.369259,139.592779 35.369176,139.592465 35.369025,139.591683 35.368627,139.590731 35.368134,139.590073 35.367794,139.589176 35.367331,139.588262 35.366863,139.587735 35.366595,139.587508 35.366492,139.587332 35.366423,139.587176 35.366366,139.586864 35.36627,139.586549 35.36619,139.586207 35.366124,139.585855 35.366079,139.58554 35.366056,139.585204 35.366045,139.584925 35.366056,139.584641 35.366077,139.584445 35.366098,139.584247 35.366129,139.583902 35.366197,139.583561 35.366286,139.583199 35.366401,139.582906 35.366519,139.582565 35.366683,139.582309 35.366816,139.582119 35.366939,139.581933 35.36708,139.581569 35.367314,139.580995 35.367738,139.58035 35.368227,139.57968 35.368715,139.579054 35.36918,139.57882 35.369336,139.578604 35.369493,139.578352 35.369657,139.57811 35.369804,139.577878 35.369934,139.577738 35.37001,139.577589 35.370082,139.577418 35.37016,139.57732 35.37024,139.577346 35.370322,139.577492 35.37064,139.577597 35.370922,139.577664 35.371296,139.577672 35.371576,139.577637 35.371872,139.577594 35.372092,139.577509 35.37234,139.577434 35.372512,139.577012 35.373505,139.576949 35.373583,139.576554 35.374506,139.576609 35.374598,139.576838 35.374662,139.577089 35.374742,139.577114 35.374811,139.577011 35.375071,139.576926 35.375085,139.576557 35.374996,139.576471 35.374952,139.576471 35.374876,139.576599 35.374564,139.576916 35.37381,139.577 35.373789,139.577313 35.373871,139.577775 35.373993,139.57802 35.374046,139.578323 35.374087,139.578673 35.37412,139.578925 35.374131,139.579275 35.374129,139.579628 35.374106,139.580017 35.374068,139.581646 35.373846,139.582008 35.373789,139.582074 35.373858,139.581963 35.374391,139.581837 35.374905,139.581765 35.375135,139.581633 35.375534,139.581499 35.375863,139.58136 35.376167,139.581301 35.376281,139.580964 35.376851,139.580326 35.377869,139.578402 35.380965,139.577251 35.382825,139.577101 35.383021,139.577006 35.383152,139.576828 35.383389,139.576656 35.383641,139.576268 35.38417,139.575843 35.384739,139.575346 35.385415,139.57515 35.385681,139.574998 35.385894,139.574891 35.386055,139.574838 35.386181,139.574775 35.386347,139.574692 35.386563,139.574638 35.386697,139.574542 35.386724,139.574012 35.386604,139.572348 35.386248,139.571268 35.386011,139.570946 35.385959,139.570646 35.385931,139.570348 35.385922,139.570073 35.385934,139.569805 35.385966,139.569487 35.38603,139.568931 35.386182,139.568571 35.386286,139.568349 35.386354,139.568135 35.386442,139.56792 35.386576,139.567769 35.386716,139.56766 35.386853,139.56657 35.388671,139.566268 35.389177,139.566173 35.389362,139.566128 35.389529,139.566101 35.389744,139.566109 35.389929,139.566138 35.390063,139.566177 35.390179,139.566254 35.390338,139.567047 35.391736,139.567153 35.391917,139.567294 35.39212,139.567409 35.392281,139.56767 35.392568,139.568151 35.393053,139.568317 35.393216,139.568846 35.393747,139.569264 35.394174,139.569335 35.394232,139.569419 35.394281,139.569508 35.394292,139.570315 35.394297,139.571173 35.3943,139.571345 35.394307,139.571454 35.394351,139.57165 35.394627,139.572146 35.395394,139.573316 35.397176,139.573368 35.39726,139.573414 35.397345,139.573465 35.397477,139.573485 35.397559,139.573544 35.398079,139.573553 35.398228,139.573516 35.398354,139.573465 35.398447,139.57337 35.398541,139.57296 35.398854,139.57273 35.399024,139.570967 35.400511,139.570926 35.400555,139.570918 35.400594,139.571046 35.400985,139.57111 35.401124,139.571253 35.401304,139.571308 35.401365,139.571382 35.401397,139.571453 35.401399,139.57153 35.401349,139.571716 35.401204,139.571767 35.401149,139.571763 35.401079,139.57172 35.401031,139.571503 35.400935,139.571434 35.40092,139.57138 35.400931,139.57135 35.400995)'),1, 'endcap=round join=round')

あれ? この円は何? と思い色々弄っていたのですが

を、変えてみました。

としてみました。

"当たり"だったようです。

ただ、これでは幅をメートルで設定できないので、メートル設定をしてみました.

GPSの誤差として50メートルくらいでいいかな、と

SELECT ST_Buffer(st_transform(st_setsrid(ST_GeomFromtext(
'linestring(139.595873 35.378154,139.595779 35.37812,139.595766 35.378045,..... ,139.57138 35.400931,139.57135 35.400995)'
),4326),3857),50,'endcap=round join=round')

 

2023年5月5日2023,江端さんの技術メモ

Posted by ebata