2022/05,江端さんの技術メモ

// go get github.com/lib/pq を忘れずに
// go run main6.go

/*
	前提は、https://ebata-books.booth.pm/items/3351408 の環境

	(1)スタート地点  139.91804656152837 36.57246810341353,
	(2)ゴール地点 139.93515104919825 36.55883778950532
	(3)とした時、この2点に最も近いノードを探し出して、
	(4)その2点をダイクストラ法で、最短距離探索をするプログラム

	分かったこと
	上記(3)のノードが、孤立ノードになっている場合があり、この場合(4)が実施できないので、
	ダイクストラに失敗する場合は、別のノードでトライアルする必要あり 
	(このコードでは展開しないけど)
*/

package main

import (
	"database/sql"
	"fmt"
	"log"

	_ "github.com/lib/pq"
)

var source int
var longitude float64
var latitude float64
var dist float64

func main() {
	db, err := sql.Open("postgres",
		"user=postgres password=password host=localhost port=15432 dbname=utsu_db sslmode=disable")
	if err != nil {
		log.Fatal("OpenError: ", err)
	}
	defer db.Close()

	// スタート地点  139.91804656152837 36.57246810341353
	rows, err := db.Query("SELECT source, x1 as longitude, y1 as latitude, ST_Distance('SRID=4326;POINT(139.91804656152837 36.57246810341353)'::GEOGRAPHY, the_geom) as dist FROM ways WHERE ST_DWithin(the_geom, ST_GeographyFromText('SRID=4326;POINT(139.91804656152837 36.57246810341353)'), 300.0) ORDER BY dist LIMIT 1")
	if err != nil {
		log.Fatal(err)
	}
	defer rows.Close()

	for rows.Next() {

		if err := rows.Scan(&source, &longitude, &latitude, &dist); err != nil {
			fmt.Println(err)
		}
		fmt.Println(source, longitude, latitude, dist)
	}

	o_source := source

	// ゴール地点 139.93515104919825 36.55883778950532
	rows, err = db.Query("SELECT source, x1 as longitude, y1 as latitude, ST_Distance('SRID=4326;POINT(139.93515104919825 36.55883778950532)'::GEOGRAPHY, the_geom) as dist FROM ways WHERE ST_DWithin(the_geom, ST_GeographyFromText('SRID=4326;POINT(139.93515104919825 36.55883778950532)'), 300.0) ORDER BY dist LIMIT 1")
	if err != nil {
		log.Fatal(err)
	}
	defer rows.Close()

	for rows.Next() {

		if err := rows.Scan(&source, &longitude, &latitude, &dist); err != nil {
			fmt.Println(err)
		}
		fmt.Println(source, longitude, latitude, dist)
	}

	d_source := source

	fmt.Println(o_source, d_source)

	// これが基本形
	// rows, err = db.Query("SELECT seq, node, edge, cost FROM pgr_dijkstra('SELECT gid as id, source, target,length as cost FROM ways',8889, 22848, false)")

	// この bigint なるものは、https://wp.kobore.net/%E6%B1%9F%E7%AB%AF%E3%81%95%E3%82%93%E3%81%AE%E6%8A%80%E8%A1%93%E3%83%A1%E3%83%A2/error-function-pgr_dijkstraunknown-start_vid-bigint-end_vid-bigint-directed-boolean-does-not-exist-line-7-pgr_dijkstra-hint-no-function-matches-the-given-name/ に記載がある
	rows, err = db.Query("SELECT seq, node, edge, cost FROM pgr_dijkstra('SELECT gid as id, source, target,length as cost FROM ways', $1::bigint , $2::bigint , false)", o_source, d_source)

	// こっちも稼動します
	// str := "SELECT seq, node, edge, cost FROM pgr_dijkstra('SELECT gid as id, source, target,length as cost FROM ways'," + strconv.Itoa(o_source) + "," + strconv.Itoa(d_source) + ", false)"
	// fmt.Println(str)
	// rows, err = db.Query(str)

	if err != nil {
		log.Fatal(err)
	}
	defer rows.Close()

	for rows.Next() {
		var seq int
		var node int
		var edge int
		var cost float64

		if err := rows.Scan(&seq, &node, &edge, &cost); err != nil {
			fmt.Println(err)
		}
		fmt.Println(seq, node, edge, cost)
	}

	if err := db.Ping(); err != nil {
		log.Fatal("PingError: ", err)
	}
}

2022/05,江端さんの技術メモ

// go get github.com/lib/pq を忘れずに
// go run main5.go
package main

import (
	"database/sql"
	"fmt"
	"log"

	_ "github.com/lib/pq"
)

func main() {
	db, err := sql.Open("postgres", "user=postgres password=password host=localhost port=15432 dbname=utsu_db sslmode=disable")
	if err != nil {
		log.Fatal("OpenError: ", err)
	}
	defer db.Close()

	rows, err := db.Query("SELECT gid,tag_id,x1 FROM ways")
	if err != nil {
		log.Fatal(err)
	}

	// 139.9182893339256, 36.573831584767085

	defer rows.Close()

	for rows.Next() {
		var gid int
		var tag_id int
		var x1 float64
		if err := rows.Scan(&gid, &tag_id, &x1); err != nil {
			fmt.Println(err)
		}
		fmt.Println(gid, tag_id, x1)
	}

	if err := db.Ping(); err != nil {
		log.Fatal("PingError: ", err)
	}
}

さて

http://www.kobore.net/diary_techno/?date=20181101

の中の、点(139.9182893339256 36.573831584767085)から半径300メートル以内の全部のノードを、近い順に出せ

を、psqlで出力させると、こんな感じになります。

Golangで実施するとなると、こんなコードになります(select文のコピペで足ります)。

// go get github.com/lib/pq を忘れずに
// go run main5.go
package main

import (
	"database/sql"
	"fmt"
	"log"

	_ "github.com/lib/pq"
)

