[Proj] Michigan Georef projection support?

Thu Jul 20 07:23:22 PDT 2006

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, warmerdam at pobox.com

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: MAPSERVER-USERS at LISTS.UMN.EDU
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)

    PhiOut = Pi / 2 - 2 * Atn(t)

'Iterative Solution for Phi Out

'    For i = 1 To 30

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

'        PhiOut = Pi / 2 - 2 * Atn(t * Temp)

'    Next i

'Single Line Solution for Phi Out

    PhiOut = PhiOut + Sin(2 * PhiOut) * (e2 / 2 + (5 * e2 ^ 2) / 24 + (e2 ^ 4) / 12 + (13 * e2 ^ 6) / 360) + Sin(4 * PhiOut) * ((7 * e2 ^ 2) / 48 + (29 * e2 ^ 4) / 240 + (811 * e2 ^ 6) / 11520) + Sin(6 * PhiOut) * ((7 * e2 ^ 4) / 120 + (81 * e2 ^ 6) / 1120) + Sin(8 * PhiOut) * ((4279 * e2 ^ 6) / 161280) 

    temp3 = (Sp * Cos(Gamma0) - Vp * Sin(Gamma0)) / Cos(B * us / A1)

    LonOut = Lon0 - Atn(temp3) / B

    ReDim Preserve CTemp(1 To 2) As Double

    CTemp(1) = r2dms(PhiOut) 

    CTemp(2) = r2dms(LonOut)

    GREF2G = CTemp

End Function

-----Original Message-----

From: Frank Warmerdam [mailto:fwarmerdam at gmail.com] On Behalf Of Frank Warmerdam

Sent: Saturday, July 15, 2006 10:01 PM

To: Clever, Max

Subject: Re: [UMN_MAPSERVER-USERS] FW: Michigan Georef Projection Problems in Proj4

Clever, Max wrote:

> Hi,

> 

>  

> 

> Did anyone see this the last time I sent it?  It relates to Mapserver as 

> well since Mapserver uses Proj4 for its projections. 

Max,

I have skimmed this material, but frankly I'm not sure what action item

there is, and I find myself with limited time for work on PROJ.4.

Is the problem that the parameters Melita offered a couple years ago

aren't accurate enough?  Or that the underlying translation from EPSG

to PROJ.4 wasn't fixed so the epsg file keeps getting regenerated

wrong?  If it is a tranlation problem, then that is basically something

I ought to fix.  But I basically need some formulation to recognise

a distinct oblique mercator case for the michigan projection from the

EPSG codes (or parameters), and what that should map to in WKT format,

and in PROJ.4 format.  I'm happy to use the ESRI WKT representation if

there isn't an obvious existing form for this special case.

If you can walk me through what should be changed, I'm willing to work

on it.

Best regards,

...

> Two years ago, I ran into a problem with the Michigan Georef Projection 

> and the way that proj identified it.  I had sent emails back and forth 

> for a while until someone sent a temporary solution of providing false 

> parameters that worked.. *for the most part*.  This temporary solution, 

> of course, did not actually solve the problem, but instead delayed the 

> fixing of the methods that proj identifies projections and translates 

> them.  For that I am sorry for not remaining vigilant in seeing a true 

> solution being devised.  But now, since I have just now installed the 

> latest version of GRASS 6.1, I have come full circle and face this 

> problem again.  To provide a quick access to the background of what has 

> already been said on this projection please note the previous emails 

> below.  I believe, at this time still, that the *omerc* projection and 

> its parameters as used by proj cannot correctly describe or transform a 

> omerc projection with a "natural origin".  From what I understand, 

> Hotine oblique mercator and Rectified Skew Orthomorphic are one and the 

> same or depend on where the skew is corrected.  Has there been a 

> solution determined for this projection?  If not, maybe a solution to 

> this problem would be to have Proj have oblique mercators split between 

> "natural origins" and cartesian center point origins.  I hope, maybe, 

> someone has been looking at this lately but I doubt it.  Any comments or 

> solutions would be very welcome.  Thanks.

-- 

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

I set the clouds in motion - turn up   | Frank Warmerdam, warmerdam at pobox.com

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

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

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/proj/attachments/20060720/850d84aa/attachment.html>