7-Parameter Molodensky Transformation: Coordinate Conversion Using Bursa-Wolf Model

Math Problem Statement

%The 7-parameter Molodensky transformation involves converting coordinates between two different geodetic systems with slightly different reference ellipsoids using parameters for translation, rotation, and scale change. Here is the code modified to use the 7-parameter Molodensky transformation: % Program to compute Transformation parameters from SL old to SLD99 clc; clear; format long;

% Computed Latitudes and longitudes values from earlier programs are used % with orthometric heights of the points

[PI, PM, PS, LI, LM, LS, HI] = textread('known-points.txt', '%f %f %f %f %f %f %f'); [PRR, LRR, Hv] = textread('Everest.txt', '%f %f %f'); % Determine the number of points NPoints = min(length(PI), length(PRR)); % Constants for WGS84 ag = 6378137.0; bg = 6356752.3142; eg = sqrt(((ag^2) - (bg^2)) / (ag^2));

for i = 1:NPoints % Convert degree, minute, second to decimal degrees PR(i) = PI(i) + (PM(i) / 60) + (PS(i) / 3600); LR(i) = LI(i) + (LM(i) / 60) + (LS(i) / 3600);

ng(i) = (eg * sind(PR(i)))^2;
vg(i) = ag / sqrt(1 - ng(i));
XI(i) = (vg(i) + HI(i)) * cosd(PR(i)) * cosd(LR(i));
YI(i) = (vg(i) + HI(i)) * cosd(PR(i)) * sind(LR(i));
ZI(i) = (vg(i) * (1 - (eg^2)) + HI(i)) * sind(PR(i));

end

% Conversion of (Everest) Lat, lon, and ellipsoidal heights to Geocentric (X, Y, Z) coordinates av = 6377276.345; bv = 6356075.413; ev = sqrt(((av^2) - (bv^2)) / (av^2));

for i = 1:NPoints Nv = av / sqrt(1 - (ev * sind(PRR(i)))^2); Xv(i) = (Nv + Hv(i)) * cosd(PRR(i)) * cosd(LRR(i)); Yv(i) = (Nv + Hv(i)) * cosd(PRR(i)) * sind(LRR(i)); Zv(i) = (Nv * (1 - (ev^2)) + Hv(i)) * sind(PRR(i)); end

C = zeros(NPoints * 3, 7); for i = 1:NPoints C(i * 3 - 2, 1) = XI(i); C(i * 3 - 2, 3) = -ZI(i); C(i * 3 - 2, 4) = YI(i); C(i * 3 - 2, 5) = 1; C(i * 3 - 1, 1) = YI(i); C(i * 3 - 1, 2) = ZI(i); C(i * 3 - 1, 4) = -XI(i); C(i * 3 - 1, 6) = 1; C(i * 3, 1) = ZI(i); C(i * 3, 2) = -YI(i); C(i * 3, 3) = XI(i); C(i * 3, 7) = 1; end

y = zeros(NPoints * 3, 1); for i = 1:NPoints y(i * 3 - 2, 1) = Xv(i) - XI(i); y(i * 3 - 1, 1) = Yv(i) - YI(i); y(i * 3, 1) = Zv(i) - ZI(i); end

