Implementation for WGS84 (lat,lon) conversion to Greek Grid (GGRS87)

In my last  post I said that I had found a site with an algorithm that would convert WGS84 (lat/lon) coordinates to the Greek Grid system (GGRS87).  I successfully converted that algorithm into a .php program and I present it here.  The lines commented out are the lines from the algorithm.  The immediately following line is my php interpretation of this.  As it is it rounds down to the whole meter.  But its accuracy appears to be of the order of one one-thousandth of a meter (millimeter).  Plenty good enough for the likes of us.


<?php

// from algorithm at: http://studentguru.gr/b/kokyri/archive/2012/02/21/convert-coordinates-from-greek-grid-ggrs87-87-espg-2100-to-latitude-longitude-wsg84

//public static double[] ToGreekGrid(double lat, double lng)
function ToGreekGrid($lat, $lng)
//        {
{
//            lat = lat * Math.PI / 180;
$lat = $lat * M_PI / 180.0;
//            lng = lng * Math.PI / 180;
$lng = $lng * M_PI / 180.0;
//            var a = 6378137;
$a = 6378137;
//            var b = 6356752.31424518;
$b = 6356752.31424518;
//            var e2 = (Math.Pow(a,2) - Math.Pow(b,2)) / Math.Pow(a, 2);
$e2 = ($a * $a - $b * $b) / ($a * $a);
//            var v = a / Math.Sqrt(1 - e2 * (Math.Sin(lat) * Math.Sin(lat)));
$v = $a / sqrt(1 - $e2 * (sin($lat) * sin($lat)));
//            var X = v * Math.Cos(lat) * Math.Cos(lng);
$X = $v * cos($lat) * cos($lng);
//            var Y = v * Math.Cos(lat) * Math.Sin(lng);
$Y = $v * cos($lat) * sin($lng);
//            var Z = (v * (1 - e2)) * Math.Sin(lat);
$Z = ($v * (1.0 - $e2)) * sin($lat);
//            var px = 199.723;
$px = 199.723;
//            var py = -74.03;
$py = -74.03;
//            var pz = -246.018;
$pz = -246.018;
//            var rw = 0;
$rw = 0;
//            var rf = 0;
$rf = 0;
//            var rk = 0;
$rk = 0;
//            var ks = 1;
$ks = 1;
//            var c1 = Math.Cos(rw);
$c1 = cos($rw);
//            var c2 = Math.Cos(rf);
$c2 = cos($rf);
//            var c3 = Math.Cos(rk);
$c3 = cos($rk);
//            var s1 = Math.Sin(rw);
$s1 = sin($rw);
//            var s2 = Math.Sin(rf);
$s2 = sin($rf);
//            var s3 = Math.Sin(rk);
$s3 = sin($rk);
//            var D11 = c2 * c3;
$D11 = $c2 * $c3;

//            var D21 = -c2 * s3;
$D21 = $c2*(-1) * $s3;

//            var D31 = Math.Sin(rf);
$D31 = sin($rf);

//            var D12 = s1 * s2 * c3 + c1 * s3;
$D12 = $s1 * $s2 * $c3 + $c1 * $s3;

//            var D22 = -s1 * s2 * s3 + c1 * c3;
$D22 = $s1*(-1) * $s2 *$s3 + $c1 * $c3;

//            var D32 = -s1 * c2;
$D32 = $s1*(-1) * $c2;

//            var D13 = -c1 * s2 * c3 + s1 * s3;
$D13 = $c1*(-1) * $s2 * $c3 + $s1 * $s3;

//            var D23 = c1 * s2 * s3 + s1 * c3;
$D23 = $c1 * $s2 * $s3 + $s1 * $c3;

//            var D33 = c1 * c2;
$D33 = $c1 * $c2;

//            var X1 = px + ks * (D11 * X + D12 * Y + D13 * Z);
$X1 = $px + $ks * ($D11 * $X + $D12 * $Y + $D13 * $Z);

//            var Y1 = py + ks * (D21 * X + D22 * Y + D23 * Z);
$Y1 = $py + $ks * ($D21 * $X + $D22 * $Y + $D23 * $Z);

//            var Z1 = pz + ks * (D31 * X + D32 * Y + D33 * Z);
$Z1 = $pz + $ks * ($D31 * $X + $D32 * $Y + $D33 * $Z);

//            lng = Math.Atan2(Y1, X1);
$lng = atan2($Y1, $X1);

//            var lat0 = Math.Atan2(Z1, Math.Sqrt(X1 * X1 + Y1 * Y1));
$lat0 = atan2($Z1, sqrt($X1 * $X1 + $Y1 * $Y1));

//            b = 6356752.31414036;
$b = 6356752.31414036;

//            e2 = (Math.Pow(a, 2) - Math.Pow(b,2)) / Math.Pow(a, 2);
$e2 = ($a * $a - $b * $b) / ($a * $a);

//            while (Math.Abs(lat - lat0) > 0.0000000001)
while (abs($lat - $lat0) > 0.0000000001)
//            {
{

//                var No = a / Math.Sqrt(1 - e2 * Math.Sin(lat0) * Math.Sin(lat0));
$No = $a / sqrt(1 - $e2 * sin($lat0) * sin($lat0));


//                var h = Math.Sqrt(X1 * X1 + Y1 * Y1) / Math.Cos(lat0) - No;
$h = sqrt($X1 * $X1 + $Y1 * $Y1) / cos($lat0) - $No;

//                lat = lat0;
$lat = $lat0;

//                lat0 = Math.Atan(Z1 / Math.Sqrt(X1 * X1 + Y1 * Y1) * (1 / (1 - e2 * No / (No+ h))));
$lat0 = atan($Z1 /sqrt($X1 * $X1 + $Y1 * $Y1) * (1 / (1 - $e2 * $No / ($No + $h))));
//            }
}
//            lng = lng - 24 * Math.PI / 180;
$lng = $lng - 24 * M_PI / 180;

//            var m0 = 0.9996;
$m0 = 0.9996;

//            var es2 = (Math.Pow(a,2) - Math.Pow(b,2)) / (Math.Pow(b,2));
$es2 = ($a * $a - $b * $b) / ($b * $b);

//            var V = Math.Sqrt(1 + es2 * Math.Cos(lat) * Math.Cos(lat));
$V = sqrt(1 + $es2 * cos($lat) * cos($lat));

//            var eta = Math.Sqrt(es2 * Math.Cos(lat) * Math.Cos(lat));
$eta = sqrt($es2 * cos($lat) * cos($lat));

//            var Bf = Math.Atan(Math.Tan(lat) / Math.Cos(V * lng) * (1 + eta * eta / 6 * (1 - 3 * Math.Sin(lat) * Math.Sin(lat)) * lng * lng * lng * lng));
                     $Bf =  atan(tan($lat) / cos($V * $lng) * (1 + $eta * $eta / 6 * (1 - 3 * sin($lat) * sin($lat)) * $lng * $lng * $lng * $lng));

//            var Vf = Math.Sqrt(1 + es2 * Math.Cos(Bf) * Math.Cos(Bf));
$Vf = sqrt(1 + $es2 * cos($Bf) * cos($Bf));

//            var etaf = Math.Sqrt(es2 * Math.Cos(Bf) * Math.Cos(Bf));
$etaf = sqrt($es2 * cos($Bf) * cos($Bf));

//            var n = (a - b) / (a + b);
$n = ($a - $b) / ($a + $b);

//            var r1 = (1 + n * n / 4 + n * n * n * n / 64) * Bf;
$r1 = (1 + $n * $n / 4 + $n * $n * $n * $n / 64) * $Bf;

//            var r2 = 3.0 / 2.0 * n * (1 - n * n / 8) * Math.Sin(2 * Bf);
$r2 = 3.0 / 2.0 * $n * (1 - $n * $n / 8) * sin(2 * $Bf);

//            var r3 = 15.0 / 16.0 * n * n * (1 - n * n / 4) * Math.Sin(4 * Bf);
$r3 = 15.0 / 16.0 * $n * $n * (1 - $n * $n / 4) * sin (4 * $Bf);

//            var r4 = 35.0 / 48.0 * n * n * n * Math.Sin(6 * Bf);
$r4 = 35.0 / 48.0 * $n * $n * $n * sin(6 * $Bf);

//            var r5 = 315.0 / 512.0 * n * n * n * n * Math.Sin(8 * Bf);
$r5 = 315.0 / 512.0 * $n * $n * $n * $n * sin(8 * $Bf);

//            var Northing = a / (1 + n) * (r1 - r2 + r3 - r4 + r5) * m0 - 0.001 + 4202812 - 4207988.1206046063;
$Northing = $a / (1 + $n) * ($r1 - $r2 + $r3 - $r4 + $r5) * $m0 - 0.001 + 4202812 - 4207988.1206046063;

//            var ys = Math.Tan(lng) * Math.Cos(Bf) / Vf * (1 + etaf * etaf * lng * lng * Math.Cos(Bf) * Math.Cos(Bf) * (etaf * etaf / 6 + lng * lng / 10));

$ys = tan($lng) * cos($Bf) / $Vf * (1 + $etaf * $etaf * $lng * $lng * cos($Bf) * cos($Bf) * ($etaf * $etaf / 6 + $lng * $lng / 10));

//            ys = Math.Log(ys + Math.Sqrt(ys * ys + 1));
$ys = log($ys + sqrt($ys * $ys + 1));

//            var Easting = m0 * Math.Pow(a,2) / b * ys + 500000;
$Easting = $m0 * pow($a, 2) / $b * $ys + 500000;

$Northing = $Northing + 5176.122; // Added correction per comments

// echo $Easting.PHP_EOL;
// echo $Northing.PHP_EOL.PHP_EOL;

return array(0 => round($Easting), 1 => round($Northing));

//        }
}