func main() {
	db, err := sql.Open("postgres",
		"user=postgres password=password host=localhost port=15432 dbname=utsu_db sslmode=disable")
	if err != nil {
		log.Fatal("OpenError: ", err)
	}
	defer db.Close()

	rows, err := db.Query("SELECT source, x1 as longitude, y1 as latitude, ST_Distance('SRID=4326;POINT(139.9182893339256 36.573831584767085)'::GEOGRAPHY, the_geom) as dist FROM ways WHERE ST_DWithin(the_geom, ST_GeographyFromText('SRID=4326;POINT(139.9182893339256 36.573831584767085)'), 300.0) ORDER BY dist")
	if err != nil {
		log.Fatal(err)
	}
	defer rows.Close()

	for rows.Next() {
		var source int
		var longitude float64
		var latitude float64
		var dist float64

		if err := rows.Scan(&source, &longitude, &latitude, &dist); err != nil {
			fmt.Println(err)
		}
		fmt.Println(source, longitude, latitude, dist)
	}

	if err := db.Ping(); err != nil {
		log.Fatal("PingError: ", err)
	}
}

出力結果は同じものになります。

 

2022/04,江端さんの技術メモ

OpenStreetMapのデータからosmiumを使って鉄道だけを抽出しPostGISに入力する方法（宇都宮周辺の鉄道データを使った例）

を、参考にさせて頂いて、試しているのですが、うちのosmiumは、

root@c517fca5dec0:/usr/bin# osmium export -f pg -o utsunomiya-railway-latest.pg utsunomiya-railway-latest.osm.pbf
Set output format with --output-format or -f to 'geojson', 'geojsonseq', or 'text'.

といって、pgファイルを作ってくれません。色々調べたのですが、どうにも上手くいかずに、別の方法を探しています。

で、ちょっと興味本位で、

root@c517fca5dec0:/usr/bin# apt install osm2pgsql

をやってみたら、コンテナの中に、サクッとインストールされたので、これを使う方向で検討しています。

------

>docker container exec -it utsunomiya_db_1 bash

で入ったのですが、osm2pgsqlとかosmiumはインストールできなかったのですが、apt update, apt install, apt upgrade, などを叩き込んで、

>apt install osm2pgsql

>apt install osmium-tool  (×osmium)

>apt install postgresql

を強行しました。

こうしたら、上記のエラー(Set output format with --output-format or -f to 'geojson', 'geojsonseq', or 'text'.)が見えなくなりましたので、osm2pgsqlの方向は見合わせることにしました。

よし、続行だ。

この後、アップデートすると動かなくなるので、アップデートを中途半端に行って(完璧にしない)で、pgファイルを作って、インストールに成功。

その後、上記参照ページのSさんから教えて頂いた

utsu_rail_db=# select geom, tags->>'railway' as railway, tags->>'name' as name from utunomiya where tags->>'name' like '%宇都宮ライトレール';

を実施しました。

| construction | (仮称)宇都宮ライトレール
0102000020E6100000070000000725CCB43D7D61408EA2186A5E474240E9CD4D40387D61406409C61C5F474240C6490625317D61407055230560474240197DBBDB2B7D6140A5AA645B61474240BB1171CE237D614050BA3EBD63474240970E8C721F7D6140A95B2CFB644742
4055FBCFF5187D6140B3171B0467474240

| construction | (仮称)宇都宮ライトレール
0102000020E610000008000000BE5B0F15197D61400686072868474240B3075A81217D6140BA50549165474240419B1C3E297D6140B53D30366347424063AD461B2C7D6140685E697462474240CC7C073F317D6140638DC4156147424009E29755337D6140ECB886BE604742
405A2A6F47387D6140E7CD3C146047424069520ABA3D7D6140B130444E5F474240

| construction | (仮称)宇都宮ライトレール
(74 rows)

と、出てきました。

ーーーようやく、宇都宮GISデータベース目処が付きました。

2022/04,江端さんの技術メモ

「OpenStreetMapのデータから鉄道だけを抽出してGeoJSONで出力する方法」を試してみた件

を前提として、osmから道路と鉄道路線のみを残して残りは消去する方法を探してみました。

どうやら、これでできるようです。

まず、宇都宮地区を抽出します。

osmium extract --bbox 139.76675736787732,36.47004171587971,140.1596765021647,36.659949299138 -o utunomiya.osm.pbf kanto-latest.osm.pbf

次に、以下を行います。

root@2e317961f1d0:/tmp# osmium tags-filter utunomiya.osm.pbf nw/railway nw/highway=motorway,service,primary,secondary,tertiary,unclassified,residential -o utunomiya-lrt-latest.osm.pbf

上記のコマンドのnw(r)は、以下のような意味になっています。

nwr は、OpenStreetMapの node, way, relationsを指定しています。

nw/(キー)=(タグ),(タグ),(タグ),(タグ),(タグ).....

という記述をします。

(キー)は、https://taginfo.openstreetmap.org/keys を参照

(タグ)は、https://taginfo.openstreetmap.org/tags を参照するといいでしょう。

https://wiki.openstreetmap.org/wiki/JA:Tag:highway%3Dresidential の、

上記をクリックすると、motorway,service,primary,secondary,tertiary,unclassified,residential 

の意味が分かります(多分)。

2022/04,江端さんの技術メモ

http://kobore.net/test1031-2.cpp

/////////////////////////////////////////
// 
//
// 適当な2つの座標を指定して、そのルートをトレースする
// 
//  ST_Distance, ST_DWithinの使い方
//



/*

gcc -g -I"D:\PostgreSQL\pg96\include" test1031-2.cpp -o test1031-2.exe -L"D:\PostgreSQL\pg96\lib" -llibpq -lwsock32 -lws2_32

*/

#include <stdio.h>
#include <stdlib.h>
#include <stdint.h>
#include <string.h>
#include <sys/types.h>
#include <math.h>
#include "libpq-fe.h"

struct PERSON{  
  double p_x; // 現在のX座標
  double p_y; // 現在のY座標
};

#define rad2deg(a) ((a)/M_PI * 180.0) /* rad を deg に換算するマクロ関数 */
#define deg2rad(a) ((a)/180.0 * M_PI) /* deg を rad に換算するマクロ関数 */

double distance_km(double a_longitude, 
				   double a_latitude, 
				   double b_longitude, 
				   double b_latitude, 
				   double *rad_up)				   				 
{
  double earth_r = 6378.137;

  double loRe = deg2rad(b_longitude - a_longitude); // 東西  経度は135度
  double laRe = deg2rad(b_latitude - a_latitude); // 南北  緯度は34度39分

  double EWD = cos(deg2rad(a_latitude))*earth_r*loRe; // 東西距離
  double NSD = earth_r*laRe; //南北距離

  double distance = sqrt(pow(NSD,2)+pow(EWD,2));  
  *rad_up = atan2(NSD, EWD);

  return distance;
}


