メルボルンの小さな地域の緯度/経度の値があります。-37.803134,145.132377 と、openstreet マップ (Osmarender Image ) からエクスポートしたフラット イメージ。画像の幅：1018、高さ：916
C++ を使用して、緯度/経度をポイントが位置を反映する X、Y 座標に変換できるようにしたいと考えています。
以下に示すようなWebで見つけたさまざまな式を使用しましたが、何も役に立ちません。
var y = ((-1 * lat) + 90) * (MAP_HEIGHT / 180);
var x = (lon + 180) * (MAP_WIDTH / 360);
これを行う方法について明確な説明があれば、非常に役立ちます。どんなコードでも大歓迎です。
これを行うには、単一の緯度/経度ペアよりも多くの情報が必要です。
この段階では、提供された情報には次の 2 つの情報が不足しています。
また、画像が南北に配置されていると仮定します。たとえば、北が左上隅を向いていません。それは物事を複雑にする傾向があります。
最も簡単な方法は、(0, 0) ピクセルと (1017, 915) ピクセルに対応する緯度/経度座標を正確に把握することです。次に、補間により、指定された緯度/経度座標に対応するピクセルを把握できます。
そのプロセスの概要を簡単に説明すると、(-37.803134, 145.132377) 緯度/経度が (0, 0) ピクセルに対応し、(1017, 915) ピクセルが緯度/経度 (- 37.798917、145.138535)。左下隅にピクセル (0, 0) を配置する通常の規則を仮定すると、これは北が画像の上にあることを意味します。
次に、ターゲット座標 (-37.801465, 145.134984) に関心がある場合は、次のように画像の対応するピクセル数を計算できます。
pixelY = ((targetLat - minLat) / (maxLat - minLat)) * (maxYPixel - minYPixel)
= ((-37.801465 - -37.803134) / (-37.798917 - -37.803134)) * (915 - 0)
= 362.138
つまり、対応するピクセルは画像の下から 362 ピクセルです。その後、水平方向のピクセル配置についてもまったく同じことができますが、代わりに経度と X ピクセルを使用します。
このパーツ((targetLat - minLat) / (maxLat - minLat))は、2 つの参照座標間の距離を計算し、最初の座標にいることを示す 0、2 番目の座標にいることを示す 1、およびその間の位置を示す間の数字を示します。したがって、たとえば、2 つの参照座標間で 25% 北にいることを示すには 0.25 が生成されます。最後のビットは、それを同等のピクセルに変換します。
チッ！
編集OK、あなたのコメントに基づいて、もう少し具体的にすることができます。左上隅を主な基準点として使用しているように見えるので、次の定義を使用します。
minLat = -37.803134
maxLat = -37.806232
MAP_HEIGHT = 916
次に、(-37.804465, 145.134984) の座標の例を使用すると、対応するピクセルの Y 座標は、左上隅を基準にして次のようになります。
pixelY = ((targetLat - minLat) / (maxLat - minLat)) * (MAP_HEIGHT - 1)
= ((-37.804465 - -37.803134) / (-37.806232 - -37.803134)) * 915
= 393.11
したがって、対応するピクセルは上から 393 ピクセル下になります。水平方向に相当するものを自分で考えてみましょう。基本的には同じです。注:-1が付いているのMAP_HEIGHTは、ゼロから開始すると、最大ピクセル数が 916 ではなく 915 になるためです。
編集:この機会に指摘したいことの 1 つは、これは概算であるということです。実際には、緯度と経度の座標と他の形式のデカルト座標の間に単純な線形関係はありません。これには、地図を作成する際に使用される投影法や、地球が完全な球体ではないという事実など、さまざまな理由があります。小さな領域では、この近似値は十分に近く、大きな違いはありませんが、大きなスケールでは不一致が明らかになる場合があります。あなたのニーズに応じて、YMMV. （以下の回答で、これが事実であることを思い出させてくれたurayに感謝します）。
測地 (lot,lan) を定義したデカルト座標 (基準点からの x,y メートル) に正確に変換したい場合は、ここのコード スニペットで行うことができます。この関数は測地座標をラジアンで受け入れ、結果は x,y
入力：
コード：
#define GD_semiMajorAxis 6378137.000000
#define GD_TranMercB 6356752.314245
#define GD_geocentF 0.003352810664
void geodeticOffsetInv( double refLat, double refLon,
double lat, double lon,
double& xOffset, double& yOffset )
{
double a = GD_semiMajorAxis;
double b = GD_TranMercB;
double f = GD_geocentF;
double L = lon-refLon
double U1 = atan((1-f) * tan(refLat));
double U2 = atan((1-f) * tan(lat));
double sinU1 = sin(U1);
double cosU1 = cos(U1);
double sinU2 = sin(U2);
double cosU2 = cos(U2);
double lambda = L;
double lambdaP;
double sinSigma;
double sigma;
double cosSigma;
double cosSqAlpha;
double cos2SigmaM;
double sinLambda;
double cosLambda;
double sinAlpha;
int iterLimit = 100;
do {
sinLambda = sin(lambda);
cosLambda = cos(lambda);
sinSigma = sqrt((cosU2*sinLambda) * (cosU2*sinLambda) +
(cosU1*sinU2-sinU1*cosU2*cosLambda) *
(cosU1*sinU2-sinU1*cosU2*cosLambda) );
if (sinSigma==0)
{
xOffset = 0.0;
yOffset = 0.0;
return ; // co-incident points
}
cosSigma = sinU1*sinU2 + cosU1*cosU2*cosLambda;
sigma = atan2(sinSigma, cosSigma);
sinAlpha = cosU1 * cosU2 * sinLambda / sinSigma;
cosSqAlpha = 1 - sinAlpha*sinAlpha;
cos2SigmaM = cosSigma - 2*sinU1*sinU2/cosSqAlpha;
if (cos2SigmaM != cos2SigmaM) //isNaN
{
cos2SigmaM = 0; // equatorial line: cosSqAlpha=0 (§6)
}
double C = f/16*cosSqAlpha*(4+f*(4-3*cosSqAlpha));
lambdaP = lambda;
lambda = L + (1-C) * f * sinAlpha *
(sigma + C*sinSigma*(cos2SigmaM+C*cosSigma*(-1+2*cos2SigmaM*cos2SigmaM)));
} while (fabs(lambda-lambdaP) > 1e-12 && --iterLimit>0);
if (iterLimit==0)
{
xOffset = 0.0;
yOffset = 0.0;
return; // formula failed to converge
}
double uSq = cosSqAlpha * (a*a - b*b) / (b*b);
double A = 1 + uSq/16384*(4096+uSq*(-768+uSq*(320-175*uSq)));
double B = uSq/1024 * (256+uSq*(-128+uSq*(74-47*uSq)));
double deltaSigma = B*sinSigma*(cos2SigmaM+B/4*(cosSigma*(-1+2*cos2SigmaM*cos2SigmaM)-
B/6*cos2SigmaM*(-3+4*sinSigma*sinSigma)*(-3+4*cos2SigmaM*cos2SigmaM)));
double s = b*A*(sigma-deltaSigma);
double bearing = atan2(cosU2*sinLambda, cosU1*sinU2-sinU1*cosU2*cosLambda);
xOffset = sin(bearing)*s;
yOffset = cos(bearing)*s;
}
地球の曲率についてはあまり気にしません。私はこれまでopenstreetmapを使用したことがありませんが、ざっと見てみると、メルカトル図法を使用しているようです。
これは単に、惑星を平らにして長方形にし、Xを経度に比例させ、Yを緯度にほぼ正確に比例させることを意味します。
したがって、先に進んでMacの簡単な数式を使用すると、非常に正確になります。あなたの緯度は、あなたが扱っている小さな地図のピクセルの価値よりはるかに少ないでしょう。ビクトリアのサイズの地図でも、2〜3％の誤差しかありません。
diverscuba23は、楕円体を選択する必要があることを指摘しました... openstreetmapはWGS84を使用しており、最新のマッピングも同様です。ただし、オーストラリアの多くの地図では、100〜200メートルほど異なる可能性のある古いAGD66が使用されていることに注意してください。