位置情報とパブリックデータ#05: (地理空間関数)バッファ範囲の特定
(TeradataのSQL実行環境に関してはこちらをご覧ください)
TeradataのGeospatial関数を引き続きご紹介。以降でst_buffer()という位置オブジェクトのバッファ(周辺)を計算する関数と、st_length()というlinestringオブジェクトの距離を計算する関数をご紹介しますが、これら関数の入出力単位は度数です。度数は60進法(度分秒)で、日本のどこにいるかでも緯度1秒あたりの距離、経度1秒あたりの距離が異なってきます。大雑把に把握するという前提で、以下のサイトにてまとめてくださってますが、1秒あたり経度方向で25m、緯度方向で31mが大まかな基準です。度数に換算するなら、3600秒を掛けた値、つまり1度は3600 * 25メートル、3600 * 31メートルということになります。
日本付近で角度一秒とは約、緯度は31m、経度は25mと覚えておきましょう
そしてこれをさらに10進法に変換します。単純に、1メートルあたりと考えると、1で割ればよいため、1メートルの入力の場合、以下の通りになります。500メートルなら分子を500とします。
経度: 1/(3600 * 25)
緯度: 1/(3600 * 36)
もう一つ考慮しなければならない点があります。関数の入力は1つの値であるため、上記の経度を採用するか、緯度を採用するかを決めなければなりません。明らかに平行に伸びる線のバッファであれば経度を、縦に伸びる線のバッファであれば緯度を採用すれば現実に近くなるでしょう。でも円やポリゴンの場合は中庸の1/(28*3600)とするか、厳しめに緯度の方を採用するか、緩めに経度の方を採用するかを決める必要があります。また、実務において継続的にこれら距離指標を見ていくのであれば、基準としてどういった計算方法を採用するのかを統一した方が、分析を実施するチームにおいて混乱もなくなりますし、分析結果を共有する先にも説明がしやすいです。なによりも、いくつかの分析テーマを実施した際に結果を同じ前提でみて比較することができます。ということでノートブックを以下に。データはマガジン内以前の記事にて利用したデータを再利用しています。
出力した結果を描画したのが以下です。若干縦長に潰れた形でバッファを作成しているのがわかります。
ライン距離の計算
ラインの距離計算に関しても単位が度数であるため、こちらは出てきた結果に対して追加の計算をしてあげる必要があります。
経度: 出力結果 * (25 * 3600)
緯度: 出力結果 * (31 * 3600)
中庸: 出力結果 * (28 * 3600)
これで大まかな計算は可能ですが、より正確に計算する必要がある場合は、linestringの座標をいったん複数のpointに分解して、それぞれに対してst_spheroidaldistance()で距離を計算し、足し合わせればよいかと思います。点ABCで構成されるlinestringであれば、AとB、BとCの距離を計算するイメージです。
もう一つのバッファ計算方法
より正確にバッファを計算したいのであれば、ちょっと面倒ですが別のアプローチを採ります。おおまかなアプローチとしては、まず特定点からの距離を北方向と東方向にとることができる関数であるst_spheroidalbuffermbr()を用います。この関数は地球の楕円球体を前提に、特定点から入力距離(メートル)に基づいた座標位置に相当する値を4つ計算してくれます。4つというのがx軸の最大値最小値、そしてy軸の最大値最小値です。この4つの値は、任意のポリゴンを囲むことができる最小の四角形を定義できることから、MBR: Minimum Bounding Rectangle (最小境界矩形)と呼ばれるそうです。st_spheroidalbuffermbr()は、楕円球体を前提に、バッファを作成し、そのポリゴンにおける4つの値を出力する関数ということになります。点を中心点としたとき、そこからのx軸、y軸における半径距離での座標が決まると、角度をずらしていくことによって中心点を中心とした同一距離のバッファを意味する座標上の点を計算することが可能です。この点を限りなくたくさん出力してつなげれば、バッファ円が作成可能となり、それをつなぎ合わせれば、線や面のバッファポリゴンも作成可能となります。今回は比較的シンプルに、36個の座標を出力して、バッファ円を作成してみます。
なお、楕円上の座標計算方法ですが、最終的には三角関数を用います。
st_spheroidalbuffermbr()において入力値は2つで、頭につける点座標と、()内にセットする距離の2つがありますが、ここに上述の値が入ります。x軸とy軸それぞれ計算し、x軸が経度、y軸が緯度を示します。そして、ここで出力された値から、x軸、y軸の最大値を取得します。これも関数が用意されており、xmax()とymax()です。ここまでを合わせると、以下の計算になります。
([st_geometry型の基準点ポイント座標データ].st_spheroidalbuffermbr(Nメートル * (25 * 3600))).xmax() as 指定距離におけるx軸の経路位置(真東)
([st_geometry型の基準点ポイント座標データ].st_spheroidalbuffermbr(Nメートル * (31 * 3600))).ymax() as 指定距離におけるy軸の経路位置(真北)
これで真東と真北における座標位置が10進法度数ベースで求められました。度数ベースでの中心点からの距離は、この値から中心点の座標を引くことによって得られます。
指定距離におけるx軸の経路位置(真東) - 基準点ポイント座標のx位置 as 経度ベースでの距離
指定距離におけるy軸の経路位置(真北) - 基準点ポイント座標のy位置 as 緯度ベースでの距離
そして任意の角度における楕円の位置は、以下にて計算できます。
経度ベースでの距離 * コサイン(ほしい角度 * π/180) + 基準点ポイント座標のx位置 as ほしい角度での座標値
経度ベースでの距離 * サイン(ほしい角度 * π/180) + 基準点ポイント座標のy位置 as ほしい角度での座標値
ここらへんの詳細は以下を:
コサインとサインの関数における入力はラジアン角度なので、パイ(π, 円周率)をかけて180で割る処理をして、コサインの関数に突っ込みます。緯度もサインに対して同じくです。円周率は、Teradataにはほかのデータベースでよくあるpi()がないため、近似値をセットします。これを360度方向に角度を変えながら計算して楕円上の座標を取得します。4つしか座標がなければダイヤモンド型ですが、座標の数を増やせば、形は円に近づいていきます。今回は10度ずつ角度をずらし、36座標取得することにします。ここまでを一式のSQLにしました。結構長くなりSQLも醜いですが、考え方さえ理解してしまえば、あとは盲目的にデータを入力すれば出力してくれます。ノートブックを以下に。データは先ほど同様、マガジン内以前の記事にて利用したデータを再利用しています。
この後、結果をグラフに出す際にはポリゴンデータを分解する必要がありますが、ポリゴンを構成する複数座標を分解抽出して数値型の列に出力する例を別途ご紹介したいので、それについてはまた改めて別の機会に。グラフ描画だけ以下に。前のst_buffer()での結果と合わせて描画した結果ですが、比較的縦横比がきれいな円でバッファを作成してくれることがわかります。
以降に、上述した2つのグラフを描画するために利用した定義用jsonファイルとデータのcsv/topojsonファイルを置いておきます。
以上です。
#geospatial #gis #teradata #sql #analytics #st_geometry