double diff_longitude(double diff_p_x, double latitude) 
{
  double earth_r = 6378.137;
  // ↓ これが正解だけど、
  double loRe = diff_p_x / earth_r / cos(deg2rad(latitude)); // 東西
  // 面倒なので、これで統一しよう(あまり差が出ないしね)
  //double loRe = diff_p_x / earth_r / cos(deg2rad(35.700759)); // 東西
  double diff_lo = rad2deg(loRe); // 東西

  return diff_lo; // 東西
}
  
double diff_latitude(double diff_p_y) 
{
  double earth_r = 6378.137;
  double laRe = diff_p_y / earth_r;  // 南北
  double diff_la = rad2deg(laRe); // 南北
  
  return diff_la; // 南北
}

int main (){

  PERSON test_person;

  // 適当な2つの座標
  double dep_x = 139.46507, dep_y = 35.59577;
  double arr_x = 139.47627, arr_y = 35.60358;

  const char *conninfo = "host=localhost user=postgres password=c-anemone dbname=city_routing";
  
  // データベースとの接続を確立する 
  PGconn *conn = PQconnectdb(conninfo);
  PGresult *res;

  // バックエンドとの接続確立に成功したかを確認する
  if (PQstatus(conn) != CONNECTION_OK){
	fprintf(stderr, "Connection to database failed: %s",
			PQerrorMessage(conn));
  }
  
  /////////////////////// startのもっとも近いノード番号を検索する ////////////////////////
  
  char stringSQL[1024] = {0};  
  
  // 半径300メートルくらいで検索
  sprintf(stringSQL,"SELECT source, ST_Distance('SRID=4326;POINT(%f %f)'::GEOGRAPHY, the_geom) as dist FROM ways WHERE ST_DWithin(the_geom, ST_GeographyFromText('SRID=4326;POINT(%f %f)'), 300.0) ORDER BY dist;", dep_x, dep_y, dep_x, dep_y);
  
  res = PQexec(conn, stringSQL);
	
  if (res == NULL){
	fprintf(stderr, "failed: %s",
			PQerrorMessage(conn));
  }
  
  // SELECTの場合戻り値は PGRES_TUPLES_OK.  これを確認しておく
  if (PQresultStatus(res) != PGRES_TUPLES_OK){
	PQerrorMessage(conn);
  }
  
  int nFields = PQnfields(res);
  
  // そして行を表示する。
  int i_count = 0;
  
  while( atof(PQgetvalue(res, i_count, 1)) == 0.0){
	i_count += 1;
  }
  
  double dist = atof(PQgetvalue(res, i_count, 1));
  int start_source = atoi(PQgetvalue(res, i_count, 0));
  
  PQclear(res); // SQL文の発行に対して必ずクリアする
  
  /////////////////////// startとstartのもっとも近いノードの距離を計測する ////////////////////////  
  
  memset( stringSQL, 0, sizeof(stringSQL)); // 念の為クリア

  sprintf(stringSQL,"SELECT ST_Distance('SRID=4326;POINT(%f %f)'::GEOGRAPHY, the_geom) from ways where source = %d;", dep_x, dep_y, start_source );

  res = PQexec(conn, stringSQL);  

  if (res == NULL){
	fprintf(stderr, "failed: %s",
			PQerrorMessage(conn));
  }
  
  // SELECTの場合戻り値は PGRES_TUPLES_OK.  これを確認しておく
  if (PQresultStatus(res) != PGRES_TUPLES_OK){
	PQerrorMessage(conn);
  }
  
  nFields = PQnfields(res);

  double dep_dist = atof(PQgetvalue(res, 0, 0));  

  PQclear(res); // SQL文の発行に対して必ずクリアする


  ////////////////////// endのもっとも近いノード番号を検索する///////////////////////////

  memset( stringSQL, 0, sizeof(stringSQL)); // 念の為クリア
  
  // 半径300メートルくらいで検索
  sprintf(stringSQL,"SELECT source, ST_Distance('SRID=4326;POINT(%f %f)'::GEOGRAPHY, the_geom) as dist FROM ways WHERE ST_DWithin(the_geom, ST_GeographyFromText('SRID=4326;POINT(%f %f)'), 300.0) ORDER BY dist;", arr_x, arr_y, arr_x, arr_y);

  res = PQexec(conn, stringSQL);
	
  if (res == NULL){
	fprintf(stderr, "failed: %s",
			PQerrorMessage(conn));
  }
  
  // SELECTの場合戻り値は PGRES_TUPLES_OK.  これを確認しておく
  if (PQresultStatus(res) != PGRES_TUPLES_OK){
	PQerrorMessage(conn);
  }
  
  nFields = PQnfields(res);
  
  // そして行を表示する。
  i_count = 0;
  
  while( atof(PQgetvalue(res, i_count, 1)) == 0.0){
	i_count += 1;
  }
  
  dist = atof(PQgetvalue(res, i_count, 1));
  int end_source = atoi(PQgetvalue(res, i_count, 0));
  
  PQclear(res); // SQL文の発行に対して必ずクリアする

  /////////////////////// endとendのもっとも近いノードの距離を計測する ////////////////////////  
  
  memset( stringSQL, 0, sizeof(stringSQL)); // 念の為クリア

  sprintf(stringSQL,"SELECT ST_Distance('SRID=4326;POINT(%f %f)'::GEOGRAPHY, the_geom) from ways where source = %d;", arr_x, arr_y, end_source );

  res = PQexec(conn, stringSQL);  

  if (res == NULL){
	fprintf(stderr, "failed: %s",
			PQerrorMessage(conn));
  }
  
  // SELECTの場合戻り値は PGRES_TUPLES_OK.  これを確認しておく
  if (PQresultStatus(res) != PGRES_TUPLES_OK){
	PQerrorMessage(conn);
  }
  
  nFields = PQnfields(res);

  double arr_dist = atof(PQgetvalue(res, 0, 0));  

  PQclear(res); // SQL文の発行に対して必ずクリアする

  printf("start_source:%d end_source:%d\n",start_source, end_source);
  printf("dep_dist:%f arr_dist:%f\n",dep_dist, arr_dist);


  
  ////////////////// start_sourceとend_sourceのノードをダイクストラで計算する //////////////////////

  printf("%-15f,%-15f\n",dep_x, dep_y); // 出発点を書き込んでおく

  memset( stringSQL, 0, sizeof(stringSQL)); // 念の為クリア

  sprintf(stringSQL, "SELECT node, edge FROM pgr_dijkstra('SELECT gid as id, source, target, cost_s As cost, reverse_cost_s AS reverse_cost FROM ways', %d, %d, true );",start_source, end_source);

  res = PQexec(conn, stringSQL);  

  if (res == NULL){
	fprintf(stderr, "failed: %s",
			PQerrorMessage(conn));
  }
  
  // SELECTの場合戻り値は PGRES_TUPLES_OK.  これを確認しておく
  if (PQresultStatus(res) != PGRES_TUPLES_OK){
	PQerrorMessage(conn);
  }
  
  nFields = PQnfields(res);

  // start_sourceとend_sourceのノードのノード番号の一つをゲットする処理
  for (int i = 0; i < PQntuples(res)-1  ; i++) {  // バス停とバス停の間のノードを全部出す
	  
	/* 
	   ここでは、両端のエッジの座標が必要になる。
	   pgr_dijkstraでは、エッジ(両端)情報(x1,y1,x2,y2)が取得できるのだけど、
	   どっちが先端でどっちが終端かが分からない。
	   
	   そこで(恐しく迂遠だけえど)、まずエッジ情報(x1,y1,x2,y2)を得てから、
	   ノード情報(x1, y1)を取得する(ちなみにノード情報のx2,y2は、
	   交差点などの場合は複数出てくるので、信用してはならない)。
	   
	   で、ノード情報のx1,y1を先端として、エッジ情報のx1,y1とx2,y2と一致していない方を終端とする
	   という処理を取っている。
	   
	   (もっと簡単な方法があったら、誰か教えて下さい)
	   
	*/
	
	double node[2] = {0.0}; // x1, y1
	double edge[4] = {0.0}; // x1, y1, x2, y2
	
	int dummy_int_1 = 0;
	int dummy_int_2 = 0;
	
	for (int j = 0; j < nFields; j++) { // (j=0:node(source)と j=1:edge(gid)の値をゲットする)
	  
	  // まずノードの方
	  if (j == 0){//(j=0:node(source)
		
		memset( stringSQL, 0, sizeof(stringSQL)); // 念の為クリア
		sprintf(stringSQL, "SELECT x1,y1 from ways where source = %s;",
				PQgetvalue(res, i, j)); // ノードの座標を得る
		
		dummy_int_1 = atof(PQgetvalue(res, i, j));
		
		PGresult *res2 = PQexec(conn, stringSQL);
		if (res2 == NULL){
		  fprintf(stderr, "failed: %s",
				  PQerrorMessage(conn));
		}
		
		
		int nFields2 = PQnfields(res2);
		
		for (int j2 = 0; j2 < nFields2; j2++) { // node(source)のx1,y1の2つの値
		  // printf("%-15s", PQgetvalue(res2, 0, j2));
		  node[j2] = atof(PQgetvalue(res2, 0, j2));  
		}
		
		PQclear(res2); // SQL文の発行に対して必ずクリアする
	  }
	  
	  //次にエッジの方
	  if (j == 1){//(j=1:edge(gid)
		
		memset( stringSQL, 0, sizeof(stringSQL)); // 念の為クリア
		
		sprintf(stringSQL, "SELECT x1,y1,x2,y2 from ways where gid = %s;",
				PQgetvalue(res, i, j)); // ノードの座標を得る
		
		dummy_int_2 = atof(PQgetvalue(res, i, j));
		
		PGresult *res2 = PQexec(conn, stringSQL);
		if (res2 == NULL){
		  fprintf(stderr, "failed: %s",
				  PQerrorMessage(conn));
		}
		
		int nFields2 = PQnfields(res2);
		
		for (int j2 = 0; j2 < nFields2; j2++) { // node(source)のx1,y1の2つの値
		  // printf("%-15s", PQgetvalue(res2, 0, j2));
		  edge[j2] = atof(PQgetvalue(res2, 0, j2));  
		}
		
		PQclear(res2); // SQL文の発行に対して必ずクリアする
		
	  }
	}
	
	double start_x, start_y, end_x, end_y;
	
	//出揃ったところで、始点と終点の判定を行う
	if ((fabs(node[0] - edge[0]) < 1e-6) && (fabs(node[1] - edge[1]) < 1e-6)){
	  start_x = edge[0];
	  start_y = edge[1];
	  end_x = edge[2];
	  end_y = edge[3];
	}else{
	  start_x = edge[2];
	  start_y = edge[3];
	  end_x = edge[0];
	  end_y = edge[1];
	}
	
	// printf("両端の確定 %f,%f,%f,%f\n", start_x, start_y, end_x, end_y);
	
	//両端の進行方向(ベクトル)を計算する
	double rad_up1;
	distance_km(start_x, start_y, end_x, end_y, &rad_up1);
	
	// 確定直後にバス位置は、出発点(start_x, start_y)に強制移動する
	test_person.p_x = start_x;
	test_person.p_y = start_y;
	
	// ここから0.1m単位でグルグル回す

	int do_flag = 0;
	do{
	  // printf("x=%-15f,y=%-15f\n",test_person.p_x, test_person.p_y);
	  printf("%-15f,%-15f\n",test_person.p_x, test_person.p_y);
	  
	  // 以下の0.1は1サンプリング0.1メートルを擬制
	  test_person.p_x += diff_longitude(0.1 * cos(rad_up1), test_person.p_y);	
	  test_person.p_y += diff_latitude(0.1 * sin(rad_up1));
	  
	  double rad0 = atan2((end_y - start_y),(end_x - start_x));
	  double rad1 = atan2((end_y - test_person.p_y),(end_x - test_person.p_x));
	  
	  // ここは、http://kobore.net/over90.jpg で解説してある
	  
	  if (fabs(rad0 - rad1) >= 3.141592654 * 0.5){
		// 終点越えの場合、終点に座標を矯正する
		test_person.p_x = end_x;
		test_person.p_y = end_y;
		
		do_flag = 1; // フラグを上げろ
	  }
	  
	}while(do_flag == 0);
	
  } // バス停とバス停の間のノードを全部出す
  
  PQclear(res); // SQL文の発行に対して必ずクリアする
  
  printf("%-15f,%-15f\n",arr_x, arr_y); // 終着点を書き込んでおく

}

