Geodetic Latitude (surface normal)

CesiumJS
1

I was looking over how Cesium calculates Geodetic Latitude (surface normal.) Well actually latitude is just one aspect this function determines, I don’t see a function that just determines GD latitude.

Ellipsoid.prototype.cartesianToCartographic

var p = this.scaleToGeodeticSurface(cartesian, cartesianToCartographicP);

calls Ellipsoid.prototype.scaleToGeodeticSurface

var n = this.geodeticSurfaceNormal(p, cartesianToCartographicN);

calls Ellipsoid.prototype.geodeticSurfaceNormal

Ellipsoid.prototype.geodeticSurfaceNormal

Is pretty basic: (positionX/radiiX^2 , positionY/radiiY^2, positionZ/radiiZ^2)

However Ellipsoid.prototype.scaleToGeodeticSurface is not so straight forward. It seems to go through a loop starting with an approximation till an epsilon threshold is passed.

Well I was trying out another method to obtain the surface normal at any height from the surface. It’s quite simple really and the results are surprisingly close to what Cesium currently calculates. In fact it’s exactly what Cesium calculates at height 0.

What I do is determine the foci points for the ellipse cross-section for the point being measured.

I create 2 vectors:

-1 going through the far foci point and the camera position

-1 going through the near foci point and the camera position

The average pitch of these 2 vectors is the surface normal pitch for that camera position.

It’s exact at height 0. At height 6 million meters the biggest difference from what Cesium calculates is 0.048 deg at 45 deg latitude. But which is accurate?

I was reading this site

http://www.algebra.com/algebra/homework/Quadratic-relations-and-conic-sections/REVIEW-of-lessons-on-ellipses.lesson

I noticed this quote

“For any ellipse’s point the normal to the ellipse at this point bisects the angle between the straight lines drawn from the ellipse foci to the point.”

If you look at the image by that quote you’ll notice it’s indicating the same thing.

Perhaps this holds true even for points away from the surface? One could test by placing the camera on the surface, determine the normal, extend the camera 6 million meters in the direction the normal, then calculate both ways and see what results.

There may be an easier method to determine height as well, but I have not tested it yet.

2

Hyper,

For details on how Cesium computes surface normals, start on Page 22 here: http://books.google.com/books?id=x3ttSs_bd_4C&lpg=PP1&pg=PA19#v=onepage&q&f=false

For a demo, run Source/Examples/Chapter02/EllipsoidSurfaceNormals from this repo: https://github.com/virtualglobebook/OpenGlobe

Patrick

3

Thanks for the links Patrick, I’ll be sure to check them out.

Here’s a SandCastle app demonstrating the foci method of latitude and height I was mentioning. The values are remarkably close to those given by Cesium.

var viewer = new Cesium.Viewer(‘cesiumContainer’);

var camera=viewer.camera;

function distP2P_2D(pointA,pointB)

{

var delta = {x:pointB.x-pointA.x,y:pointB.y-pointA.y};

return Math.sqrt(Math.pow(delta.x,2) + Math.pow(delta.y,2));

}

viewer.clock.onTick.addEventListener(function(clock)

{

var buf="";

var mj = 6378137;var mi = 6356752.3142451793; //major,minor

var mjs = Math.pow(mj,2);var mis = Math.pow(mi,2);

var fd = Math.sqrt(mjs - mis); //foci distance

var fl={x:-fd,y:0};var fr={x:fd,y:0}; //foci left and right

var horiz = Math.sqrt(Math.pow(camera.position.x,2) + Math.pow(camera.position.y,2));

var p={x:horiz,y:camera.position.z};

//determine height from surface

var surf_rad = 2*mj;

var elip_rad = distP2P_2D(fl,p) + distP2P_2D(fr,p);

var elip_height = (elip_rad-surf_rad)/2;

buf+=" foci height "+elip_height.toFixed(3);

buf+=" Cesium height "+camera.positionCartographic.height.toFixed(3);

var heightdif=Math.abs(camera.positionCartographic.height-elip_height);

buf+=" foci_Cesium dif: "+heightdif.toFixed(3);

buf+=" % dif "+(heightdif/camera.positionCartographic.height*100).toFixed(3);

buf+="
";

//determine latitude

var from_fl = new Cesium.Cartesian2();

Cesium.Cartesian2.subtract(p,fl,from_fl);

var from_fl_ang = Math.atan2(from_fl.y,from_fl.x);

var from_fr = new Cesium.Cartesian2();

Cesium.Cartesian2.subtract(p,fr,from_fr);

var from_fr_ang = Math.atan2(from_fr.y,from_fr.x);

var between=(from_fl_ang+from_fr_ang)/2;

buf+=" foci latitude: "+Cesium.Math.toDegrees(between).toFixed(3);

buf+=" Cesium latitude: "+Cesium.Math.toDegrees(camera.positionCartographic.latitude).toFixed(3);

var latdif=Math.abs(camera.positionCartographic.latitude-between);

buf+=" foci_Cesium dif: "+Cesium.Math.toDegrees(latdif).toFixed(3);

buf+=" % dif "+(latdif/camera.positionCartographic.latitude*100).toFixed(3);

document.getElementById(“toolbar”).innerHTML=buf;

});

