Michigan Georef projection support?

Previous Topic Next Topic
 
classic Classic list List threaded Threaded
2 messages Options
Reply | Threaded
Open this post in threaded view
|

Michigan Georef projection support?

Clever, Max

Hello Proj Developers,

 

I've read a few posts and it looks like Gerald has already developed support for RSO's.  In the previous emails shown below, I have code written for Hotine Oblique Mercators that uses the Michigan Georef Parameters.  Hopefully the code you have written is more accurate.  I could not get the Northing accurate enough.  If you want the parameters for Michigan Georef they are available online if you google "Michigan Georef Projection".  This projection is very important for any Michigan Proj users since most freely available data is provided in Michigan Georef.  I think the last time I saw it in a proj format it was listed as a regular omerc which is not correct even though many sites describe it that way. 

 

Anyhow, I was just wondering what the status was on the development of this projection

 

-Max Clever

 

 

 

 

 

----------------------------------------------------------------------------------------

Max,

 

OK, well implementing new projection methods is definately Gerald's area, not mine.  So if you can interest him in it, then things should

be good.   I just not competent to do it.

 

Best regards,

--

---------------------------------------+--------------------------------

---------------------------------------+------

I set the clouds in motion - turn up   | Frank Warmerdam, [hidden email]

light and sound - activate the windows | http://pobox.com/~warmerdam

and watch the world go round - Rush    | President OSGF, http://osgeo.org

 

-----Original Message-----
From: Clever, Max
Sent: Tuesday, July 18, 2006 9:26 AM
To: 'Frank Warmerdam'
Cc: [hidden email]
Subject: RE: [UMN_MAPSERVER-USERS] FW: Michigan Georef Projection Problems in Proj4

 

Frank,

 

First, to answer your question, the parameters Melita offered are not accurate enough and are not the true parameters of the Michigan Georef Projection.  Also, the translation from EPSG to PROJ4 is not correct.  I don't think Proj 4 supports a point azimuth method transformation for Hotine oblique Mercator (Rectified Skew Orthomorphic Projection)

 

Here is a link that contains the formulas for oblique mercator and Hotine Oblique Mercator (Michigan Georef is Hotine Oblique).

 

 http://www.posc.org/Epicentre.2_2/DataModel/ExamplesofUsage/eu_cs34i.html

 

please note though, that in the forward case, I believe the small "v" should be computed using Log not natural log. 

 

Also, if you want it for reference, below is the code of a couple of VB forms that I wrote 3 years ago with help from my professor in college to perform the forward and reverse cases.  It is only accurate to about 7 cm in the Northing and 3 mm in the Easting when converting a point from Michigan Georef to Geographic and back.  I think there must be some lack of precision in the code that is the reason for the accuracy problem.  I hope this helps.

 

 

Const a = 6378137                         ' semi major axis of the GRS80 ellipsoid

Const e2 = 0.0066943800229                ' eccentricity^2 of the GRS80 ellipsoid

Dim Pi As Double

Dim E As Double                           ' excentricity of the GRS80 ellipsoid

 