上記のプログラムと同じような振舞をするプログラムをGolangで書いてみました。

// go get github.com/lib/pq を忘れずに
// go run main7.go

/*
	前提は、https://ebata-books.booth.pm/items/3351408 の環境

	(1)スタート地点  139.91804656152837 36.57246810341353,
	(2)ゴール地点 139.93515104919825 36.55883778950532
	(3)とした時、この2点に最も近いノードを探し出して、
	(4)その2点をダイクストラ法で、最短距離探索をするプログラム

	分かったこと
	上記(3)のノードが、孤立ノードになっている場合があり、この場合(4)が実施できないので、
	ダイクストラに失敗する場合は、別のノードでトライアルする必要あり
	(このコードでは展開しないけど)
*/

package main

import (
	"database/sql"
	"fmt"
	"log"
	"math"

	_ "github.com/lib/pq"
)

var source int
var longitude float64
var latitude float64
var dist float64

func rad2deg(a float64) float64 {
	return a / math.Pi * 180.0
}

func deg2rad(a float64) float64 {
	return a / 180.0 * math.Pi
}

func distance_km(a_longitude, a_latitude, b_longitude, b_latitude float64) (float64, float64) {
	earth_r := 6378.137

	loRe := deg2rad(b_longitude - a_longitude) // 東西  経度は135度
	laRe := deg2rad(b_latitude - a_latitude)   // 南北  緯度は34度39分

	EWD := math.Cos(deg2rad(a_latitude)) * earth_r * loRe // 東西距離
	NSD := earth_r * laRe                                 //南北距離

	distance_km := math.Sqrt(math.Pow(NSD, 2) + math.Pow(EWD, 2))
	rad_up := math.Atan2(NSD, EWD)

	return distance_km, rad_up
}

