2

My coordinates are projected using the following projection:

proj= {"UTMZone32", {"GridOrigin" -> {500000, 0}, "CentralScaleFactor" -> 0.9996}};

Now I wish to calculate the distance between two points (ignoring elevation), e.g.

p1= GeoGridPosition[{359577, 5.51291*10^6,0}, proj]
p2= GeoGridPosition[{509108, 5.972*10^6,0}, proj]

When I try GeoDistance

GeoDistance[p1,p2]

it fails with the error message

GeoDistance::invparam: "Invalid parameters \!\(\"GeoGridPosition[{359577, 5.51291*^6, 0},
{\\\"UTMZone32\\\", {\\\"GridOrigin\\\" -> {500000, 0}, 
 \\\"CentralScaleFactor\\\" -> 0.9996}}]\"\). "

Also, the GeoPositionXYZ function, as in

GeoPositionXYZ[p1]

fails with the error messages

ToString::nonopt: Options expected (instead of InputForm) beyond position 2 in 
ToString[None,{GridOrigin->{500000,0},CentralScaleFactor->0.9996},InputForm]. 
An option must be a rule or a list of rules. >>

GeoGridPosition::invparam: "Invalid parameters ToString[!(None, {
 "GridOrigin" -> {500000, 0}, "CentralScaleFactor" -> 0.9996`}, InputForm)]."


GeoPositionXYZ::invcoord: "!("GeoPosition[GeoGridPosition[{359577, 5.51291*^6, 0}, 
{\"UTMZone32\", {\"GridOrigin\" -> {500000, 0}, \"CentralScaleFactor\" -> 
0.9996}}]]") is not a valid coordinate specification."

Both functions work, however, when I switch proj to the string UTMZone32.

Do I need to get the full projection specification to work?

EDIT: After some further googling, I realized that in UTM coordinates the distance between two points is simply 

Norm[{p1[[1,1;;2]]-p2[[1,1;;2]]}]

so I would answer my own question with no.

Kuba
  • 136,707
  • 13
  • 279
  • 740
Karsten W.
  • 1,383
  • 8
  • 21
  • 1
    I'm not sure about this, but I think that you can't use {"UTMZone32", {"GridOrigin" -> {500000, 0}, "CentralScaleFactor" -> 0.9996}} as a projection. "UTMZone32" is a defined projection on its own: GeoProjectionData["UTMZone32"] ==> {"TransverseMercator", {"Centering" -> {0, 9}, "CentralScaleFactor" -> 1, "GridOrigin" -> {0, 0}, "ReferenceModel" -> "WGS84"}} Given the centering and scaling you want perhaps you could use proj = {"TransverseMercator", {"GridOrigin" -> {500000, 0}, "CentralScaleFactor" -> 0.9996, "Centering" -> {0, 9}, "ReferenceModel" -> "WGS84"}} – Sjoerd C. de Vries Jan 23 '14 at 20:06

1 Answers1

3

As @Sjoerd states in the comments, your projection system (UTMZone32) has a defined set of parameters. You can check these using GeoProjectionData:

GeoProjectionData["UTMZone32"]

{"TransverseMercator", {"Centering" -> {0, 9}, "CentralScaleFactor" -> 0.9996, "GridOrigin" -> {500000, 0}, "ReferenceModel" -> "WGS84"}}

These coincide with the ones you are trying to set.

To define your own projection system similar to UTM (based on Transverse Mercator), you can simply specify those in GeoGridPosition:

GeoGridPosition[{1000000, 1000000}, 
 {"TransverseMercator", {"Centering" -> {0, 0}, "CentralScaleFactor" -> 0.95,
 "GridOrigin" -> {500000, 0}, "ReferenceModel" -> "WGS84"}}]

This now can be easily converted to LatitudeLongitude.

So, since this is a projected coordinate system, and as you state at the end of the question, can be easily calculated using Norm or EuclideanDistance or whatever:

Norm[{359577, 5.51291*10^6, 0} - {509108, 5.972*10^6, 0}]
EuclideanDistance[{359577, 5.51291*10^6, 0}, {509108, 5.972*10^6, 0}]

482828.

482828.

But we can also use the built-in GeoDistance which in v10 returns a Quantity:

pos1 = GeoGridPosition[{359577, 5.51291*10^6, 0}, "UTMZone32"];
pos2 = GeoGridPosition[{509108, 5.972*10^6, 0}, "UTMZone32"];
GeoDistance[pos1, pos2]~UnitConvert~"Meters"

482985. m

Sadly, they're 157 meters apart.

kale
  • 10,922
  • 1
  • 32
  • 69
  • These two locations are over 480 km apart and the farther apart two points are in any map projection, the larger the difference between the Euclidean distance and the geodesic distance given by GeoDistance. If the two locations were closer together, the difference in the distances would be much smaller. So I am quibbling with the use of the word "Sadly". That's simply what happens with all 3D-to-2D map projections. – JimB Aug 02 '16 at 14:28
  • 1
    @JimBaldwin Haha. Exactly... Sadly errors arise when you try to view a sphere in 2D... Not sadly Wolfram's algorithm are suspect... – kale Aug 03 '16 at 23:23