/// main

//ToGreekGrid (37.97412198000, 22.59898503000);

list ($Easting, $Northing) = ToGreekGrid(37.97412198000, 22.59898503000);

echo PHP_EOL.PHP_EOL;
echo $Easting.PHP_EOL;
echo $Northing;

?>


Commentary

You'll need a php interpreter to run this, of course.  If you want to copy this stand-alone program from Google Drive then click on this.

The main routine does nothing but call the routine ToGreekGrid and print out the Easting and Northing results.  You'll modify that, of  course.  I see no barrier now to my adding the GGRS87 coordinates to every entry  in the Mycenaean Atlas Project.

If you try this then let me know what you think.

Yannis Lolos and the GGRS87

I was fortunate to receive Yannis Lolos’ Land of Sikyon a few days ago.  This is an investigation of the territory formerly dominated by Sicyon in the western Corinthia.  Not only is the book a valuable reference but it’s a joy just as a physical object – a beautifully made book.

One of its features is a gazetteer of all the find spots in the territory of Sicyon from the Neolithic to the Early Modern.  Each place is described and precisely located.  Only one problem.  Dr. Lolos’ does not give the coordinates in ordinary lat/lon form (technically the WGS84 projection) but in the Greek Grid system (GGRS87/EGSA87) which will be incomprehensible to most readers.  

