11
$\begingroup$

How to use the UTM coordinate system (northing, easting) in Mathematica using the built-in functions (GeoPositionENU, GeoGridPosition, etc.). There are UTM Zones implemented in Mathematica, e.g.:

GeoProjectionData["UTMZone33"]

gives

{"TransverseMercator", {"Centering" -> {0, 15},   "CentralScaleFactor" -> 1, "GridOrigin" -> {0, 0},   "ReferenceModel" -> "WGS84"}}

From N45, E15 one should get Easting 500000 and Northing 4982950.4 (zone 33), but I can't reproduce it. 

GeoGridPosition[GeoPosition[{45, 15, 0}, "WGS84"], "UTMZone33"]

gives

GeoGridPosition[{0., 4.98494*10^6, 0}, "UTMZone33"]

(Easting offset 500000 is understandable - default, but Northing is way off.)

Any help would be appreciated.

Improve this question
$\endgroup$
4
  • $\begingroup$ I've wanted to know this for a while too. Right now, I just use the function included in this blog post to go from LatLon to UTM: blog.wolfram.com/2009/04/17/mapping-gps-data $\endgroup$
    – kale
    Commented Jan 8, 2013 at 14:56
  • 2
    $\begingroup$ The output of GeoProjectionData is not correct for the UTM system: the value of GridOrigin is wrong (it should be {500000,0} and so is the value of CentralScaleFactor (it should be 0.9996). $\endgroup$
    – whuber
    Commented Jan 8, 2013 at 17:04
  • $\begingroup$ @kale: the short version of the function would therefore be: ll2utm[coord_] := GeoGridPosition[ GeoPosition[coord, "WGS84Original"], {"UTMZone33", "CentralScaleFactor" -> 0.9996, "GridOrigin" -> {500000, 0}}] [[1, {1, 2}]] and then ll2utm[{45, 15}] // Round gives {500000, 4982950}. $\endgroup$
    – Boocko
    Commented Jan 9, 2013 at 11:08
  • $\begingroup$ @Boocko, Except, you would need to include the logic to figure out what zone you were in and build a string "UTMZone"<>ToString@zone. $\endgroup$
    – kale
    Commented Jan 9, 2013 at 15:19

2 Answers 2

Reset to default
6
$\begingroup$

As in comments by WReach, the correct answer is

GeoGridPosition[ 
  GeoPosition[{45, 15, 0}, "WGS84"], 
  {"UTMZone33", "CentralScaleFactor" -> 0.9996, "GridOrigin" -> {500000, 0}}
]

The problem lies with the scale factor: 0.9996 and the grid origin

If you are familiar with projection systems used in GIS systems, you can check out http://spatialreference.org/ref/epsg/32633/

Details are here:

PROJCS["WGS 84 / UTM zone 33N",
    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.01745329251994328,
            AUTHORITY["EPSG","9122"]],
        AUTHORITY["EPSG","4326"]],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]],
    PROJECTION["Transverse_Mercator"],
    PARAMETER["latitude_of_origin",0],
    PARAMETER["central_meridian",15],
    PARAMETER["scale_factor",0.9996],
    PARAMETER["false_easting",500000],
    PARAMETER["false_northing",0],
    AUTHORITY["EPSG","32633"],
    AXIS["Easting",EAST],
    AXIS["Northing",NORTH]

More explanations of the UTM system can be found here: "The central meridian in each UTM zone has a scale factor of 0.9996"

Improve this answer
$\endgroup$
5
  • 2
    $\begingroup$ I feel this falls a bit short as an answer. Could you elaborate somewhat? $\endgroup$
    – Sjoerd C. de Vries
    Commented Jan 8, 2013 at 16:29
  • 1
    $\begingroup$ I agree with your diagnosis, although the easting should also be 500000 as the OP observes. You might want to mention that the desired result can be obtained using the expression GeoGridPosition[ GeoPosition[{45, 15, 0}, "WGS84"], {"UTMZone33", "CentralScaleFactor" -> 0.9996, "GridOrigin" -> {500000, 0}}]. $\endgroup$
    – WReach
    Commented Jan 8, 2013 at 19:18
  • $\begingroup$ Thank you, this works. But to convert back, the simple GeoPosition on previous output should work, but I get errors. Looks like the GeoPosition is not happy with {"UTMZone33",...}. $\endgroup$
    – Boocko
    Commented Jan 9, 2013 at 10:51
  • 1
    $\begingroup$ works in this way: GeoProjectionData["UTMZone33"] gives {"TransverseMercator", {"Centering" -> {0, 15}, "CentralScaleFactor" -> 1, "GridOrigin" -> {0, 0}, "ReferenceModel" -> "WGS84"}}, then use this: GeoPosition@ GeoGridPosition[ GeoPosition[{45, 15, 0}, "WGS84"], {"TransverseMercator", "Centering" -> {0, 15}, "CentralScaleFactor" -> 0.9996, "GridOrigin" -> {500000, 0}, "ReferenceModel" -> "WGS84"}] $\endgroup$
    – Tuku
    Commented Jan 9, 2013 at 15:30
  • $\begingroup$ @Tuku that works flawlessly. Thanks. $\endgroup$
    – Boocko
    Commented Jan 10, 2013 at 10:27
0
$\begingroup$

The settings appear to have been "fixed" (at least for Mathematica version 12).

Using Los Angeles, CA as an example:

(* Convert from LatLong to UTM *)
(* You need to know what UTM Zone you want *)
lat = 34.052234;
long = -118.243685;
elev = 86.654;
utmZone = "UTMZone11";
utm = GeoGridPosition[GeoPosition[{lat, long, elev}], utmZone]
(* GeoGridPosition[{385215.34335191007`,3.7686452431931966`*^6,86.654`},"UTMZone11"] *)

(* Convert from UTM to LatLong *)
utmZone = "UTMZone11";
utmEasting = 385215.34335191007;
utmNorthing = 3768645.2431931966;
elev = 86.654;
GeoPosition[GeoGridPosition[{utmEasting, utmNorthing, elev}, utmZone]]
(* GeoPosition[{34.052234000000006`,-118.243685`,86.654`}] *)
Improve this answer
$\endgroup$

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.