WaBis

walter.bislins.ch

Globe Arc Projection to Flat Earth

Friday, October 19, 2018 - 01:51 | Author: wabis | Topics: FlatEarth, Mathematics, Geometry, Interactive
This App shows how a circle given by its center in latitude φ and longitude λ and its radius s along the surface of the globe earth looks like, if its curve is projected from the globe to the flat earth. I show the formulas used, derived by myself, and the source code of the working App.

The App

Problem

There is a problem if the circle crosses the south pole, because on the globe the south pole is a singe point while on the flat earth the south pole is the whole circumference of the earth. A single point can not be projected to a circle, so the corresponding projection is undefined. If the circle crosses near the south pole the resolution of the projection gets problematic, because for the south pole region there are only a few points of the circle left that have to be distributed around half the flat earth.

Formulas used

It is not practical to give a single formula for the projection of a circle point on the globe into a corresponding circle point on the flat earth, not to speak of a formula that transforms a whole circle equation. The formula would be simply too long and unstructured, if possible at all. So the calculation is divided into meaningfull parts.

Given are the center of the circle in latitude φ and longitude λ and the readius s measured along the surface of the globe. A single point on the circle is defined by the angle α.

The following principle is applied to derive the formulas: A point of a circle around the north pole is computed. This point is then moved along the surface of the globe to the location, defined by the lat/long of the center of the circle, by rotating the globe in two steps accordingly. This results in a point on the globe given in X/Y/Z coordinates. These coordinates are then converted back into lat/long of the point. This lat/long of the point are the same on the flat earth. Finally the lat/long of the flat earth are converted into X/Y coordinates of the flat earth image, so it can be drawn with a graphics programm.

Computing a Circle Point on the Globe

First we need some intermediate values to simplify the equations:

(1)
(2)
with
where'
' =' 'radius of the earth
' =' 'angle that specifies a specific point of the circle
' =' 'Latitude of the circle center
' =' 'Longitude of the circle center
' =' 'radius of the circle along the surface of the globe
' =' 'function that returns 1 for a positive argument and -1 for a negative argument

The intermediate variable c,u,v and are now used in the following formulas:

(3)
where'
' =' 'coordinates of the arc point in the globe centered corrdinate system
' =' 'Latitude of the circle center
' =' 'Longitude of the circle center

Computing Lat/Long of the Circle Point

(4)
with
where'
' =' 'Latitude of the arc point
' =' 'Longitude of the arc point
' =' 'function that returns 1 for a positive argument and -1 for a negative argument

Latitude and Longitude of the arc point are the same for the globe and the flat earth.

If q = 0, which is the case if the point (x,y,z) is at one of the poles, the longitude λ is undefined and set to 0 in the App.

Converting Lat/Long into Cartesian Coordinates

The cartesian corrdinate system in the App is such that points down and points to the right.

(5)
with
where'
' =' 'cartesian coordinates of the arc point on the flat earth
' =' 'Latitude of the arc point
' =' 'Longitude of the arc point
' =' 'radius of the earth

Code


#INCLUDE JsGraph.inc
#INCLUDE ControlPanel.inc
#INCLUDE EarthMap.inc

var pi = Math.PI;
var sin = Math.sin;
var cos = Math.cos;
var acos = Math.acos;
function rad(x) { return x * pi / 180; }
function deg(x) { return x * 180 / pi; }
function len(x,y) { return Math.sqrt( x*x + y*y ) }
function sign(x) { return x >= 0 ? 1 : -1; }