What to do?

I found an easy-to-use site that will allow you to convert coordinate pairs from the GGRS87 projection to the WGS84 (what humans use).  It is here: 

http://www.themos-sa.gr/proj/index.php

To use this site you enter the GGRS coordinates in the upper right-hand box separated by commas.  You press the convert button and the lat/lon pair shows up in the upper left-hand box.  

So let's use a little Mycenaean settlement called 'Klimenti, Mastrogianni'  (Lolos [2011], p. 425, 'HS-11', 'Klimenti, Mastrogianni') as an example.  This is a small area (about 0.7 ha.) for which Dr. Lolos gives us the four corner points in GGRS 87 easting/northing pairs.  The first coordinate is 376799, 4203582.  You enter it into the conversion interface like this:

The conversion screen.  Enter the GGRS pair in the  right upper window (center),
press the 'Convert' button, and read the lat/lon (WGS84) pair in the upper left window.

The illustration shows you the conversion process for that point.  It maps to point 1 in this Google Earth image:

Four corner points for the site of Klimenti, Mastrogianni (HS-11).
Coordinates are numbered in the order in which Lolos presents them on p. 425.

Our first coordinate is the white arrow numbered '1' in the image above.   The set of four coordinates is this:

Easting	Northing	Lat	Lon
376799	4203582	37.97412198	22.59898503
376648	4203590	37.97417358	22.59726491
376653	4203631	37.97454367	22.59731479
376800	4203636	37.97460867	22.59898716

The 'file' button on this interface allows you to save a set of converted lat/lon pairs that have been converted to GGRS87 but I don't see that there's a way to create a file of GGRS87 points and have them converted all at once to WGS84.

I find an algorithm for converting from WGS84 to the Greek system at this web page:

http://studentguru.gr/b/kokyri/archive/2012/02/21/convert-coordinates-from-greek-grid-ggrs87-87-espg-2100-to-latitude-longitude-wsg84

This would have to be coded if I wanted, for example, to add the GGRS87 coordinate pairs to the 2000+ points in the Mycenaean Atlas Project, for example.  

I'll think about it.

Status of the Mycenaean Atlas Project:  www.helladic.info

I uploaded a new DB to helladic.info.  The specifications are: 

Database Update: Rev. 36. (MAP_Rev_0.036__09_13_17_Test)

Of the 1500 points in Dr. Murray's DB that I have targetted for initial integration I have completed about 1200.  There are currently 2055 separate sites in the Mycenaean Atlas Project DB.

Anyone who would like to have a copy of the MAP database can send an e-mail to bobconsoli 'at' gmail.com or leave a comment on any of my posts.  To run the MAP database requires a SQL server running on your desktop computer.   MySQL is such a server and it is powerful, industry-standard, and free.   Anyone who has a version of the DB from more than a month ago should request a new one.

I can and will make .kml or .kmz files, which can be opened directly in Google Earth, available to those who would like them.

  

I can also create .csv files for people who would like to import Mycenaean Atlas Project data into Google Earth but would like it in tabular form.  I do recommend that you use the website at www.helladic.info      Most of your needs can probably be met on that site and you can generate .kmz files for maps that you create.

If you like these posts then please follow me on Twitter (Squinchpix) or on Google+   (Robert Consoli)

Friends don't let friends use Facebook.

Bibliography

Lolos [2011]:      Lolos, Yannis A.   Land of Sikyon: Archaeology and History of a Greek City-State. Hesperia supplements, 39.   Princeton:  American School of Classical Studies at Athens, 2011.