Public Function G2GREF(CoordIN() As Double) As Double()

    Dim M1 As Double ' Paramter

    Dim m2 As Double ' Paramter

    Dim t1 As Double ' Paramter

    Dim t2 As Double ' Paramter

    Dim t0 As Double ' Paramter

    Dim n As Double ' Paramter

    Dim F As Double ' Paramter

    Dim Eb As Double ' Paramter

    Dim Nb As Double ' Paramter

    Dim R As Double ' Paramter

   

    MS_GeoRef = "Projection Oblique Mercator; Datum:  NAD83; Ellipsoid: GRS80 " & vbCrLf & _

        "Standard Units: Meters" & vbCrLf & _

        "Origin = 86° 00' 00 W and 45° 18' 33 N" & vbCrLf & _

        "Scale factor at projection's center: k= 0.9996" & vbCrLf & _

        "Azimuth at center of projection: 337° 15' 20" & vbCrLf & _

        "False Easting: 2546731.496, False Northing: -4354009.816"

       

    Lon1 = -86#                 ' Longitude of projection's origin: 86° 00' 00" W

    Lon1 = Lon1 * Pi / 180#

    Lat1 = 45.30916666667       ' Latitude of projection's origin: 45° 18' 33" N

    Lat1 = Lat1 * Pi / 180#

    Az = 337.255555555556       ' Azimuth at center of projection: 337° 15' 20

    Az = Az * Pi / 180#

    SF = 0.9996                 ' Scale factor at projection's center

    Eb = 2546731.496            ' False Easting ( Eb = 528600.24)

    Nb = -4354009.816           ' False Northing (Nb = 499839.834)

   

    B = Sqr(1 + e2 * Cos(Lat1) ^ 4 / (1 - e2))

    A1 = a * B * SF * Sqr(1 - e2) / (1 - e2 * (Sin(Lat1) ^ 2))

    Temp = ((1 - E * Sin(Lat1)) / (1 + E * Sin(Lat1))) ^ (0.5 * E)

    t0 = Tan(Pi / 4 - (Lat1) / 2) / Temp

    D = B * Sqr(1 - e2) / (Cos(Lat1) * Sqr(1 - e2 * Sin(Lat1) ^ 2))

    F = D + Lat1 / Abs(Lat1) * Sqr(D ^ 2 - 1)       ' Eq. 4.110

    E1 = F * t0 ^ B                                 ' Eq. 4.111

    G = 0.5 * (F - 1 / F)                           ' Eq. 4.112

    Gamma0 = Isin(Sin(Az) / D)                      ' Eq. 4.113

    Lon0 = Lon1 - Isin(G * Tan(Gamma0)) / B  ' Lambda 0 at Eq. 4.114

    U0 = (Lat1 / Abs(Lat1)) * (A1 / B) * Atn(Sqr(D ^ 2 - 1) / Cos(Az))

    V0 = 0

   

    LatIN = DMS2R(CoordIN(1))

    LonIN = DMS2R(CoordIN(2))

   

    Temp = ((1 - E * Sin(LatIN)) / (1 + E * Sin(LatIN))) ^ (0.5 * E)

    t = Tan(Pi / 4 - (LatIN) / 2) / Temp            ' Eq. 4.117

    Q = E1 / (t ^ B)

    S = 0.5 * (Q - 1 / Q)

    Tc1 = 0.5 * (Q + 1 / Q)

    V1 = Sin(B * (LonIN - Lon0))

    U1 = (-V1 * Cos(Gamma0) + S * Sin(Gamma0)) / Tc1

    vl = A1 * Log((1 - U1) / (1 + U1)) / (2 * B)

    temp1 = (S * Cos(Gamma0) + V1 * Sin(Gamma0)) / Cos(B * (LonIN - Lon0))

    ul = (A1 / B) * Atn(temp1)

    x = vl * Cos(Az) + ul * Sin(Az) + Eb

    y = ul * Cos(Az) - vl * Sin(Az) + Nb

   

    ReDim Preserve CTemp(1 To 2) As Double

   

    CTemp(1) = x

    CTemp(2) = y

    G2GREF = CTemp

       

End Function

 

 

Public Function GREF2G(CoordIN() As Double) As Double()

    Dim M1 As Double ' Paramter

    Dim m2 As Double ' Paramter

    Dim t1 As Double ' Paramter

    Dim t2 As Double ' Paramter

    Dim t0 As Double ' Paramter

    Dim n As Double ' Paramter

    Dim F As Double ' Paramter

    Dim Eb As Double ' Paramter

    Dim Nb As Double ' Paramter

    Dim PhiOut As Double ' Paramter

    Dim LonOut As Double ' Paramter

   

    MS_GeoRef = "Projection Oblique Mercator; Datum:  NAD83; Ellipsoid: GRS80 " & vbCrLf & _

        "Standard Units: Meters" & vbCrLf & _

        "Origin = 86° 00' 00 W and 45° 18' 33 N" & vbCrLf & _

        "Scale factor at projection's center: k= 0.9996" & vbCrLf & _

        "Azimuth at center of projection: 337° 15' 20" & vbCrLf & _

        "False Easting: 2546731.496, False Northing: -4354009.816"

       

    Lon1 = -86#                 ' Longitude of projection's origin: 86° 00' 00" W

    Lon1 = Lon1 * Pi / 180#

    Lat1 = 45.309166666667       ' Latitude of projection's origin: 45° 18' 33" N

    Lat1 = Lat1 * Pi / 180#

    Az = 337.255555555556       ' Azimuth at center of projection: 337° 15' 20

    Az = Az * Pi / 180#

    SF = 0.9996                 ' Scale factor at projection's center

    Eb = 2546731.496            ' False Easting ( Eb = 528600.24)

    Nb = -4354009.816           ' False Northing (Nb = 499839.834)

   

   

    ttemp = e2 * (Cos(Lat1) ^ 4) / (1 - e2)

    B = Sqr(1 + ttemp)

    A1 = a * B * SF * Sqr(1 - e2) / (1 - e2 * (Sin(Lat1)) ^ 2)

    Temp = ((1 - E * Sin(Lat1)) / (1 + E * Sin(Lat1))) ^ ( 0.5 * E)

    t0 = Tan(Pi / 4 - (Lat1) / 2) / Temp

    ttemp1 = Cos(Lat1) * Sqr(1 - e2 * (Sin(Lat1)) ^ 2)

    D = B * Sqr(1 - e2) / ttemp1

   

    F = D + Lat1 / Abs(Lat1) * Sqr(D ^ 2 - 1)       ' Eq. 4.110

    E1 = F * t0 ^ B                                 ' Eq. 4.111

    G = 0.5 * (F - 1 / F)                           ' Eq. 4.112

    Gamma0 = Isin(Sin(Az) / D)                      ' Eq. 4.113

    Lon0 = Lon1 - Isin(G * Tan(Gamma0)) / B  ' Lambda 0 at Eq. 4.114

   

    xIN = CoordIN(1)

    yIN = CoordIN(2)

 

