Chuyển đổi độ dài/vĩ độ thành các tọa độ X/Y

Question

Tôi có giá trị Lat/Long của một khu vực nhỏ ở Melbourne; -37.803134,145.132377 và cũng là một hình ảnh phẳng mà tôi xuất từ ​​bản đồ openstreet (Osmarender Image). Chiều rộng hình ảnh: 1018 và Chiều cao: 916

Tôi muốn chuyển đổi, sử dụng C++, độ dài Lat/Long thành X, Y, nơi điểm sẽ phản ánh vị trí.

Tôi đã sử dụng nhiều công thức khác nhau mà tôi tìm thấy trên web như được đưa ra dưới đây nhưng không có gì hữu ích.

var y = ((-1 * lat) + 90) * (MAP_HEIGHT/180); 
var x = (lon + 180) * (MAP_WIDTH/360); 

Sẽ rất hữu ích nếu có ai giải thích rõ ràng về cách thực hiện việc này. Bất kỳ mã nào cũng sẽ được đánh giá cao.

Answer 1

Bạn cần thêm thông tin hơn chỉ là một cặp lat/lon đơn để có thể thực hiện việc này.

Ở giai đoạn này, các thông tin mà bạn đã cung cấp là mất tích hai điều:

• cách lớn diện tích làm bìa hình ảnh của bạn (về lat/lon)? Dựa trên những gì bạn đã cung cấp, tôi không biết liệu hình ảnh có hiển thị diện tích rộng một mét hay rộng một cây số hay không.
• tọa độ trên hình ảnh của bạn trong tọa độ tham chiếu của bạn (-37.803134, 145.132377) là gì? Nó là một trong những góc? Ở giữa một nơi nào đó?

Tôi cũng sẽ giả định rằng hình ảnh của bạn được căn chỉnh về phía bắc/nam - ví dụ: nó không có hướng bắc về phía góc trên cùng bên trái. Điều đó sẽ có xu hướng làm phức tạp mọi thứ.

Cách tiếp cận dễ nhất là tìm ra chính xác tọa độ vĩ độ/vĩ độ tương ứng với điểm ảnh (0, 0) và pixel (1017, 915). Sau đó, bạn có thể tìm ra pixel tương ứng với toạ độ vĩ độ/kinh độ đã cho qua interpolation.

Để phác thảo ngắn gọn quá trình đó, hãy tưởng tượng rằng vĩ độ (/ 37,803134, 145.132377) của bạn tương ứng với pixel (0, 0) và bạn đã phát hiện ra điểm ảnh (1017, 915) tương ứng với lat/lon (-37,798917, 145.138535). Giả sử quy ước thông thường với pixel (0, 0) ở góc dưới cùng bên trái, điều này có nghĩa là phía bắc nằm trong ảnh.

Sau đó, nếu bạn quan tâm đến mục tiêu phối hợp (-37,801465, 145,134984), bạn có thể tìm ra các số tương ứng của điểm ảnh lên hình ảnh nó là như sau:

pixelY = ((targetLat - minLat)/(maxLat - minLat)) * (maxYPixel - minYPixel) 
     = ((-37.801465 - -37.803134)/(-37.798917 - -37.803134)) * (915 - 0) 
     = 362.138 

Đó là, tương ứng pixel là 362 pixel từ dưới cùng của hình ảnh. Sau đó, bạn có thể thực hiện tương tự cho vị trí pixel ngang, nhưng thay vào đó, sử dụng kinh độ và pixel X.

Phần ((targetLat - minLat)/(maxLat - minLat)) tính toán khoảng cách giữa hai toạ độ tham chiếu và cho 0 để cho biết bạn đang ở vị trí đầu tiên, 1 để chỉ số thứ hai và số ở giữa để cho biết vị trí ở giữa. Vì vậy, ví dụ, nó sẽ tạo ra 0,25 để chỉ ra rằng bạn là 25% về phía bắc giữa hai tọa độ tham chiếu. Bit cuối cùng chuyển đổi nó thành các pixel tương đương.

HTH!

EDIT Ok, dựa trên nhận xét của bạn, tôi có thể cụ thể hơn một chút.Cho rằng bạn đang sử dụng góc trên bên trái là điểm tham chiếu chính của bạn, tôi sẽ sử dụng các định nghĩa sau đây:

minLat = -37.803134 
maxLat = -37.806232 
MAP_HEIGHT = 916 

Sau đó, nếu chúng tôi sử dụng một ví dụ toạ độ (-37,804465, 145,134984), Y tọa độ của pixel tương ứng, so với góc trên cùng bên trái, là:

pixelY = ((targetLat - minLat)/(maxLat - minLat)) * (MAP_HEIGHT - 1) 
     = ((-37.804465 - -37.803134)/(-37.806232 - -37.803134)) * 915 
     = 393.11 

Do đó, điểm ảnh tương ứng là 393 pixel từ trên xuống. Tôi sẽ cho phép bạn làm việc ra tương đương ngang cho chính mình - về cơ bản là giống nhau. LƯU Ý Các -1 với MAP_HEIGHT là bởi vì nếu bạn bắt đầu từ số không, số điểm ảnh tối đa là 915, không 916.

EDIT: Một điều tôi muốn tận dụng cơ hội để chỉ ra là điều này là một xấp xỉ. Trong thực tế, không có mối quan hệ tuyến tính đơn giản giữa vĩ độ và kinh độ và các dạng khác của tọa độ Descartes vì ​​một số lý do bao gồm các dự báo được sử dụng trong khi tạo bản đồ và thực tế là Trái Đất không phải là một hình cầu hoàn hảo. Trong các khu vực nhỏ, xấp xỉ này đủ gần để không tạo ra sự khác biệt đáng kể, nhưng về những sự khác biệt về quy mô lớn hơn có thể trở nên rõ ràng. Tùy thuộc vào nhu cầu của bạn, YMMV. (Tôi cảm ơn uray, câu trả lời dưới đây nhắc tôi rằng đây là trường hợp).

Answer 2

nếu bạn đang tìm kiếm chuyển đổi chính xác trắc địa (lô, lan) thành tọa độ được xác định của bạn (x, y mét từ điểm tham chiếu), bạn có thể thực hiện với đoạn mã của tôi tại đây, chức năng này sẽ chấp nhận tọa độ trắc địa radian và đầu ra kết quả trong x, y

đầu vào:

• refLat, refLon: trắc địa phối hợp mà bạn định nghĩa là 0,0 trong Descartes phối hợp (đơn vị là trong radian)
• lat , lon: geodetic.210 phối hợp mà bạn muốn tính toán Descartes của mình phối hợp (đơn vị là trong radian)
• xoffset, yOffset: kết quả trong Descartes tọa độ x, y (đơn vị là tính bằng mét)

mã:

#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; 
} 

Answer 3

Tôi sẽ không lo lắng quá nhiều về độ cong của trái đất. Tôi đã không sử dụng openstreetmap trước đây, nhưng tôi đã có một cái nhìn nhanh chóng và nó xuất hiện rằng họ sử dụng một chiếu Mercator.

Điều đơn giản có nghĩa là chúng đã làm phẳng hành tinh cho bạn lên một hình chữ nhật, làm cho X tỷ lệ thuận với Longitude và Y gần như chính xác tỷ lệ với Latitude.

Vì vậy, bạn có thể tiếp tục và sử dụng các công thức đơn giản của Mac và bạn sẽ rất chính xác. Latitude của bạn sẽ được giảm giá ít hơn nhiều so với giá trị của pixel đối với bản đồ nhỏ mà bạn đang xử lý. Ngay cả trên bản đồ, kích thước của Victoria bạn sẽ chỉ nhận được một lỗi 2-3%.

diverscuba23 chỉ ra rằng bạn phải chọn một ellipsoid ... openstreetmap sử dụng WGS84 và bản đồ hiện đại nhất cũng vậy. Tuy nhiên, hãy cẩn thận rằng nhiều bản đồ ở Úc sử dụng AGD66 cũ hơn có thể khác nhau từ 100-200 mét.?

Answer 4

double testClass::getX(double lon, int width) 
{ 
    // width is map width 
    double x = fmod((width*(180+lon)/360), (width +(width/2))); 

    return x; 
} 

double testClass::getY(double lat, int height, int width) 
{ 
    // height and width are map height and width 
    double PI = 3.14159265359; 
    double latRad = lat*PI/180; 

    // get y value 
    double mercN = log(tan((PI/4)+(latRad/2))); 
    double y  = (height/2)-(width*mercN/(2*PI)); 
    return y; 
}