Convert ETRS89_Portugal_TM06 coordinate to Google Map Coordinate

Question

I am trying convert an coordinate (point = {-26841,7714704367;42953,3965652831}) from shapefile with follow prj file:

PROJCS["ETRS89_Portugal_TM06",GEOGCS["GCS_ETRS_1989",DATUM["D_ETRS_1989",SPHEROID["GRS_1980",6378137,298.257222101]],PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",39.66825833333333],PARAMETER["central_meridian",-8.133108333333334],PARAMETER["scale_factor",1],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["Meter",1]]

To a Google Maps coordenate. I am reading about it and I found that I need to transform this coordenate into EPSG 3857. I am using the follow code to do that:

const string latlongwkt = "PROJCS[\"Popular Visualisation CRS / Mercator\", GEOGCS[\"Popular Visualisation CRS\", DATUM[\"WGS84\", SPHEROID[\"WGS84\", 6378137.0, 298.257223563, AUTHORITY[\"EPSG\",\"7059\"]], AUTHORITY[\"EPSG\",\"6055\"]], PRIMEM[\"Greenwich\", 0, AUTHORITY[\"EPSG\", \"8901\"]], UNIT[\"degree\", 0.0174532925199433, AUTHORITY[\"EPSG\", \"9102\"]], AXIS[\"E\", EAST], AXIS[\"N\", NORTH], AUTHORITY[\"EPSG\",\"4055\"]], PROJECTION[\"Mercator\"], PARAMETER[\"semi_minor\",6378137], PARAMETER[\"False_Easting\", 0], PARAMETER[\"False_Northing\", 0], PARAMETER[\"Central_Meridian\", 0], PARAMETER[\"Latitude_of_origin\", 0], UNIT[\"metre\", 1, AUTHORITY[\"EPSG\", \"9001\"]], AXIS[\"East\", EAST], AXIS[\"North\", NORTH], AUTHORITY[\"EPSG\",\"3785\"]]";
IProjectedCoordinateSystem googleProjectedCoordinates = CoordinateSystemWktReader.Parse(latlongwkt) as IProjectedCoordinateSystem;

var wktstring = "PROJCS[\"ETRS89_Portugal_TM06\",GEOGCS[\"GCS_ETRS_1989\",DATUM[\"D_ETRS_1989\",SPHEROID[\"GRS_1980\",6378137,298.257222101]],PRIMEM[\"Greenwich\",0],UNIT[\"Degree\",0.017453292519943295]],PROJECTION[\"Transverse_Mercator\"],PARAMETER[\"latitude_of_origin\",39.66825833333333],PARAMETER[\"central_meridian\",-8.133108333333334],PARAMETER[\"scale_factor\",1],PARAMETER[\"false_easting\",0],PARAMETER[\"false_northing\",0],UNIT[\"Meter\",1]]";
IProjectedCoordinateSystem sourceCS = CoordinateSystemWktReader.Parse(wktstring) as IProjectedCoordinateSystem;

CoordinateTransformationFactory ctFactory = new CoordinateTransformationFactory();
ICoordinateTransformation transformer1 = ctFactory.CreateFromCoordinateSystems(sourceCS, googleProjectedCoordinates);

double[] toPoint = transformer1.MathTransform.Transform(fromPoint);

But it still give me a wrong coordinate.

I also tried to replace IProjectedCoordinateSystem googleProjectedCoordinates to:ICoordinateSystem geographicCoordinates = ProjNet.CoordinateSystems.GeographicCoordinateSystem.WGS84

And I got and lat/log coordinate but with wrong location in Google Map.

I already spent a lot of time and I'm stuck on it. Could anyone give me some help? I really can't figure out what is wrong and why it didn't work.

I am using the Proj.Net (http://projnet.codeplex.com/) to apply the projection.

I am using the wkt (that I think is compatible with Google Maps coordinates, correct me if I am wrong) provided in this article (https://alastaira.wordpress.com/2011/01/23/the-google-maps-bing-maps-spherical-mercator-projection/)

---

I made some changes in code, I changed proj.net to DotSpatial that it is based in Proj.4 and it was recently updated (Thu Apr 28, 2016 at 8:00 AM).

So the code is the follow:

Shapefile nycBoroughs = Shapefile.OpenFile(filename); nycBoroughs.Reproject(KnownCoordinateSystems.Geographic.Europe.ETRS1989);

and now I got the coordinates as you said: coordenate = {(-8,44768822459892, 40,0546873641873)}

I also apply the formulas and I got the coordinates also as you said: x = -940392.35141178139, y = 4873892.4828391289

But now I can't use this coordinates directly in Google Maps, I need to convert again to ETRS1989 to get the coordinates in degree and use it with google maps?

Sorry, probably it is a simple question, but I am a programmer and i don't have enough know-how in system coordinates to understand it easily.

Answer 1

The Proj.Net homepage notes that the code is rather old, from 2009. The development of Proj.4 has improved since then, including proper use of the Google Mercator projection.

If you can't use the Proj.4/GDAL libraries, I suggest to reproject from TM06 to ETRS89/WGS84 coordinates (should be around (-8.4477 40.0547)),then convert to Google mercator using the formulas in the linked article:

double x = lon * 20037508.34 / 180;
double y = Math.Log(Math.Tan((90 + lat) * Math.PI / 360)) / (Math.PI / 180) * 20037508.34 / 180; 

Note that Math.Log is the logarithmus naturalis, not decimalis.

The result should be (-940392 4873892).

Answer 2

Sorry I was lazy and I didn't read the end of the document. The code that I need to write to get the correct conversion is (could be useful in future for someone):

Shapefile nycBoroughs = Shapefile.OpenFile(filename);
nycBoroughs.Reproject(KnownCoordinateSystems.Geographic.Europe.ETRS1989);

foreach (IFeature feature in nycBoroughs.Features) {
  DataRow attributes = feature.DataRow;

    foreach (Coordinate coordenate in feature.Coordinates) {
      double x = coordenate.X * 20037508.34 / 180;
      double y = (Math.Log(Math.Tan((90 + coordenate.Y) * Math.PI / 360)) / (Math.PI / 180)) * 20037508.34 / 180;

      double lon = (x / 20037508.34) * 180;
      double lat = (y / 20037508.34) * 180;

      lat = 180/Math.PI * (2 * Math.Atan(Math.Exp(lat * Math.PI / 180)) - Math.PI / 2);

      //Convert Points
      locationCollection.Add(new Microsoft.Maps.MapControl.WPF.Location(lat, lon));
...

I hope it is useful in future for someone.

Thank you for your help AndreJ. I wouldn't be able to do this without your help.

Answer 3

If you are using Javascript use this library https://github.com/proj4js/proj4js.

proj4.defs('ETRS89_Portugal_TM06', `PROJCS["ETRS89_Portugal_TM06",GEOGCS["GCS_ETRS_1989",DATUM["D_ETRS_1989",SPHEROID["GRS_1980",6378137,298.257222101]],PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",39.66825833333333],PARAMETER["central_meridian",-8.133108333333334],PARAMETER["scale_factor",1],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["Meter",1]]`)

proj4('ETRS89_Portugal_TM06',[-71,41]);

If you are extracting the zip file directly from official sources you will have a .prj with the projection text you should use.