var Model = {
  // input
  ArcCenterLat: 0,
  ArcCenterLong: 0,
  ArcSurfRadius: 2500000,
  ArcStartAngle: 0,
  ArcEndAngle: 360,
  ArcDeltaAngle: 0.1,

  // output as f( ArcCenter, ArcSurfRadius, ArcAngle )
  Xglobe: 0,
  Yglobe: 0,
  Zglobe: 0,
  LatArcPoint: 0,
  LongArcPoint: 0,
  Xfe: 0,
  Yfe: 0,

  // constants
  R: 6371000,
  Re: Math.PI/2 * 6371000,

  PointGlobeToFe: function( angle, radius ) {

    // comp globe coordiantes X,Y,Z of arc point
    var alpha = rad( angle + this.ArcCenterLong );
    var gamma = radius / this.R;
    var rPrime = this.R * sin( gamma );
    var a = rPrime * sin( alpha );
    var b = this.R * cos( gamma );
    var c = len( a, b );
    var beta = sign(a) * acos( b / c );
    var latPrime = rad( 90 - this.ArcCenterLat );
    var u = rPrime * cos( alpha );
    var v = c * sin( beta + latPrime );
    var lambda = rad( this.ArcCenterLong );
    this.Xglobe = u * cos( lambda ) + v * sin( lambda );
    this.Yglobe = -u * sin( lambda ) + v * cos( lambda );
    this.Zglobe = c * cos( beta + latPrime );

    // comp lat/long of arc point on globe = lat/long of arc on FE
    this.LatArcPoint = 90 - deg( acos( this.Zglobe / this.R ) );
    var rPrimeArcPoint = len( this.Xglobe, this.Yglobe );
    if (rPrimeArcPoint == 0) {
      // LongArcPoint is undefined at the poles
      this.LongArcPoint = 0;
    } else {
      this.LongArcPoint = deg( sign( this.Yglobe ) * acos( this.Xglobe / rPrimeArcPoint ) );
    }

    // comp arc point lat/long to FE coordinates X/Y
    var xy = this.LatLongToFeXY( this.LatArcPoint, this.LongArcPoint );
    this.Xfe = xy[0];
    this.Yfe = xy[1];
  },

  LatLongToFeXY: function( lat, long ) {
    var d = this.Re * ( 1 - lat / 90 );
    var x = d * cos( rad( long ) );
    var y = d * sin( rad( long ) );
    return [ x, -y ];
  },

  DrawFE: function( g ) {

    EarthMap.Radius = 2 * this.Re;
    EarthMap.SetWaterColor( '#d3e2f5' );
    EarthMap.SetLakeColor( '#d3e2f5', '#8cbe5d' );
    EarthMap.SetContinentColor( null, '#c6dfaf', '#8cbe5d' );
    EarthMap.SetLandColor( 'Antarctica', '#eee', '#ccc' );
    EarthMap.DrawFlatEarth( g );

    g.SetAlpha( 0.3 );
    g.SetLineAttr( 'gray', 1 );
    EarthMap.DrawFlatEarthGrid( g, 15, 15 );
    g.SetLineAttr( 'black', 1 );
    EarthMap.DrawFlatEarthEquator( g );
    EarthMap.DrawFlatEarthBorder( g );
    EarthMap.DrawFlatEarthMeridian( g );
    g.SetAlpha( 1 );

  },

  Draw: function( g ) {

    g.Reset();
    g.SetAngleMeasure( 'deg' );
    g.MapWindow( 0, 0, 4.05*this.Re, 4.05*this.Re );

    // draw flat earth 
    this.DrawFE( g );

    // draw circle center
    g.SetMarkerAttr( 'Plus', 12, 'red', 'red', 2 );
    this.PointGlobeToFe( 0, 0 )
    g.Marker( this.Xfe, this.Yfe );

    g.SetLineAttr( 'blue', 2 );
    var startAngle = this.ArcStartAngle;
    var endAngle = this.ArcEndAngle;
    if (startAngle > endAngle) {
      var tmp = startAngle;
      startAngle = endAngle;
      endAngle = tmp;
    }
    if (endAngle - startAngle > 360) endAngle = startAngle + 360;
    g.NewPoly();
    for (var angle = startAngle; angle <= endAngle; angle += this.ArcDeltaAngle) {
      this.PointGlobeToFe( angle, this.ArcSurfRadius );
      g.AddPointToPoly( this.Xfe, this.Yfe );
    }
    g.DrawPoly( 1 );
  },
};

var graph = NewGraph2D( {
  Width: '100%',
  Height: '66%',
  DrawFunc: function(g){ Model.Draw( g ); },
  AutoReset: true,
  AutoClear: true,
} );

function UpdateAll() {
  ControlPanels.Update();
  Model.Draw( graph );
}


ControlPanels.NewSliderPanel( { 
  Name: 'SliderPanel1',
  ModelRef: 'Model',
  NCols: 1, 
  ValuePos: 'left',
  OnModelChange: UpdateAll, 
  Format: 'std', 
  Digits: 5,
  ReadOnly: false, 
  PanelFormat: 'InputNormalWidth'

} ).AddValueSliderField( {
  Name: 'ArcCenterLat',
  Color: 'blue',
  Min: -90,
  Max: 90,
  Units: '°',

} ).AddValueSliderField( {
  Name: 'ArcCenterLong',
  Color: 'blue',
  Min: -180,
  Max: 180,
  Units: '°',

} ).AddValueSliderField( {
  Name: 'ArcSurfRadius',
  Color: 'red',
  Min: 0,
  Max: 2*Model.Re,
  Units: 'km',
  Mult: 1000,

} ).AddValueSliderField( {
  Name: 'ArcStartAngle',
  Color: 'green',
  Min: -360,
  Max: 360,
  Units: '°',

} ).AddValueSliderField( {
  Name: 'ArcEndAngle',
  Color: 'green',
  Min: -360,
  Max: 360,
  Units: '°',

} ).Render();

More Page Infos / Sitemap
Created Friday, October 19, 2018
Scroll to Top of Page
Changed Thursday, August 22, 2019