<HTML dir=ltr><HEAD><TITLE>[Proj] Using Proj.4 to compress a 10-parameter Molodensky-Badekas datum shift into 7 parameters</TITLE>
<META http-equiv=Content-Type content="text/html; charset=unicode">
<META content="MSHTML 6.00.6000.17080" name=GENERATOR></HEAD>
<BODY>
<DIV id=idOWAReplyText72307 dir=ltr>
<DIV dir=ltr><FONT face="Times New Roman" color=#000000 size=2>The "evaluation point" is the Datum Origin, a rather important item.&nbsp; In regard to the "unexplained" difference in ellipsoid height results, my guess is the approximations made by the inverse algorithm from Geocentric coordinates.</FONT></DIV>
<DIV dir=ltr><FONT size=2></FONT>&nbsp;</DIV>
<DIV dir=ltr><FONT size=2>Personally, I prefer the elegance of incorporating the Datum Origin.&nbsp; </FONT></DIV>
<DIV dir=ltr><FONT size=2></FONT>&nbsp;</DIV>
<DIV dir=ltr><FONT size=2>I have also been told by Dr. Muneendra Kumar that Badekas had nothing to do with Molodensky's original&nbsp;method.&nbsp; I believe they were in school together at the same time, and Prof. Molodensky preceeded them on another continent by some years if not decades.</FONT></DIV>
<DIV dir=ltr><FONT face="Times New Roman" color=#000000 size=2></FONT>&nbsp;</DIV></DIV>
<DIV id=idSignature85798 dir=ltr>
<DIV><FONT face="Times New Roman" color=#000000 size=2><SPAN style="FONT-SIZE: 10pt">
<DIV class=MsoNormal style="MARGIN: 0in 0in 0pt"><SPAN lang=DE style="FONT-SIZE: 10pt; mso-ansi-language: DE">Clifford J. Mugnier, C.P., C.M.S.<?xml:namespace prefix = o ns = "urn:schemas-microsoft-com:office:office" /><o:p></o:p></SPAN></DIV>
<DIV class=MsoNormal style="MARGIN: 0in 0in 0pt"><SPAN style="FONT-SIZE: 10pt">Chief of Geodesy,<o:p></o:p></SPAN></DIV>
<DIV class=MsoNormal style="MARGIN: 0in 0in 0pt"><B><SPAN style="FONT-SIZE: 10pt; FONT-VARIANT: small-caps">Center for GeoInformatics<o:p></o:p></SPAN></B></DIV>
<DIV class=MsoNormal style="MARGIN: 0in 0in 0pt"><SPAN style="FONT-SIZE: 10pt">Department of Civil Engineering <o:p></o:p></SPAN></DIV>
<DIV class=MsoNormal style="MARGIN: 0in 0in 0pt"><SPAN style="FONT-SIZE: 10pt">Patrick F. Taylor Hall 3223A<o:p></o:p></SPAN></DIV>
<DIV class=MsoNormal style="MARGIN: 0in 0in 0pt"><B><SPAN style="FONT-SIZE: 10pt">LOUISIANA STATE UNIVERSITY <o:p></o:p></SPAN></B></DIV>
<DIV class=MsoNormal style="MARGIN: 0in 0in 0pt"><SPAN style="FONT-SIZE: 10pt">Baton Rouge, LA<SPAN style="mso-spacerun: yes">&nbsp; </SPAN>70803<o:p></o:p></SPAN></DIV>
<DIV class=MsoNormal style="MARGIN: 0in 0in 0pt"><SPAN style="FONT-SIZE: 10pt">Voice and Facsimile:<SPAN style="mso-spacerun: yes">&nbsp; </SPAN>(225) 578-8536 [Academic] <o:p></o:p></SPAN></DIV>
<DIV class=MsoNormal style="MARGIN: 0in 0in 0pt"><SPAN style="FONT-SIZE: 10pt">Voice and Facsimile:<SPAN style="mso-spacerun: yes">&nbsp; </SPAN>(225) 578-4578 [Research] <o:p></o:p></SPAN></DIV>
<DIV class=MsoNormal style="MARGIN: 0in 0in 0pt"><SPAN style="FONT-SIZE: 10pt">Cell: (225) 238-8975 [Academic &amp; Research]<o:p></o:p></SPAN></DIV>
<DIV class=MsoNormal style="MARGIN: 0in 0in 0pt"><SPAN style="FONT-SIZE: 10pt">Honorary Life Member of the <o:p></o:p></SPAN></DIV>
<DIV class=MsoNormal style="MARGIN: 0in 0in 0pt"><SPAN style="FONT-SIZE: 10pt">Louisiana Society of Professional Surveyors <o:p></o:p></SPAN></DIV>
<DIV class=MsoNormal style="MARGIN: 0in 0in 0pt"><SPAN style="FONT-SIZE: 10pt">Fellow Emeritus of the ASPRS <o:p></o:p></SPAN></DIV>
<DIV class=MsoNormal style="MARGIN: 0in 0in 0pt"><SPAN style="FONT-SIZE: 10pt">Member of the Americas Petroleum Survey Group<o:p></o:p></SPAN></DIV><BR></SPAN></FONT></DIV></DIV>
<DIV dir=ltr><BR>
<HR tabIndex=-1>
<FONT face=Tahoma size=2><B>From:</B> proj-bounces@lists.maptools.org on behalf of Mikael Rittri<BR><B>Sent:</B> Tue 12-Oct-10 03:45<BR><B>To:</B> PROJ.4 and general Projections Discussions<BR><B>Subject:</B> [Proj] Using Proj.4 to compress a 10-parameter Molodensky-Badekas datum shift into 7 parameters<BR></FONT><BR></DIV>
<DIV><BR>
<P><FONT size=2>A Molodensky-Badekas datum shift is represented by 10 parameters.<BR>But it can be compressed to a 7-parameter representation.<BR><BR>As Noel Zinn has explained (<A href="http://lists.maptools.org/pipermail/proj/2010-October/005424.html">http://lists.maptools.org/pipermail/proj/2010-October/005424.html</A>),<BR>the 10-parameter representation has advantages over the 7-parameter<BR>one, for people who derive, analyze and compare datum shifts.<BR><BR>But for people who just use a datum shift to transform coordinates,<BR>the 7-parameter representation gives the same results. In this<BR>limited sense, the compression is lossless, which is useful if<BR>your software (like Proj.4) does not support 10-parameter M-B.<BR><BR>I learned how to do the compression from Melita Kennedy (thanks<BR>again!).&nbsp; And I have now realized that one can persuade Proj.4 to<BR>do most of the work of the compression process.&nbsp; Since I thought<BR>this is rather nice, I'd like to explain it by an example.<BR><BR>Example: M-B datum shift EPSG:1096, "La Canoa to WGS 84 (2)",<BR>which has the same parameters as EPSG:1771, "La Canoa to REGVEN (1)".<BR>&nbsp;&nbsp;&nbsp;<BR>&nbsp;&nbsp;&nbsp; Source ellipsoid: International 1924.<BR>&nbsp;&nbsp;&nbsp; Target ellipsoid: WGS 84.<BR>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;<BR>&nbsp;&nbsp;&nbsp; dX = -270.933 m<BR>&nbsp;&nbsp;&nbsp; dY =&nbsp; 115.599 m<BR>&nbsp;&nbsp;&nbsp; dZ = -360.226 m<BR>&nbsp;&nbsp;&nbsp; RX = -5.266 arc sec<BR>&nbsp;&nbsp;&nbsp; RY = -1.238 arc sec<BR>&nbsp;&nbsp;&nbsp; RZ =&nbsp; 2.381 arc sec<BR>&nbsp;&nbsp;&nbsp; dS = -5.109 ppm<BR>&nbsp;&nbsp;&nbsp; eX =&nbsp; 2464351.59 m<BR>&nbsp;&nbsp;&nbsp; eY = -5783466.61 m<BR>&nbsp;&nbsp;&nbsp; eZ =&nbsp;&nbsp; 974809.81 m<BR><BR>where eX,eY,eZ are the geocentric coordinates of the so-called<BR>evaluation point.<BR><BR>In principle, the equivalent 7-parameter datum shift should have<BR>the same rotations RX, RY, RZ and the same scale difference dS.<BR>But there is a catch: the sign convention for the rotations.<BR>M-B uses the same sign convention as Coordinate Frame Rotation,<BR>while cs2cs uses the opposite convention of the Position Vector<BR>Transform, so we must negate the three rotations.&nbsp; So, in cs2cs<BR>syntax, the 7-parameter datum shift we seek has the form<BR><BR>&nbsp;&nbsp;&nbsp; +towgs84=?,?,?,5.266,1.238,-2.381,-5.109<BR>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;<BR>To find the question marks, we have to do three steps.&nbsp;<BR><BR>i) Make a 7-parameter "something" by temporarily replacing<BR>&nbsp;&nbsp; the question marks by zeros.&nbsp; It is not really a datum<BR>&nbsp;&nbsp; shift, but we might call it a 7-parameter compressor.<BR>&nbsp;&nbsp; Use it in a cs2cs command like this<BR>&nbsp;&nbsp;<BR>&nbsp;&nbsp; &gt;cs2cs -f "%.3f" +proj=geocent +towgs84=0,0,0,5.226,1.238,-2.381,-5.109 +to +datum=WGS84 +proj=geocent<BR><BR>&nbsp;&nbsp; As input, use the evaluation point eX,eY,eZ of the M-B datum shift:<BR><BR>&nbsp;&nbsp; 2464351.59&nbsp;&nbsp; -5783466.61&nbsp;&nbsp; 974809.81<BR><BR>&nbsp;&nbsp; We get the output<BR><BR>&nbsp;&nbsp; 2464278.090&nbsp; -5783490.207&nbsp; 974643.507<BR><BR>ii) Now subtract the output row from the input row, and we get<BR><BR>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; 73.500&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; 23.597&nbsp;&nbsp;&nbsp;&nbsp; 166.303<BR><BR>iii) Then add the original dX,dY,dZ from the M-B datum shift,<BR>&nbsp;&nbsp;&nbsp;&nbsp; which were -270.933,115.599,-360.226, and we get<BR><BR>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; -197.433&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; 139.196&nbsp;&nbsp;&nbsp; -193.923<BR>&nbsp;&nbsp;<BR>The last three values are the question marks in our 7-parameter<BR>datum shift.<BR><BR>To test it, we use the test point given in section 2.4.4.1,<BR>pages 114-115, EPSG Guidance Note 7.2, <A href="http://www.epsg.org/guides/docs/G7-2.pdf">http://www.epsg.org/guides/docs/G7-2.pdf</A><BR><BR>&gt;cs2cs +ellps=intl +proj=longlat +towgs84=-197.433,139.196,-193.923,5.226,1.238,-2.381,-5.109 +to +datum=WGS84 +proj=longlat<BR>66d4'48.091"W&nbsp;&nbsp; 9d35'0.386"N&nbsp;&nbsp;&nbsp; 201.46<BR>66d4'54.705"W&nbsp;&nbsp; 9d34'49.001"N&nbsp;&nbsp; 180.499<BR><BR>According to GN 7.2, the result should have been&nbsp;<BR><BR>66°04'54.705"W&nbsp; 9°34'49.001"N&nbsp;&nbsp; 180.51 m<BR><BR>so we have achieved a lossless compression - apart from the<BR>unexplained difference of about 1.1 cm in the ellipoidal height.&nbsp;<BR>&nbsp;<BR>Regards,<BR>Mikael Rittri<BR>Carmenta AB<BR>Sweden<BR>www.carmenta.com<BR>_______________________________________________<BR>Proj mailing list<BR>Proj@lists.maptools.org<BR><A href="http://lists.maptools.org/mailman/listinfo/proj">http://lists.maptools.org/mailman/listinfo/proj</A><BR></FONT></P></DIV></BODY></HTML>