func diff_longitude(diff_p_x, latitude float64) float64 {

	earth_r := 6378.137
	// ↓ これが正解だけど、
	loRe := diff_p_x / earth_r / math.Cos(deg2rad(latitude)) // 東西
	// 面倒なので、これで統一しよう(あまり差が出ないしね)
	//loRe := diff_p_x / earth_r / math.Cos(deg2rad(35.700759)) // 東西
	diff_lo := rad2deg(loRe) // 東西

	return diff_lo // 東西
}

func diff_latitude(diff_p_y float64) float64 {
	earth_r := 6378.137
	laRe := diff_p_y / earth_r // 南北
	diff_la := rad2deg(laRe)   // 南北

	return diff_la // 南北
}

func sub_main() {

	x1, y1 := 139.91804656152837, 36.57246810341353
	x2, y2 := 139.93515104919825, 36.55883778950532

	d_km, rad_up := distance_km(x1, y1, x2, y2)

	//x1 += diff_longitude(0.00566*math.Cos(rad_up), y1)
	//y1 += diff_latitude(0.00566 * math.Sin(rad_up))

	x1 += diff_longitude(1.0*math.Cos(rad_up), y1)
	y1 += diff_latitude(1.0 * math.Sin(rad_up))

	fmt.Println(d_km, rad_up)
	fmt.Println(x1, y1)

	d_km, rad_up = distance_km(x1, y1, x2, y2)

	//x1 += diff_longitude(0.00566*math.Cos(rad_up), y1)
	//y1 += diff_latitude(0.00566 * math.Sin(rad_up))

	x1 += diff_longitude(1.0*math.Cos(rad_up), y1)
	y1 += diff_latitude(1.0 * math.Sin(rad_up))

	fmt.Println(d_km, rad_up)
	fmt.Println(x1, y1)

}