% Bursa wolf 7 parameters are computed using least square principal. S = (inv(C' * C)) * (C' * y); Rx = ((180 / pi) * 3600) * S(2); Ry = ((180 / pi) * 3600) * S(3); Rz = ((180 / pi) * 3600) * S(4); Dx = S(5); Dy = S(6); Dz = S(7); Sf = 1000000 * S(1);

disp(S); M = [Dx; Dy; Dz; Rx; Ry; Rz; Sf]; fprintf('%.6f \n%.6f \n%.6f \n%.6f \n%.6f \n%.6f \n%.6f\n', Dx, Dy, Dz, Rx, Ry, Rz, Sf);

Solution

The given MATLAB code implements the 7-parameter Molodensky-Badekas transformation to convert geodetic coordinates (latitude, longitude, and height) from one reference ellipsoid to another. The transformation uses the Bursa-Wolf model, which involves translation, rotation, and scaling parameters.

Code Explanation:

Data Import:
- The code reads known latitude, longitude, and height values (PI, PM, PS, LI, LM, LS, HI) from known-points.txt.
- It also reads the Everest ellipsoid values (PRR, LRR, Hv) from Everest.txt.
WGS84 Ellipsoid Constants:
- ag and bg are the semi-major and semi-minor axes of the WGS84 ellipsoid.
- eg is the first eccentricity of the WGS84 ellipsoid, calculated from these axes.
Conversion to Cartesian (X, Y, Z):
- The geodetic latitude and longitude in degrees, minutes, and seconds (DMS) are converted into decimal degrees.
- Then, the geodetic coordinates (lat, lon, height) are transformed into Earth-Centered Earth-Fixed (ECEF) Cartesian coordinates (X, Y, Z) using the ellipsoid parameters of WGS84 and Everest.
Matrix Setup for Least-Squares Solution:
- The transformation matrix C and the residual vector y are constructed. Each point contributes 3 equations corresponding to differences in the X, Y, Z coordinates between the two reference systems.
7-Parameter Solution (Bursa-Wolf Parameters):
- The least-squares method is used to compute the transformation parameters (S):
  - S(1): Scale factor.
  - S(2), S(3), S(4): Rotation parameters (in radians, converted to arc-seconds).
  - S(5), S(6), S(7): Translation parameters (in meters).
- The rotations are converted from radians to arc-seconds, and the scale factor is adjusted by multiplying by $10^6$ .
Output:
- The 7 parameters are displayed:
  - Translations $D_x$ , $D_y$ , $D_z$ (in meters).
  - Rotations $R_x$ , $R_y$ , $R_z$ (in arc-seconds).
  - Scale factor $S_f$ .

Summary of Output Variables:

$Dx, Dy, Dz$ : Translation parameters (meters).
$Rx, Ry, Rz$ : Rotation parameters (arc-seconds).
$Sf$ : Scale factor (ppm).

Example Output:

After running the script with proper input data, the output would display the 7 transformation parameters in the following format:
```
Dx
Dy
Dz
Rx
Ry
Rz
Sf
```

Questions:

What are the specific differences between the WGS84 and Everest ellipsoids?
How does the Molodensky transformation differ from other geodetic transformations?
Can this script be adapted for other ellipsoid pairs, such as GRS80?
What are the typical applications of the 7-parameter transformation in geodesy?
How does the least-squares method ensure the accuracy of the computed parameters?

Tip:

Ensure the input geodetic data files (known-points.txt and Everest.txt) are formatted correctly and contain consistent units (degrees, meters).

Ask a new question for Free

By Image

Drop file here or Click Here to upload

Math Problem Analysis

Mathematical Concepts

Geodesy
Coordinate Transformation
Least-Squares Method

Formulas

First eccentricity formula: e = sqrt((a^2 - b^2) / a^2)
ECEF (Earth-Centered, Earth-Fixed) Cartesian coordinate conversion formula
Least-squares solution formula: S = (C' * C)^(-1) * C' * y
Rotation to arc-seconds: R = (180 / pi) * 3600 * S

Theorems

Bursa-Wolf 7-parameter transformation
Molodensky-Badekas transformation

Suitable Grade Level

Undergraduate Geodesy or Geomatics Engineering

Related Recommendation

MATLAB Geodetic Coordinate Transformation using Bursa-Wolf Model

Coordinate Transformation from WGS84 to SLD99 Using Latitude and Longitude in DMS Format

MATLAB Code for Veis Transformation Parameters and Standard Deviations

Compute (U,V,W) Using Clarke 1880 Ellipsoidal Parameters and Find Foot Point

Coordinate Transformation from WGS84 to SLD99 System Using MATLAB