'Actual Computations for Reverse case Hotine Oblique Mercator

 

    xr = xIN - Eb

    yr = yIN - Nb

   

    vs = xr * Cos(Az) - yr * Sin(Az)

    us = yr * Cos(Az) + xr * Sin(Az)

   

    temp1 = -B * vs / A1

    Qp = (2.71828182845905) ^ temp1

    Sp = 0.5 * (Qp - 1 / Qp)

    Tp = 0.5 * (Qp + 1 / Qp)

   

    Vp = Sin(B * us / A1)

    Up = (Vp * Cos(Gamma0) + Sp * Sin(Gamma0)) / Tp

 

    temp2 = (1 + Up) / (1 - Up)

    t = (E1 / Sqr(temp2)) ^ (1 / B)

Threaded

Open this post in threaded view
|

Re: Michigan Georef projection support?

OvV_HN
From: "Clever, Max" <Maxc-spicergroup.com>
Date: Thu, 20 Jul 2006 10:23:22 -0400
Subject: [Proj] Michigan Georef projection support?

> I've read a few posts and it looks like Gerald has already developed support
> for RSO's.  In the previous emails shown below, I have code written for Hotine
> Oblique Mercators that uses the Michigan Georef Parameters. Hopefully the code
> you have written is more accurate.  I could not get the Northing accurate
> enough.  If you want the parameters for Michigan Georef they are available
> online if you google "Michigan Georef Projection".

I still haven't seen any (near) official test points (lat,lon to x,y) for
Michigan GeoRef. Not in the maptools postings, not on the internet.
If nobody will produce some serious test points, no problem can be
diagnosed.

The test points M. Kennedy gave on 2004-12-20 in a posting can be understood
as follows.
The parameters are for my code, a slightly improved version of libproj4's
omerc with additions from the EPSG/OGP Guidance Note.

lat = 43; lon = -86;
// GRS80 ellipsoid
type = 5; // Hotine-RSO according to EPSG
lat0 = 45.30916666667; // latitude of projection's origin: 45d18m33s N
lonc = -86; // longitude of projection's origin: 86d00m00s W
alfac = 337.255555555556; // azimuth at center of projection: 337d15m20s
gamma = ""; // let the code calculate gamma
k0 = 0.9996; // scale factor at projection's center
x0 = 2546731.496; // false easting
y0 = -4354009.816; // false northing
norot = 0; no_off = 0;
omerc(type,lat,lon,lat0,lonc,alfac,gamma,k0,x0,y0,norot,no_off)
Result: x,y = 499864.500077, 272108.859263

Some people get mighty nervous seeing all those decimals, but for testing
purposes they are necessary.

M. Kennedy, 2004-12-20:
"
and now the results from the ESRI Projection Engine
cymru{melita}: forward91 102123
Projection Engine Version 9.0 (Dec 15 2004)
499864.50  272108.86
"

And now for something completely different:
(only the differences with above are mentioned):

type = 1; // comparable with libproj4 omerc or the USA-type 'Hotine'
no_off = 1;
Result: 499864.500077, 272108.859263

As long as no official test points are available I would suggest:
Do the MI GeoRef with (lib)proj omerc and the no_offset flag set.




_______________________________________________
Proj mailing list
[hidden email]
http://lists.maptools.org/mailman/listinfo/proj