``

4

I created a test case and determined that Cesium’s method is accurate. Though with some changes to my foci method I believe there’s alternate method to calculate height and latitude accurately. If I manage to find a way I’ll post it here.

At the un-scaled point both methods yield an accurate 62.022 latitude and 0 height.

At the scaled point

-my current foci method yields 61.988 latitude and 2994666.895 height

-while Cesium yields an accurate 62.022 latitude and an accurate 3000000.000 height

some eval code used to obtained the points

//get a point on the surface

var mj = 6378137;var mi = 6356752.3142451793;

var mjs = Math.pow(mj,2);var mis = Math.pow(mi,2);

x=3000000;

//x = Math.sqrt( mjs - mjs*(Math.pow(y,2)/mis));

y = Math.sqrt( mis - mis*(Math.pow(x,2)/mjs));

console.log(“xy”,x,y);

//normal

var nx = x/mjs;var ny = y/mis

var mag = Math.sqrt(Math.pow(nx,2)+Math.pow(ny,2));

nx = nx/mag;ny = ny/mag;

console.log(“nxny”,nx,ny);

//point scaled along height

var scale = 3000000;

var sx = x + nx * scale;

var sy = y + ny * scale;

console.log(“sxsy”,sx,sy);

``

When transferring these points to the SandCastle app

x to x

y to z

set y to 0

5

I had hypothesized that if you were to extrude each point on the surface of an ellipse a fixed distance in the direction of their respective normals that you’d end up with a larger ellipse. I’ve tested this hypothesis using Earth’s ellipse of longitude and extruding 3 mega meters outward, just sampling 2 points. If you extrude 3,000,000 meters from Earth’s surface the resulting major radius is 6378137 + 3000000 = 9378137, and this squared is 87949453590769. Knowing a point on the ellipse and the major radius, the minor radius can be computed.

Unfortunately the resulting minor radius of each point do not match each-other, even though each point is at the same height from the reference ellipse. Here’s a test Sandcastle app demonstrating this, you can switch points with by swapping 0 and 1 in the if statement. I’m a bit surprised that extruding along heights doesn’t form an ellipse.

var viewer = new Cesium.Viewer(‘cesiumContainer’);

var camera=viewer.camera;

function distP2P_2D(pointA,pointB)

{

var delta = {x:pointB.x-pointA.x,y:pointB.y-pointA.y};

return Math.sqrt(Math.pow(delta.x,2) + Math.pow(delta.y,2));

}

