EllipsoidLambertConformalConicProjection

Sep 17, 2009 at 7:18 PM

Hi,

 The current projection classes are using Spheres without considering the flattening. I've created a class EllipsoidLambertConformalConicProjection for my own use that deals with flattening because we are using state plane coordinates system (roughly 60-70% of US State Plane CS is in LCC projection) . I am posting here and hopefully someone may find it useful, even better check into SVN so everyone can benefit. 

new class: EllipsoidLambertConformalConicProjection.cs

//------------------------------------------------------------------------------
// References: //http://pubs.er.usgs.gov/djvu/PP/PP_1395.pdf
// http://www.posc.org/Epicentre.2_2/DataModel/ExamplesofUsage/eu_cs34.html
// for NCSP83: LatLng(-80.5666 35.102363)<  === > XY (1531463.95, 495879.744);

//------------------------------------------------------------------------------
using System;
using System.Collections.Generic;
using System.Diagnostics;
using System.Globalization;

namespace SQLSpatialTools {
  internal sealed class EllipsoidLambertConformalConicProjection : Projection {
    // Angles are in degrees.
    // ivf: inverse flattening 1/f=a/(a-b)
    // latitude0: latitude for the origin of cartesian coordinates
    // fi1: standard parallel
    // fi2: standard parallel
    // semimajor:  major axis in meters. 
    // unit: per meter.
    // fe: false easting
    // fn: false northing
    // longitude0: central meridian longitude
    public EllipsoidLambertConformalConicProjection(Dictionary<String, double> parameters)
      : base(parameters) {
      double phi1 = InputLatitude("fi1");
      double phi2 = InputLatitude("fi2");
      double phiF = InputLatitude("latitude0");
      this._lamdaF = this.CentralLongitudeRad;// longitude0
      this._a = parameters["semimajor"] / parameters["unit"];
      double f_i = parameters["ivf"];
      this._FE = parameters["fe"];
      this._FN = parameters["fn"];
      if (Math.Abs(phi1 - phi2) <= MathX.Tolerance) {
        throw new ArgumentException(Resource.Fi1AndFi2MustBeDifferent);
      }
      double f = 1.0 / f_i; //e: eccentricity of the ellipsoid where e^2  =  2f - f^2 
      double es = 2 * f - f * f;
      this._e = Math.Sqrt(es);
      double m1 = this.Calc_m(phi1, es);
      double m2 = this.Calc_m(phi2, es);
      double tF = this.Calc_t(phiF, this._e);
      double t1 = this.Calc_t(phi1, this._e);
      double t2 = this.Calc_t(phi2, this._e);
      this._n = Math.Log(m1 / m2) / Math.Log(t1 / t2);
      this._F = m1 / (this._n * Math.Pow(t1, this._n));
      this._rF = this.Calc_r(this._a, this._F, tF, this._n);
    }
    private double Calc_m(double phi, double es) {
      double sinphi = Math.Sin(phi);
      return Math.Cos(phi) / Math.Sqrt(1 - es * sinphi * sinphi);
    }
    private double Calc_t(double phi, double e) {
      double esinphi = e * Math.Sin(phi);
      return Math.Tan(Math.PI / 4 - phi / 2) / Math.Pow((1 - esinphi) / (1 + esinphi), e / 2);
    }
    private double Calc_r(double a, double F, double t, double n) {
      return a * F * Math.Pow(t, n);
    }
    private double Calc_phi(double t_i, double e, double phi) {
      var esinphi = e * Math.Sin(phi);
      return Math.PI / 2 - 2 * Math.Atan(t_i * Math.Pow((1 - esinphi) / (1 + esinphi), e / 2));
    }
    private double Solve_phi(double t_i, double e, double init) {
      // iteration
      int i = 0;
      double phi = init;
      double newphi = this.Calc_phi(t_i, e, phi);//this.
      while (Math.Abs(newphi - phi) > 0.000000001 && i < 10) {
        i++;
        phi = newphi;
        newphi = this.Calc_phi(t_i, e, phi);
      }
      return newphi;
    }

    protected internal override void Project(double latitude, double longitude, out double x, out double y) {
      // note latitude and longitude in radius
      double t = this.Calc_t(latitude, this._e);
      double r = this.Calc_r(this._a, this._F, t, this._n);
      double theta = this._n * (longitude);// already normalized in SqlProjection.ProjectPoint - this._lamdaF);
      x = this._FE + r * Math.Sin(theta);
      y = this._FN + this._rF - r * Math.Cos(theta);
    }

    protected internal override void Unproject(double x, double y, out double latitude, out double longitude) {
      double theta_i = Math.Atan((x - this._FE) / (this._rF - (y - this._FN)));
      double r_i = (this._n > 0 ? 1 : -1) * Math.Sqrt((x - this._FE) * (x - this._FE) + (this._rF - (y - this._FN)) * (this._rF - (y - this._FN)));
      double t_i = Math.Pow((r_i / (this._a * this._F)), 1 / this._n);
      latitude = this.Solve_phi(t_i, this._e, 0);
      longitude = theta_i / this._n;// already normalized in SqlProjection.ProjectPoint +this._lamdaF;
    }

    private readonly double _n;
    private readonly double _rF;
    private readonly double _F;
    private readonly double _e;

    private readonly double _a; // major axis
    private readonly double _lamdaF; // central_meridian
    private readonly double _FE; // fasle easting in unit
    private readonly double _FN; // false northing

  }
}

Modify SqlProjection.cs, adding:

    [SqlMethod(IsDeterministic = true, IsPrecise = false)]
    public static SqlProjection EllipsoidLambertConformalConicProjection(double longitude0, double latitude, double fi1, double fi2,
      double semimajor, double ivf, double fe, double fn, double unit) {
      Dictionary<String, double> parameters = new Dictionary<string, double>();
      parameters["longitude0"] = longitude0;
      parameters["latitude0"] = latitude;
      parameters["fi1"] = fi1;
      parameters["fi2"] = fi2;
      parameters["semimajor"] = semimajor;
      parameters["ivf"] = ivf;
      parameters["fe"] = fe;
      parameters["fn"] = fn;
      parameters["unit"] = unit;
      return new SqlProjection(new EllipsoidLambertConformalConicProjection(parameters));
    }

After compile and unregister / register, run sample with North Carolina State Plane NAD83 feet:

declare @ncsp Projection
set @ncsp = Projection::EllipsoidLambertConformalConicProjection(-79.0, 33.75, 34.33333333333334, 36.16666666666666, 6378137.0, 298.257222101, 2000000.002616666, 0, 0.3048006096012192)
select @ncsp.Project('POINT (-80.5666 35.102363)').ToString()
-- >> POINT (1531463.9435256813 495879.7353486158)
select @ncsp.UnprojectWithSRID('POINT (1531463.95 495879.744)', 4326).ToString()
-- >> POINT (-80.566599978812476 35.102363023975727)
select @ncsp.Unproject(@ncsp.Project('LINESTRING (-80 35, -81 36)')).ToString()
-- >> LINESTRING (-80 34.999999999925748, -81 35.999999999931944)
select @ncsp.ToString()
-- >> SQLSpatialTools.EllipsoidLambertConformalConicProjection longitude0=-79;latitude0=33.75;fi1=34.333333333333343;fi2=36.166666666666657;semimajor=6378137;ivf=298.257222101;fe=2000000.002616666;fn=0;unit=0.30480060960121919

 

Thanks.