func main() {
	db, err := sql.Open("postgres",
		"user=postgres password=password host=localhost port=15432 dbname=utsu_db sslmode=disable")
	if err != nil {
		log.Fatal("OpenError: ", err)
	}
	defer db.Close()

	// スタート地点  139.91804656152837 36.57246810341353
	rows, err := db.Query("SELECT source, x1 as longitude, y1 as latitude, ST_Distance('SRID=4326;POINT(139.91804656152837 36.57246810341353)'::GEOGRAPHY, the_geom) as dist FROM ways WHERE ST_DWithin(the_geom, ST_GeographyFromText('SRID=4326;POINT(139.91804656152837 36.57246810341353)'), 300.0) ORDER BY dist LIMIT 1")
	if err != nil {
		log.Fatal(err)
	}
	defer rows.Close()

	for rows.Next() {

		if err := rows.Scan(&source, &longitude, &latitude, &dist); err != nil {
			fmt.Println(err)
		}
		fmt.Println(source, longitude, latitude, dist)
	}

	o_source := source

	// ゴール地点 139.93515104919825 36.55883778950532
	rows, err = db.Query("SELECT source, x1 as longitude, y1 as latitude, ST_Distance('SRID=4326;POINT(139.93515104919825 36.55883778950532)'::GEOGRAPHY, the_geom) as dist FROM ways WHERE ST_DWithin(the_geom, ST_GeographyFromText('SRID=4326;POINT(139.93515104919825 36.55883778950532)'), 300.0) ORDER BY dist LIMIT 1")
	if err != nil {
		log.Fatal(err)
	}
	defer rows.Close()

	for rows.Next() {

		if err := rows.Scan(&source, &longitude, &latitude, &dist); err != nil {
			fmt.Println(err)
		}
		fmt.Println(source, longitude, latitude, dist)
	}

	d_source := source

	fmt.Println(o_source, d_source)

	//rows, err = db.Query("select x1,y1 from ways where source in (SELECT node as source FROM pgr_dijkstra('SELECT gid as id, source, target,length_m as cost FROM ways', $1::bigint , $2::bigint , false))", o_source, d_source)
	rows, err = db.Query("SELECT node as source FROM pgr_dijkstra('SELECT gid as id, source, target,length_m as cost FROM ways', $1::bigint , $2::bigint , false)", o_source, d_source)
	if err != nil {
		log.Fatal(err)
	}
	defer rows.Close()

	fmt.Println("2")

	x1, y1 := -1.0, -1.0
	_x1, _y1, _x2, _y2 := -1.0, -1.0, -1.0, -1.0
	px, py := -1.0, -1.0
	flag := 0
	f_flag := 0

	for rows.Next() {

		if err := rows.Scan(&source); err != nil {
			fmt.Println(err)
		}

		rows2, err := db.Query("SELECT source, x1, y1 from ways where source = $1::bigint limit 1", source)
		if err != nil {
			log.Fatal(err)
		}

		for rows2.Next() {

			if f_flag == 0 { //初回だけ2居合入力
				if err := rows2.Scan(&source, &x1, &y1); err != nil {
					fmt.Println(err)
				}
				_x1, _y1 = x1, y1
				//fmt.Println(_x1, _y1)
				f_flag = 1
				continue
			}

			if err := rows2.Scan(&source, &x1, &y1); err != nil {
				fmt.Println(err)
			}
			_x2, _y2 = x1, y1

			//fmt.Println(_x2, ",", _y2)
			_, rad_up := distance_km(_x1, _y1, _x2, _y2)
			//fmt.Println("_x1,_y1,_x2,_y2, rad_up:", _x1, _y1, _x2, _y2, rad_up)

			px, py = _x1, _y1

			for {

				//平均時速20km サンプリング時間は、3.6秒  平均速度は 5.56m 20m/3.6秒
				// 悩むところだな、サンプリング時間を1秒にするか3.6秒にするか → 1秒にするか。
				// 距離端数は切り捨てにして、実質は、平均時速を20km以下としよう
				// (カーブでは減速する、と)

				px += diff_longitude(0.00556*math.Cos(rad_up), py)
				py += diff_latitude(0.00556 * math.Sin(rad_up))

				//double rad0 = atan2((end_y - start_y),(end_x - start_x));
				//double rad1 = atan2((end_y - test_person.p_y),(end_x - test_person.p_x));

				rad0 := math.Atan2((_y2 - _y1), (_x2 - _x1))
				rad1 := math.Atan2((_y2 - py), (_x2 - px))
				// ここは、http://kobore.net/over90.jpg で解説してある

				if math.Abs(rad0-rad1) >= math.Pi*0.5 {
					// 終点越えの場合、終点に座標を矯正する
					px, py = _x2, _y2
					flag = 1 // フラグを上げろ
				}

				fmt.Println(px, ",", py)

				if flag == 1 {
					flag = 0
					_x1, _y1 = _x2, _y2
					break
				}
			}
		}
	}

	//平均時速20km サンプリング時間は、3.6秒  平均速度は 5.56m 20m/3.6秒
	// 悩むところだな、サンプリング時間を1秒にするか3.6秒にするか → 1秒にするか。
	// 距離端数は切り捨てにして、実質は、平均時速を20km以下としよう

	if err := db.Ping(); err != nil {
		log.Fatal("PingError: ", err)
	}
}

2022/04,江端さんの技術メモ

# 誤差 移動 直行 90度 越えた 垂直

2022/04,江端さんの技術メモ

「OpenStreetMapから鉄道路線データを落とせない」と愚痴っていたら、社内の人から「OSMに鉄道データ入っていますよ」と言われて、愕然とした件

という訳で、今日は、OpenStreetMapX.jl をダウンロードして色々やっていたのですが、よく分からりませんでした。というか、何ができるのか、が今一つ肚に落ちてこないのです。

で、こちらの路線は諦めて、とにかく鉄道データが取れる手段を探しました。

OpenStreetMapのデータから鉄道だけを抽出してGeoJSONで出力する方法

を試すことにしました。
から、osm.pbfファイルを落してきました。

osmium を使う方法がでていたのですが、これが私のWindows10環境では実現できませんでした。(正確に言うと、私の環境には"cl.exe"というコマンドがないので、コンパイルができない。Visual Studio C++(?)を入れるとできるらしいのですが、これ以上、パソコンの環境を汚したくなったのです)。
で、色々試みたのですが(この間5～6時間)、上手くいかず、dockerでないかなにないかなーと探していたら、見付けました。
1つ目のWindowsのコマンドプロンプト(MSYS64では上手く動かない)かったので、
C:\Users\ebata>docker run -t -i stefda/osmium-tool /bin/bash
で、いきなりDockerのコンテナの中に入りました。

/bin , /usr/bin, /usr/local/bin のどこかで当るだろうと探りを入れていたら、

root@e87f198ccc07:/usr/bin# ls o*
objcopy objdump od odbcinst ogdi-config openssl osage osmium
root@e87f198ccc07:/usr/bin# osmium
Usage: osmium COMMAND [ARG...]

/usr/binで当たりました。

で、作業用のディレクトリを探していたら、/tmpがあったので、ここにファイルを持ち込んで作業することにしました

2つ目の、別のWindowsのコマンドプロンプトを上げて、今動かしているコンテナを調べました

C:\Users\ebata>docker ps
CONTAINER ID IMAGE COMMAND CREATED STATUS PORTS NAMES
e87f198ccc07 stefda/osmium-tool "/bin/bash" 12 minutes ago Up 12 minutes wonderful_lumiere

で、ローカルからDockerコンテナの/tmpに直接コピーしました。

C:\Users\ebata>docker cp kanto-latest.osm.pbf wonderful_lumiere:/tmp

一つ目のコマンドプロンプトに戻って、ファイルが可能されているか確認しました。

root@e87f198ccc07:/# cd tmp
root@e87f198ccc07:/tmp# ls
kanto-latest.osm.pbf

ちゃんと入っていました。パスも通っているようですので、

root@e87f198ccc07:/tmp# which osmium
/usr/bin/osmium

そのまま、/tmpの中で作業しました

root@e87f198ccc07:/tmp# osmium tags-filter kanto-latest.osm.pbf w/railway -o kanto-railway-latest.osm.pbf
[======================================================================] 100%
root@e87f198ccc07:/tmp# ls
kanto-latest.osm.pbf kanto-railway-latest.osm.pbf
root@e87f198ccc07:/tmp# osmium export kanto-railway-latest.osm.pbf -o kanto-railway-latest.json
root@e87f198ccc07:/tmp# ls
kanto-latest.osm.pbf kanto-railway-latest.json kanto-railway-latest.osm.pbf

2つ目のコマンドプロンプトを立ち上げて、dockerコンテナから、ローカルに、完成したjsonファイルを送り出します。

