WebGLを使って画像を正距方位図法で投影変換する処理をまとめたので、そのメモ。

概要

以下のような処理。

緯度経度座標系上の矩形に対応した画像を元に正距方位図法に投影変換して表示。
画像の範囲は任意に指定可能。
複数枚の分割された画像にも対応。
投影変換の中心位置及び表示範囲も任意に指定可能。
楕円体ではなく、完全な球体のモデルでの正距方位図法の投影変換を使用。

表示例

表示例は以下のとおり。
https://tomosn.github.io/simple-raster-projection/map-aeqd-v1.html

ソースコード

GitHubレポジトリ

ソースコード全体は以下のGitHubレポジトリを参照。
WebGLのshaderはhtml中に記述。
simple-raster-projection

座標系

内部では以下のように３段階の座標系を定義して使用している。

緯度経度座標系
データの座標系。画像の領域はこの座標系上で定義する。
経度は $[-180,+180]$、緯度は $[-90,+90]$ の範囲とする（単位が度の場合）。
XY座標系
緯度経度の座標を投影変換した後の座標系。
画面の表示範囲はこの座標系上の矩形で設定する。
正距方位図法の場合は $[(-\pi,-\pi)-(+\pi,+\pi)]$ の範囲をとる。但し、この値は投影図法によって異なる。
画面座標系
画面（スクリーン）に対応した座標系。

shader

shaderの内容は以下。いずれも map-aeqd-v1.html からの抜粋。

vertex shader

precision highp float;
attribute float aCoordX;
attribute float aCoordY;
attribute vec2 aViewCoord;     //  投影後のView上のXY座標系での範囲に対応
varying vec2 vViewCoord;

void main() {
  gl_Position = vec4(aCoordX, aCoordY, 1.0, 1.0);
  vViewCoord = aViewCoord;
}

vertex shaderではたいしたことはおこなっていない。
aViewCoordはXY座標系上の表示範囲に対応する。

fragment shader

precision highp float;
varying vec2 vViewCoord;    //  表示中のView上の座標値
uniform vec2 uDataCoord1;   //  画像データのSouthWestの座標値
uniform vec2 uDataCoord2;   //  画像データのNorthEastの座標値
uniform vec2 uProjCenter;   //  投影中心
uniform sampler2D uTexture;

const float pi = 3.14159265;
const float epsilon = 0.00000001;
const float xyRadius = pi;

//   AEQD逆変換
vec2 proj_inverse(vec2 center, vec2 xy) {
  float sinPhi0 = sin(center.y);
  float cosPhi0 = cos(center.y);

  float c_rh = length(xy);

  if ( c_rh < epsilon ) {
    return center;
  }
  if ( c_rh - epsilon > pi ) {
    c_rh = pi;
  }
  float cos_c = cos(c_rh);
  float sin_c = sin(c_rh);

  float phi = asin( clamp( cos_c * sinPhi0 + xy.y * sin_c * cosPhi0 / c_rh, -1.0, 1.0 ) );
  float lam = mod( center.x + atan( xy.x * sin_c, c_rh * cosPhi0 * cos_c - xy.y * sinPhi0 * sin_c ) + pi, 2.0 * pi ) - pi;

  return vec2(lam, phi);
}

//  XY座標の値が範囲内か否か
float validate_xy(vec2 xy) {
  return 1.0 - step( xyRadius, length(xy) );
}

void main() {
  float inXY = validate_xy(vViewCoord);
  bool invalid = true;
  if ( 0.0 < inXY ) {
    vec2 lp = proj_inverse(uProjCenter, vViewCoord);                 //  緯度経度
    vec2 ts = (lp - uDataCoord1) / (uDataCoord2 - uDataCoord1);      //  変換元画像上の座標
    if ( 0.0 <= ts.x && 0.0 <= ts.y && ts.x <= 1.0 && ts.y <= 1.0) {
      gl_FragColor = texture2D(uTexture, ts);
      invalid = false;
    }
  }
  if (invalid) {
    discard;
  }
}

fragment shaderが投影変換の主要な処理を担っている。
uDataCoord1、uDataCoord2が画像データの緯度経度の範囲に対応する（shader中では単位としてradianを使用）。uProjCenterが投影中心の座標に対応する。

画像の投影変換は正距方位図法の逆変換を使用しており、それが proj_inverse(..)の関数にあたる。また、 validate_xy(..) の関数で範囲外の部分（投影後の円の外側かどうか）を判定している。

他の投影図法

正距方位図法と似たような投影変換の場合、fragment shaderの proj_inverse(..) と validate_xy(..) を変更すれば対応できる。
例えばランベルト正積方位図法の場合、fragment shader中の xyRadius と proj_inverse(..) を以下のように変更するだけで対応可能。

fragment shaderの抜粋

const float xyRadius = 2.0;

//   LAEA逆変換
vec2 proj_inverse(vec2 center, vec2 xy)
{
  float sinPhi0 = sin(center.y);
  float cosPhi0 = cos(center.y);

  float rho = length(xy);

  if ( rho < epsilon ) {
    return center;
  }
  if ( rho - epsilon > xyRadius ) {
    rho = xyRadius;
  }

  float c_rh = 2.0 * asin( clamp( rho/2.0, -1.0, 1.0) );

  float cos_c = cos(c_rh);
  float sin_c = sin(c_rh);

  float phi = asin( clamp( cos_c * sinPhi0 + xy.y * sin_c * cosPhi0 / rho, -1.0, 1.0 ) );
  float lam = mod( center.x + atan( xy.x * sin_c, rho * cosPhi0 * cos_c - xy.y * sinPhi0 * sin_c ) + pi, 2.0 * pi ) - pi;

  return vec2(lam, phi);
}

その他

WebGLはまだまだ理解が足りてないので、おかしなところがあるかもしれない。
緯線・経線の描画については、点列を生成して、正方向の変換をかけてプロットして繋いでいけば可能。但し、以下のような難点あり。
滑らかな曲線を描画しようとすると密な点列を生成する必要がある。
外周部（対蹠点）付近では、いくら点列を密に生成したところで大きな飛びが発生してしまう。
上記の結果、滑らかな緯線・経線を描画しようとすると、処理が結構重くなる。
上記の方法での緯線・経線の描画は既に別プログラムで実現済みだが、抜本的に方式を変更することを考え中。