[Proj] Michigan Georef Projection Problems

Clever, Max Maxc at spicergroup.com
Sat Jul 22 23:44:49 PDT 2006


Melita,
 
are you saying that since December of 2004 Proj has been corrected for the transformation of Michigan Georef?  I have downloaded what I think is the most up-to-date version of Proj4 when I installed GRASS cvs 6.1, and it still did not correctly project my data that was in Michigan Georef.  
 
Otherwise, just to refresh the subject, please read what you wrote 2 years ago on December 20th, 2004 to Gerald Evenden and others, where you basically provided manipulated parameters for the EPSG of Michigan Georef in order to make it translate to approximately the correct location.
 
Also, from what I have read in an email from Anthony Dunk dated May 18, 2006, the Alaska Zone 1, as defined in its EPSG, is not being transformed properly either.  It sounded like someone got the correct answer when they set a flag called "no_ofs".  Is this a setting that can be added to the EPSG file?  
 
-----Original Message----- 
From: proj-bounces at lists.maptools.org on behalf of proj-request at lists.maptools.org 
Sent: Sat 7/22/2006 12:00 PM 
To: proj at lists.maptools.org 
Cc: 
Subject: Proj Digest, Vol 26, Issue 13



	Send Proj mailing list submissions to
	        proj at lists.maptools.org
	
	To subscribe or unsubscribe via the World Wide Web, visit
	        http://lists.maptools.org/mailman/listinfo/proj
	or, via email, send a message with subject or body 'help' to
	        proj-request at lists.maptools.org
	
	You can reach the person managing the list at
	        proj-owner at lists.maptools.org
	
	When replying, please edit your Subject line so it is more specific
	than "Re: Contents of Proj digest..."
	
	
	Today's Topics:
	
	   1. Re: Michigan GeoRef (Melita Kennedy)
	   2. Michigan Georef Projection Problems (Clever, Max)
	   3. Re: FW: [Proj] First attempt of porting PROJ.4 to Windows CE
	      (Mateusz Loskot)
	   4. Michael P Finn of the USGS/ National Geospatial Technical
	      Operations Center is out of the office. (Michael P Finn)
	
	
	----------------------------------------------------------------------
	
	Message: 1
	Date: Fri, 21 Jul 2006 14:30:26 -0700 (GMT-07:00)
	From: Melita Kennedy <mkennedy2 at earthlink.net>
	Subject: [Proj] Re: Michigan GeoRef
	To: proj at lists.maptools.org
	Message-ID:
	        <4898027.1153517426540.JavaMail.root at elwamui-karabash.atl.sa.earthlink.net>
	       
	Content-Type: text/plain; charset=ISO-8859-1
	
	Michigan GeoRef uses the same projection as State Plane Alaska zone 1.
	I don't know which projection this is in PROJ. We shouldn't need to implement
	another implementation of Hotine oblique Mercator to support Michigan GeoRef.
	
	These are the official parameters from the Michigan DNR:
	
	Projection: Oblique Mercator
	Datum: NAD83
	Ellipsoid: GRS80
	Standard Units: Meters
	Scale factor at projection's center: 0.9996
	Longitude of projection's origin: 86° 00' 00" W
	Latitude of projection's origin: 45° 18' 33" N
	Azimuth at center of projection: 337.25556
	False Easting: 2546731.496
	False Northing: -4354009.816
	
	Please note the azimuth value. The value in decimal degrees is not equal to the original
	definition of 337d 15m 20s. Michigan GeoRef was originally defined to work with ArcInfo
	workstation which had a limitation in how many decimal places were allowed in the azimuth
	parameter. We've changed the ArcGIS version of the definition to match this--when we
	added the definition to the projection engine library, we used the DMS value--because
	Michigan clients were complaining that data wasn't lining up properly.
	
	Melita
	
	--
	Melita Kennedy
	ESRI Product Specialist
	mkennedy at esri.com (but posting from mkennedy2 at earthlink.net)
	
	-----Original Message-----
	>Date: Thu, 20 Jul 2006 10:23:22 -0400
	>From: "Clever, Max" <Maxc at spicergroup.com>
	>Subject: [Proj] Michigan Georef projection support?
	>To: <proj at lists.maptools.org>
	>
	>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
	
	
	
	
	------------------------------
	
	Message: 2
	Date: Fri, 21 Jul 2006 23:19:01 -0400
	From: "Clever, Max" <Maxc at spicergroup.com>
	Subject: [Proj] Michigan Georef Projection Problems
	To: <proj at lists.maptools.org>
	Message-ID:
	        <B5DFCBECD06C8246A4B1279EC6487B740837CE at sag-x225-exchg0.spicergroup.com>
	       
	Content-Type: text/plain;       charset="utf-8"
	
	>From reading the previous posts on proj's support for RSO's it sounds like it is already in the works.  Just to make sure though, I just want to let everyone know that the Michigan Georef Projection is not supported by proj.  The only proj definition that I have seen for it is the one described under the ESRI listing that describes it as an <omerc> which is not really correct. 
	
	Michigan Georef is a Hotine Oblique Mercator.  The official parameters are listed at this website http://www.michigan.gov/documents/DNR_Map_Proj_and_MI_Georef_Info_20889_7.pdf  The parameters as given at this site require the computation of gamma, hopefully that helps you figure out whether to have gamma as a mandatory or optional parameter requirement.
	
	I appreciate any help that could be given.  Please see the history of emails below for reference to this problem.  Thanks.
	
	-Max
	
	
	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 <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          ' excentricity^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 <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 <http://pobox.com/~warmerdam>
	and watch the world go round - Rush    | President OSGF, http://osgeo.org <http://osgeo.org>
	
	
	
	
	------------------------------
	
	Message: 3
	Date: Sat, 22 Jul 2006 13:22:24 +0200
	From: Mateusz Loskot <mateusz at loskot.net>
	Subject: Re: FW: [Proj] First attempt of porting PROJ.4 to Windows CE
	To: "PROJ.4 and general Projections Discussions"
	        <proj at lists.maptools.org>
	Message-ID: <44C20A70.2020600 at loskot.net>
	Content-Type: text/plain; charset=ISO-8859-1
	
	Wynand Kruger wrote:
	> Where can we find 'proj4-wince-port.zip'?
	
	Hi Wynand,
	
	I replied to your first e-mail in you're asking about it.
	Please, take a look at:
	http://lists.maptools.org/pipermail/proj/2006-July/002406.html
	
	Best regards
	--
	Mateusz Loskot
	http://mateusz.loskot.net
	
	
	------------------------------
	
	Message: 4
	Date: Sat, 22 Jul 2006 10:01:25 -0500
	From: Michael P Finn <mfinn at usgs.gov>
	Subject: [Proj] Michael P Finn of the USGS/ National Geospatial
	        Technical Operations Center is out of the office.
	To: "PROJ.4 and general Projections Discussions"
	        <proj at lists.maptools.org>
	Message-ID:
	        <OFDAAE9684.8189DBFE-ON862571B3.005286F5-862571B3.005286F5 at usgs.gov>
	Content-Type: text/plain; charset=US-ASCII
	
	
	I will be out of the office starting  07/22/2006 and will not return until
	07/31/2006.
	
	I will be unable to check e-mail regularly while I am out. If I can't
	respond to your message while I am out, then I will when I return. If you
	have anything that cannot wait, please contact my colleague, Mr. Timothy
	Krupinski, at tkrupinski at usgs.gov.     Mike
	
	
	
	------------------------------
	
	_______________________________________________
	Proj mailing list
	Proj at lists.maptools.org
	http://lists.maptools.org/mailman/listinfo/proj
	
	End of Proj Digest, Vol 26, Issue 13
	************************************
	
	

-------------- next part --------------
A non-text attachment was scrubbed...
Name: winmail.dat
Type: application/ms-tnef
Size: 25878 bytes
Desc: not available
URL: <http://lists.osgeo.org/pipermail/proj/attachments/20060723/c0d2430e/attachment.bin>


More information about the Proj mailing list