C:\Users\ebata\docker-osmium-tool>docker cp wonderful_lumiere:/tmp/kanto-railway-latest.json .

ローカルの方に変換されたファイルが戻ってきました。
C:\Users\ebata>ls
 Dockerfile   README.md   kanto-latest.osm.pbf   kanto-railway-latest.json   work
さて、これを、ドラッグして、QGISに放り込んで出てきた図です。

見事に鉄道だけが抽出されています。

しかし、重要なのはそこではありません。

芳賀・宇都宮LRT が、すでに表示されている ―― これは、今の私にとって、とてつもないインパクトなのです。

さて、今回はLRTを中心とした宇都宮市をカバーする地域の交通機関を切り出します。

【食べログ作れる? 】 OpenStreetMap から Osmiumで飲食店の位置情報を取得してGeoJSONに出力してみたを参照して、

$ osmium extract --bbox 左上の経度,左上の緯度,右下の経度, 右下の緯度 -o 出力したいファイル名.pbf 元のファイル名.osm.pbf

から、

$ osmium extract --bbox 139.76675736787732, 36.659949299138, 140.1596765021647, 36.47004171587971 -o utunomiya-lrt.osm.pbf kanto-latest.osm.pbf

として切り出してみます。

root@80e1a2d0dc02:/tmp# osmium extract --bbox 139.76675736787732, 36.659949299138, 140.1596765021647, 36.47004171587971 -o utunomiya-lrt.osm.pbf kanto-latest.osm.pbf
Error parsing command line: too many positional options have been specified on the command line

あれ、コンマと空白がダメだったかな？

root@80e1a2d0dc02:/tmp# osmium extract --bbox 139.76675736787732,36.659949299138,140.1596765021647,36.47004171587971 -o
utunomiya-lrt.osm.pbf kanto-latest.osm.pbf
Need LEFT < RIGHT and BOTTOM < TOP in --box/-b option.

うむ、 では、36.659949299138 と 36.47004171587971 を入れかえて、再度挑戦

root@80e1a2d0dc02:/tmp# osmium extract --bbox 139.76675736787732,36.47004171587971,140.1596765021647,36.659949299138 -o utunomiya-lrt.osm.pbf kanto-latest.osm.pbf
[======================================================================] 100%

成功したっぽいです。

root@80e1a2d0dc02:/tmp# osmium tags-filter utunomiya-lrt.osm.pbf w/railway -o utunomiya-railway-latest.osm.pbf
[======================================================================] 100%

root@80e1a2d0dc02:/tmp# osmium export utunomiya-railway-latest.osm.pbf -o utunomiya-railway-latest.json

では、ここから、2つ目のコマンドプロンプトから、

C:\Users\ebata\docker-osmium-tool>docker cp serene_shirley:/tmp/utunomiya-railway-latest.json .
C:\Users\ebata\docker-osmium-tool>docker cp serene_shirley:/tmp/utunomiya-lrt.osm.pbf .

としました。

utunomiya-railway-latest.json の表示結果はこちらです。

データ量も、関東全域と比較すると、11%程度になりました。

ところで、osmファイルと、osm.pbfファイルの違いって何だろう？

 

2022/04,江端さんの技術メモ

「OpenStreetMapから鉄道路線データを落とせない」と愚痴っていたら、社内の人から「OSMに鉄道データ入っていますよ」と言われて、愕然としました。

江端:「私は、鉄道のOSMデータの取り方に辿り付けていません。いつも、ここからOSMデータをダウンロードしていますが、地理情報しか取れていない(ように思えます)。お暇な時に御教授頂けましたら幸いです。」

同僚:「そこから取れるデータの中に、鉄道路線のデータも入っています」

江端: ―― はい？

江端:「かなりずうずうしいお願いをしていることは分かっているのです、このタグをPostGISのpostgresqlに落す方法とかご存知だったりしますか？ 多分、この辺(Configuration file = /usr/local/share/osm2pgrouting/mapconfig_for_cars.xml)です。」

同僚:「私は、JuliaのOpenStreetMapX.jlを少し改造して、直接OpenStreetMapのosmファイルをパースしていました」

https://github.com/pszufe/OpenStreetMapX.jl

で、現在、readme.mdを読んでいます。

  • Open Street Mapデータの空間解析・シミュレーション・可視化用パッケージ(プロット機能は別パッケージで提供されます)
  • このパッケージの目的は、都市のマルチエージェントモデリングとシミュレーションのためのバックボーンを提供することです。
  • このパッケージは *.osm と *.pbf (@blegat によって寄贈) ファイルを解析し、メタデータに沿った Graphs.jl 表現を生成することができます。

うむ、この連休の使い方としては、有意義な時間を過せそうです。

ところで、Juliaって、コンピュータ言語かな。これ以上、言語は、覚えたくないんだけど ―― 仕方ないですね。

2022/04,江端さんの技術メモ

以前、こちらで、redisのブロードキャスト(Pub/Sub)の方法について2つほど紹介しましたが、試した結果、こっちの方が安定して調子が良くて、現在、こちら(redigo)を使っています

redisを前提として、golangでPubSubを実現するプログラム

Redigoを使う(6)　パブリッシュ／サブスクライブ

で紹介されていた、サンプルプログラムを使って、構造体のデータを丸ごと送信しようとしたのですが、Golangの厳しい型チェックに掴まって、なかなか上手く動かすことができません。

ただ、構造体をJSON形式にすれば、成功することは分かっています。

golang内でredis経由でJSONを飛す時、golangから直接JavaScriptへJSONを飛す時の覚え書き

しかし、今の段階で構造体をJSONに変更すると、その影響が、プログラム全体に波及し、作業が膨大になるので、これはしたくありませんでした。

もしかしたら、構造体のブロードキャストも、JSONの時と同じように、"json.Unmarshal"、 "json.Marshal" を使えばいけるかな？ と思ってやってみたら、あっさりと成功しました。

上記の記事のサンプルプログラムを、構造体データ送付用に改造したものを開示しておきます。

■パブリッシャ(発行元)側

// go run pub.go
// goga\1-9-6\others\pub.go
package main

import (
	"encoding/json"
	"fmt"

	"github.com/gomodule/redigo/redis"
)