viewer.clock.onTick.addEventListener(function(clock){

var CC3=Cesium.Cartesian3;

CC3.negate(camera.position,camera.direction);

CC3.cross({x:0,y:0,z:1},camera.position,camera.right);

CC3.cross(camera.right,camera.direction,camera.up);

camera.frustum.fov=40/180*Math.PI;

if(0)

{//3000000 height above ref

camera.position.x=4407381.873624045;

camera.position.y=0;

camera.position.z=8259075.2164053;

var mjs = 87949453590769;

var mis = 87548861332798.006201037300231926;

}

else

{//3000000 height above ref

camera.position.x=5877592.495525159;

camera.position.y=0;

camera.position.z=7291100.110070896;

var mjs = 87949453590769;

var mis = 87548898296883.406616257801218887;

}

var buf="";

var fd = Math.sqrt(mjs - mis); //foci distance

var fl={x:-fd,y:0};var fr={x:fd,y:0}; //foci left and right

var horiz = Math.sqrt(Math.pow(camera.position.x,2) + Math.pow(camera.position.y,2));

var p={x:horiz,y:camera.position.z};

//determine height from surface

var surf_rad = 2*6378137;

var elip_rad = distP2P_2D(fl,p) + distP2P_2D(fr,p);

var elip_height = (elip_rad-surf_rad)/2;

buf+=" foci height "+elip_height.toFixed(3);

buf+=" Cesium height "+camera.positionCartographic.height.toFixed(3);

var heightdif=Math.abs(camera.positionCartographic.height-elip_height);

buf+=" foci_Cesium dif: "+heightdif.toFixed(3);

buf+=" % dif "+(heightdif/camera.positionCartographic.height*100).toFixed(3);

buf+="
";

//determine latitude

var from_fl = new Cesium.Cartesian2();

Cesium.Cartesian2.subtract(p,fl,from_fl);

var from_fl_ang = Math.atan2(from_fl.y,from_fl.x);

var from_fr = new Cesium.Cartesian2();

Cesium.Cartesian2.subtract(p,fr,from_fr);

var from_fr_ang = Math.atan2(from_fr.y,from_fr.x);

var between=(from_fl_ang+from_fr_ang)/2;

buf+=" foci latitude: "+Cesium.Math.toDegrees(between).toFixed(3);

buf+=" Cesium latitude: "+Cesium.Math.toDegrees(camera.positionCartographic.latitude).toFixed(3);

var latdif=Math.abs(camera.positionCartographic.latitude-between);

buf+=" foci_Cesium dif: "+Cesium.Math.toDegrees(latdif).toFixed(3);

buf+=" % dif "+(latdif/camera.positionCartographic.latitude*100).toFixed(3);

document.getElementById(“toolbar”).innerHTML=buf;

});

``

Number of observed iterations in the do-while loop of Ellipsoid.prototype.scaleToGeodeticSurface

2-3 low altitude

3 mid altitude

3-4 high altitude

I had thought there would be a simpler and faster method, but as of yet I have not found one.

6

It seems to me that you can combine the equation of an ellipse and the point-slope form of a line to determine the intersect point, and hence the geodetic height and latitude.

knowns:

Cx, Cy = camera position

As = earth major axis squared

Bs = earth minor axis squared

unknowns:

Sx, Sy = underlying surface position

equations:

m = (Sx/As)/(Sy/Bs) //slope of line (normal at surface point)

Cy-Sy = m(Cx-Sx) //equation of line (point slope form)

Sx^2 / As + Sy^2 / Bs = 1 //equation of earth longitude ellipse

objective: solve for either Sx or Sy

rearrangements of the linear equation

(Cy-Sy) = m(Cx-Sx) //point slope form of line

(Cy-Sy) = (Sx/As)/(Sy/Bs)(Cx-Sx) //substitute m

(Cy-Sy)/(Cx-Sx) = (Sx/As)/(Sy/Bs) //div BS by (Cx-Sx)

(Cy-Sy)(Sy/Bs) = (Cx-Sx)(Sx/As) //criss-cross

(CySy/Bs)-(Sy^2/Bs) = (CxSx/As)-(Sx^2/As) //mul nomials

rearrangements of the elliptical equation

(Sx^2/As)+(Sy^2/Bs) = 1 //ellipse equation

(Sx^2/As) = 1 - (Sy^2/Bs) //sub BS by (Sy^2/Bs) //branchA

(Sy^2/Bs) = 1 - (Sx^2/As) //sub BS by (Sx^2/As) //branchB

Sx^2 = As(1 - (Sy^2/Bs)) //mul BS by As //branchA

Sy^2 = Bs(1 - (Sx^2/Bs)) //mul BS by Bs //branchB

Sx = sqrt(As(1 - (Sy^2/Bs))) //mul BS by As //branchA

Sy = sqrt(Bs(1 - (Sx^2/Bs))) //mul BS by Bs //branchB

Then substitute either elliptical branch A or B into the linear equation then either isolate or put into a quadratic form. That sqrt is a hassle though. If you do manage to isolate it, the other side of the equation can become huge when squaring both sides.

If there’s any faulty logic here please let me know.