2

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.

6
  • 1
    In your code, you used EPSG:3785 instead of EPSG:3875. These deliver different coordinates. Please add which software/libs you are using.
    – AndreJ
    Commented May 29, 2017 at 5:25
  • Hi AndreJ, First of all i want to thank you for your reply. About the libs, I am using the Proj.Net to make the projections (I added it to the post). About the wkt code, I am very confuse about it. Could you tell me what wkt code should I use to convert the coordinates to lat/long compatible with Google Maps? I am reading about it to help a friend with programming code, but I am stuck and I am not able to find the projection to convert the coordinates in lat/long compatible with Google Maps.
    – André
    Commented May 29, 2017 at 22:29
  • Hi AndreJ, to make sure that the projection is ok, I made an test using mygeodata.cloud/cs2cs and I used the follow coordinates to test the projection: point = {-940392,351542685;4873892,4833824} (in ETRS89_Portugal_TM06 - EPSG:3763) I Select the option to convert to (ETRS89 - ETRS89 (EPSG:4258))
    – André
    Commented May 30, 2017 at 18:29
  • I trust cs2cs more than Proj.Net. See my answer.
    – AndreJ
    Commented May 30, 2017 at 18:34
  • Hi AndreJ, to make sure that the projection is ok, I made a test using the website mygeodata.cloud/cs2cs and I used the follow coordinates to test the projection: point = {-940392,351542685;4873892,4833824} (in ETRS89_Portugal_TM06 - EPSG:3763). I Select the option to convert to (ETRS89 (EPSG:4258)), and the result is: -35.038184551;80.3124394135 that is wrong. The EPSG:4258 is the system coordinate that you are talking? I checked the website epsg.io/4258 and it seems correct. Sorry i am not able to edit the previous comment.
    – André
    Commented May 30, 2017 at 18:36

3 Answers 3

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).

1
  • can you include the code to convert from TM06 to ETRS89/WGS84? I really need the full function cause this is a bit new for me Commented Apr 24, 2021 at 21:26
1

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.

1

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.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.