type Ch5_info struct {
	Bus_num int     // バスの番号
	CC      int     //  1: 座標情報 2: 停車位置情報
	Present int     // 停車位置番号
	Lat     float64 //
	Lon     float64 //
}

func main() {
	// 接続
	conn, err := redis.Dial("tcp", "localhost:6379")
	if err != nil {
		panic(err)
	}
	defer conn.Close()

	c5i := new(Ch5_info)

	c5i.Bus_num = 1
	c5i.CC = 2
	c5i.Present = 23
	c5i.Lat = 12.34
	c5i.Lon = 56.78

	json_c5i, _ := json.Marshal(c5i)

	// パブリッシュ
	r, err := redis.Int(conn.Do("PUBLISH", "channel_1", json_c5i))
	if err != nil {
		panic(err)
	}
	fmt.Println(r)
}

■サブスクライブ(購読者)側

// go run sub.go
// goga\1-9-6\others\sub.go

package main

import (
	"encoding/json"
	"fmt"

	"github.com/gomodule/redigo/redis"
)

type Ch5_info struct {
	Bus_num int     // バスの番号
	CC      int     //  1: 座標情報 2: 停車位置情報
	Present int     // 停車位置番号
	Lat     float64 //
	Lon     float64 //
}

func main() {
	// 接続
	conn, err := redis.Dial("tcp", "localhost:6379")
	if err != nil {
		panic(err)
	}
	defer conn.Close()

	psc := redis.PubSubConn{Conn: conn}
	psc.Subscribe("channel_1", "channel_2", "channel_3")
	for {
		switch v := psc.Receive().(type) {
		case redis.Message:
			fmt.Printf("%s: message: %s\n", v.Channel, v.Data)

			c5i := new(Ch5_info)

			_ = json.Unmarshal(v.Data, &c5i)

			// 試しに2つほど出力してみる
			fmt.Println(c5i.Bus_num)
			fmt.Println(c5i.Lon)

		case redis.Subscription:
			fmt.Printf("%s: %s %d\n", v.Channel, v.Kind, v.Count)
		case error:
			return
		}
	}
}

私のような悩み持っている人、世界中に沢山いました。
Golangは、C/C++みたいに、自己責任で自由にキャストできないんですよね。
memset()を使って、大量の配列の変数を一気に変更させてしまう、ということもできないようで、Golang不便だなぁ、って思います。

2022/04,江端さんの技術メモ

Golang で複数のgoroutineで、1つのchannelを使い回す場合(要求応答通信)をする際の注意点について

■Case 1

package main

import (
	"fmt"
	"sync"
)

//var Ch1 chan interface{}
var Ch1 chan interface{}

var mutex sync.Mutex
var wg sync.WaitGroup

func funcA(num int) {
	fmt.Println("Staring funcA No.", num)
	for i := 0; i < 10; i++ {
		mutex.Lock()
		fmt.Println("funcA:", num, " send", i+num*10)
		Ch1 <- i + num*10
		k := <-Ch1
		fmt.Println("				returned:", k)
		mutex.Unlock()
	}
	wg.Done()
}

func echo(num int) {
	fmt.Println("Staring echo No.", num)

	for {
		i := <-Ch1
		Ch1 <- i
	}

}

func main() {

	Ch1 = make(chan interface{}) // チャネルの初期化

	wg.Add(1)
	go funcA(1)

	wg.Add(1)
	go funcA(2)

	wg.Add(1)
	go echo(1)

	wg.Wait()
}

結果

funcA: 1  send 15
                                returned: 15
funcA: 1  send 16
                                returned: 16
funcA: 1  send 17
                                returned: 17
funcA: 1  send 18
                                returned: 18
funcA: 1  send 19
                                returned: 19
fatal error: all goroutines are asleep - deadlock!

goroutine 1 [semacquire]:
sync.runtime_Semacquire(0x0)
        c:/go/src/runtime/sema.go:56 +0x25
sync.(*WaitGroup).Wait(0x483fc0)
        c:/go/src/sync/waitgroup.go:130 +0x71
main.main()
        C:/Users/ebata/goga/0-14/main.go:52 +0x131

goroutine 8 [chan receive]:
main.echo(0x0)
        C:/Users/ebata/goga/0-14/main.go:32 +0x8a
created by main.main
        C:/Users/ebata/goga/0-14/main.go:50 +0x125
exit status 2

■Case 2

# 最後に、close(Ch1)を入れてみた

        wg.Add(1)
	go echo(1)

	wg.Wait()
	close(Ch1)
}

結果

funcA: 2 send 25
                                returned: 25
funcA: 2  send 26
                                returned: 26
funcA: 2  send 27
                                returned: 27
funcA: 2  send 28
                                returned: 28
funcA: 2  send 29
                                returned: 29
fatal error: all goroutines are asleep - deadlock!

goroutine 1 [semacquire]:
sync.runtime_Semacquire(0x0)
        c:/go/src/runtime/sema.go:56 +0x25
sync.(*WaitGroup).Wait(0x1e3fc0)
        c:/go/src/sync/waitgroup.go:130 +0x71
main.main()
        C:/Users/ebata/goga/0-14/main.go:52 +0x131

goroutine 8 [chan receive]:
main.echo(0x0)
        C:/Users/ebata/goga/0-14/main.go:32 +0x8a
created by main.main
        C:/Users/ebata/goga/0-14/main.go:50 +0x125
exit status 2

■Case 3

//wg.Add(1)
	go echo(1)

	wg.Wait()
	close(Ch1)
}

結果

funcA: 2  send 29
                                returned: 29
funcA: 1  send 17
                                returned: 17
funcA: 1  send 18
                                returned: 18
funcA: 1  send 19
                                returned: 19

C:\Users\ebata\goga\0-14>

結論

(1)channelを応答通信で使いたいのであれば、送信側(funcA)の中身をMutexでロックする必要がある。受信側は、一つだけならロックの必要はない
(2)送信側のgoroutineが消えると、受信側(echo)が『相手がいなくて寂しくなって』エラーを吐くので、chose(Ch1)で、きっちりchannelを殺しておくこと
(3)受信側(echo)が無限待ち"for()"であるなら、wg.Add(1)を発行してはならない。なぜならechoは、原則として終了しないgoroutineだから。

―― しかし、そこまで構わなくてもいいんじゃなかな、と思う。Golangは、相当"過保護"に作られていると思う。

以上