4

Need to convert coordinates from lat/long to EPSG 3857 x/y

Example, using online converters

latlon = {33.971722106263044, -84.45613689691451};
epsgConverter[{lat_String, lon_String}, from_ : "4326", 
  to_ : "3857"] := {"x", "y"} /. 
  Import["https://epsg.io/trans?x=" <> lon <> "&y=" <> lat <> 
    "&s_srs=" <> from <> "&t_srs=" <> to]
epsgConverter[latlon]
(*{"-9401610.05", "4025002.66"}*)

What is the appropriate Wolfram Mathematica Function to provide this functionality? GeoPositionXYZ provides a different set of numbers. Tried all the datums with no avail.

GeoPositionXYZ[GeoPosition[latlon]][[1, 2 ;;]]
(*{-5.27024*10^6, 3.54385*10^6}*)
Zviovich
  • 9,308
  • 1
  • 30
  • 52
  • GeodesyData["Datum"] gives the possible formats {"BTS84", "BTS85", "BTS86", "BTS87", "ETRF00", "ETRF05", "ETRF14", "ETRF89", "ETRF90", "ETRF91", "ETRF92", "ETRF93", "ETRF94", "ETRF96", "ETRF97", "ETRS89", "EURM", "GDA94", "IGS2000", "ITRF0", "ITRF00", "ITRF05", "ITRF08", "ITRF14", "ITRF88", "ITRF89", "ITRF90", "ITRF91", "ITRF92", "ITRF93", "ITRF94", "ITRF96", "ITRF97", "NAD27", "NAD831986", "NAD83CORS93", "NAD83CORS94", "NAD83CORS96", "NAD83HARN", "OGB7", "PZ9002", "PZ90Original", "SIRGAS", "WGS72", "WGS84G730", "WGS84G873", "WGS84Original"}. – Roman Feb 11 '22 at 17:15
  • 2
    You can construct the transformation using GeodesyData to create a new projection, but I have found it much easier just to call epsg.io as you have here. – Carl Lange Feb 11 '22 at 17:58
  • (An entertaining, if stupid, way around this gap in WL is to use Predict on lists of coordinates and correctly projected variants and saving the PredictorFunction rather than painstakingly constructing the projection.) – Carl Lange Feb 11 '22 at 18:29
  • Tried that last night, none worked. GeoPositionXYZ[GeoPosition[latlon], #] & /@ GeodesyData["Datum"] – Zviovich Feb 11 '22 at 18:34
  • @CarlLange I was trying to figure out the transformation last night, and of course, did the same as u. call epsg.io and be done with it. The thing is that EPSG 3857 seems to be used extensively. I would have expected a simple way to be available in the WL – Zviovich Feb 11 '22 at 18:37
  • Sorry, I should have said GeoProjectionData rather than GeodesyData: GeoGridPosition[GeoPosition@latlon, {"Mercator", {"Centering" -> {0, 0}, "CentralScaleFactor" -> 1, "GridOrigin" -> {0, 0}, "ReferenceModel" -> "WGS84"}}] is very close – Carl Lange Feb 11 '22 at 18:47

1 Answers1

4

You can use GeoProjectionData for this:

GeoGridPosition[GeoPosition@latlon,
 {"Mercator", {"Centering" -> {0, 0}, "CentralScaleFactor" -> 1, 
   "GridOrigin" -> {0, 0}, "ReferenceModel" -> {6378137, 6378137}}}]

outputs:

GeoGridPosition[{-9.40161*10^6, 4.02501*10^6}, {"Mercator", 
  "ReferenceModel" -> {6378137, 6378137}}]

You can see more or less how I got there by reading the WKT definition of EPSG:3857:

PROJCS["WGS 84 / Pseudo-Mercator",
    GEOGCS["WGS 84",
        DATUM["WGS_1984",
            SPHEROID["WGS 84",6378137,298.257223563,
                AUTHORITY["EPSG","7030"]],
            AUTHORITY["EPSG","6326"]],
        PRIMEM["Greenwich",0,
            AUTHORITY["EPSG","8901"]],
        UNIT["degree",0.0174532925199433,
            AUTHORITY["EPSG","9122"]],
        AUTHORITY["EPSG","4326"]],
    PROJECTION["Mercator_1SP"],
    PARAMETER["central_meridian",0],
    PARAMETER["scale_factor",1],
    PARAMETER["false_easting",0],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]],
    AXIS["X",EAST],
    AXIS["Y",NORTH],
    EXTENSION["PROJ4","+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext  +no_defs"],
    AUTHORITY["EPSG","3857"]]
Carl Lange
  • 13,065
  • 1
  • 36
  • 70
  • Just curious about your deciphering: is there a more general way to import the EPSG WKT and convert that to a MMA GeoProjection? If such functionality is missing, what is the starting point to construct one? – sunt05 Apr 25 '22 at 13:49
  • Sadly, not yet. I think building a parser isn't necessarily that hard but it's quite easy to hit the limits of what MMA is willing to do without having to write a bunch of workarounds - the particular projection in this question was simple enough to not hit these limits, although even then it was a struggle. It would be great if WKT support were implemented fully in WL, for geometry and projections, but I have been holding my breath for quite a long time :) – Carl Lange Jun 28 '22 at 19:54