How can I output the coordinates of a city in coordinates of the Mercator projection?

Question

On Mathematica version 8, I would like to use CountryData and CityData to mark Paris on a polygon map of France in the Mercator projection. In the example below, you will see that I'm able to mark Paris on a default/equirectangular projection map of France, but not on a Mercator projection.

First, I'll use CountryData to obtain coordinate data for France:

franceCoords1 = CountryData["France", "Polygon"][[1]];
franceCoords2 = CountryData["France", {"Polygon", "Equirectangular"}][[1]];
franceCoords3 = Map[Reverse, CountryData["France", "Coordinates"], {2}];
franceCoords4 = CountryData["France", {"Polygon", "Mercator"}][[1]];

The first three lists are actually identical:

(franceCoords1 == franceCoords2) && (franceCoords2 == franceCoords3)
franceCoords1 == franceCoords4

True

False

The first three are the default/equirectangular projection, while franceCoords4 is the Mercator projection.

I can plot all these coordinates by using Graphics wrapped around Polygon objects:

Grid[{
  Map[
   Graphics[Polygon[#], Frame -> True, ImageSize -> 200] &
   , {franceCoords1, franceCoords2, franceCoords3, franceCoords4}]
  }]

It's pretty clear that the first three are the same (default/equirectangular) projection, whereas franceCoords4 is a different (i.e., Mercator) projection.

Now, I can obtain the coordinates of Paris by using CityData:

parisCoords = Reverse[CityData[{"Paris", "France"}, "Coordinates"]]

{2.34, 48.86}

(CityData gives coordinates as {latitude, longitude} ordered pairs, whereas a "conventional" world map is {x, y} = {longitude, latitude}, so in parisCoords I have reversed the coordinates given by CityData.)

I can mark Paris on the maps by using Point[parisCoords]:

Grid[{
  Map[
   Graphics[{
      Polygon[#],
      Red, PointSize[Large], Point[parisCoords]
      }, Frame -> True, ImageSize -> 200] &
   , {franceCoords1, franceCoords2, franceCoords3, franceCoords4}]
  }]

It looks like the location of Paris is correct in franceCoords1, franceCoords2, and franceCoords3 ... but not in franceCoords4. This is because the maps given by franceCoords1, franceCoords2, and franceCoords3 are given in default/equirectangular coordinates, whereas the map given by franceCoords4 is given by Mercator projection coordinates.

So, my question is, how can I convert the default/equirectangular coordinates of Paris to coordinates in the Mercator projection?

Probably doing so would be simpler in Mathematica version 10, which I understand introduced new geography features, but I don't have access to version 10. I would like to accomplish it in version 8, please.

Is there any built-in way for CityData to output the coordinates of Paris in the Mercator projection, or to convert the default/equirectangular coordinates to coordinates in the Mercator projection?

Answer 1

Ok, I had to tidy this answer up:

Functional way to change coordinates in Mathematica is to use combination of GeoPosition and GeoGridPosition:

GeoGridPosition[GeoPosition[{lat, long}], (*projection*) ]

---

The answer should be done here but there is some unconsistency for older than V10 versions.

There is different scalling applied to directly downloaded data and different when you manually convert from equirectangular coordinates:

    (*direct download world in mercator projection*)
worldMercV9 = CountryData["World", {"SchematicPolygon", "Mercator"}];

    (*manual conversion*)
worldV9 = CountryData["World", "SchematicPolygon"];
worldFromEqToMercV9 = Map[
   First @ GeoGridPosition[GeoPosition[Reverse[#]], "Mercator"] &, 
   worldV9, 
   {3}]



Graphics[#, Frame -> True] & /@ {worldMercV9, worldFromEqToMercV9} // Row

so yeah, good luck.

• You can't directly download Paris coordinates in mercator projection so you have to take care about rescalling by yourself.

• Or you can go other way. Download everything in equirectangular projection and then proceed with changing coordinates, as kguler has demonstrated.

---

franceMERC = CountryData["France", {"Polygon", "Mercator"}];

parisEq = CityData["Paris", "Coordinates"] (* => {48.86, 2.34}*)

parisMerc = Composition[
             Rescale[# - {Pi, 0}, {-Pi, Pi}, {-180., 180}] &,
             First,
             GeoGridPosition[#, "Mercator"] &,
             GeoPosition
            ][parisEq]    (* => {2.34, 56.1549}*)

Graphics[{
  franceMERC,
  AbsolutePointSize@12, Red, Point@parisMerc
  },
  Frame -> True]

---

Fortunately, in V10 it is consistent.

Answer 2

Perhaps

mpF = GeoGridPosition[GeoPosition[#], "Mercator"][[1]] &;

g1 = Graphics[{Gray, CountryData["France", "Polygon"], 
               PointSize[Large], Red, Point[parisCoords]}, ImageSize -> 300];
g2 = Graphics[{Gray, Polygon[mpF /@ CountryData["France", "Coordinates"]], Red, 
               PointSize[Large], Point@mpF[Reverse@parisCoords]}, ImageSize -> 300];

Row[{g1, g2}]