# How to ContourPlot a function of the coordinates on the Earth's surface on a map projection

There are two similar questions (this and this), but one is too specific and the other is apparently not clear and has no answers. I'll try to make the question more complete.

Suppose you have a vector field depending on the coordinates x,y,z of a point on the Earth's surface [Edit: I found a mistake on this expression but to keep discussions consistent I'll leave as is.]

field = With[{rl$= Transpose[{{1.45031*10^7, -1.46446*10^11}, {-3.86086*10^8, 2.25113*10^8}, {1.00319*10^8, 3.47297*10^10}} - 1. {##1}]}, Total[({4.90254*10^12, 1.32705*10^20} rl$)/rl$^3] - 1. {-0.00569904, -0.0000209977, 0.00135953}] &  With some trick I can represent the Norm of this vector field on a map with the equirectangular projection. If I well understood the topic: [Edit: I found another supposed mistake in the original code following. From what I now understand on a map usaully we represent the geodetic latitude. But the well known parametrization of an oblate spheroid (the one I originally used) refers to the geocentric latitude. The difference between the two expression is not more than about 0.2° for GRS80, but to be rigorous we need to take into account this difference. This is the reason for the function toGeocentric I added. Feel free to comment on this to confirm or to contradict my understanding of this topic.] {a, b} = GeodesyData["GRS80", #] & /@ {"SemimajorAxis", "SemiminorAxis"} // UnitConvert // QuantityMagnitude toGeocentric = λ \[Function] Evaluate[GeodesyData[ "GRS80", {"GeocentricLatitude", λ}] /. λ -> λ °] gg = GeoGraphics[GeoRange -> "World", GeoProjection -> "Equirectangular", GeoGridLines -> Quantity[15, "AngularDegrees"], ImageSize -> 800] cp = ContourPlot[ Norm[field[a Cos[toGeocentric @θ °] Cos[φ °], a Cos[toGeocentric @θ °] Sin[φ °], b Sin[toGeocentric @θ °]]], {φ, -180, 180}, {θ, -90, 90}, ColorFunction -> (Opacity[0.6, ColorData["TemperatureMap"][#]] &)] Show[Graphics @@ gg, cp]  It is the "right" way to do that (I'm a bit concerned by Graphics@@gg for example, I read that on another question)? It is possible to use many projections (ideally all map projections supported by Mathematica) in a consistent and general way without having to know and to input all the necessary projection transformations but taking advantage of the cartography framework in Mathematica? I'm interested for example to "LambertAzimuthal", "Bonne", "Robinson"... • "different projections" - so as long as you know the explicit expressions for$\theta$and$\varphi$in terms of the other projection's coordinates, you just plug things in, yes? – J. M.'s ennui May 9 '15 at 22:07 • Yes, I suppose. But, If possible, I prefer a way that can handle ideally all (more common, world's) map projections supported by Mathematica without even having to know and to input all the details of these projection transformations. – unlikely May 9 '15 at 22:19 ## 2 Answers You can use your ContourPlot, you just need to wrap coordinates with GeoPosition, note that you have to flip order. (if x is longitude and y is latitude, because GeoPosition assumes first is latitude and so on.) cp is Graphics[GraphicsComplex[coordinates, primitives]...], it is convenient to use this form. We can apply GeoPosition in one place and reduce amout of calls to Wolfram server: MapAt[GeoPosition[#[[;; , {2, 1}]]] &, cp[[1]], {1}]  #[[;;,{2,1}]] means from {{x1, y1}, {x2, y2}...} we get {{y1, x1}, {y2, x2}... So: GeoGraphics[MapAt[GeoPosition[#[[;; , {2, 1}]]] &, cp[[1]], {1}], GeoRange -> "World", GeoProjection -> "Robinson", GeoGridLines -> Quantity[15, "AngularDegrees"], ImageSize -> 800]  And with those GeoPositions GeoGraphics will convert everything automatically. I think this answer of mine is related because I'm using projections without knowing their forms: Mollweide maps in Mathematica In general GeoGraphics are a little crippled because they are not handling plot range, cutting primitives on edge and similar cases. There is also a bug which makes that even harder: Compability between Graphics coordinates and shifted projections in GeoGraphics So very often you will have to handle coordinates carefully as it won't work automatically for every projection: It works for most projections, old bug was fixed but there are new issues I need to investigate for problematic cases: {a, b} = GeodesyData["GRS80", #] & /@ {"SemimajorAxis", "SemiminorAxis"} // UnitConvert // QuantityMagnitude cp = ContourPlot[ Norm[field[a Cos[θ Degree] Cos[φ Degree], a Cos[θ Degree] Sin[φ Degree], b Sin[θ Degree]]], {φ, -180, 180}, {θ, -90, 90}, ColorFunction -> (Opacity[0.6, Blend["TemperatureMap", #]] &), ContourStyle -> Directive[Thick, Dashed], Contours -> 4]; data = MapAt[ GeoPosition[#[[A ll, {2, 1}]]] &, Cases[cp, _GraphicsComplex, Infinity][[1]], {1} ] GeoGraphics[data, GeoRange -> "World", GeoProjection -> #, GeoGridLines -> Quantity[15, "AngularDegrees"], ImageSize -> {300, 300}, PlotLabel -> #, LabelStyle -> 18 ] & /@ GeoProjectionData[][[;; 3]]  • Thanks. This is the approach I was searching. Two possible issue: the first is that many projections are supported but some aren't. The second is: we use ContourPlot to adaptively sample the function on an Equirectangular projection but later we use these samples on another projection plane: this can eventually produce non-smooth contours with more complex functions? I understand you already did much for me, so if you can just give me a hint... – unlikely May 10 '15 at 18:58 • He's saying that the adaptive algorithm of ContourPlot[] gave a result adapted to equirectangular coordinates; if you then transform those sample points to, say, sinusoidal coordinates, the result might look rougher than it would have been if the algorithm adapted to the new coordinates. – J. M.'s ennui May 10 '15 at 19:48 • @J. M. I see, well, I don't know any good solution if we don't have those transformations. – Kuba May 10 '15 at 19:55 • One could try, to set MaxRecursion->1 and add many PlotPoints to not waste time on features that may not be relevant in target coordinates – Kuba May 10 '15 at 19:58 • @J. M. Thanks for translating – unlikely May 10 '15 at 21:43 Trying to handle the two concerns expressed in a comment to the Kuba answer, I built my own method. Indeed, I have investigated and tested 6 methods but the one I consider a good compromise between supported projections, (supposed) accuracy and speed is the following. Options[GeoContourPlot] = Join[ Options[GeoGraphics], FilterRules[Options[ContourPlot], Complement @@ (Map[First]@*Options /@ {ContourPlot, GeoGraphics})] ]; GeoContourPlot[fun_, opts : OptionsPattern[]] := Module[{gg, gp, geoModel, geoProjection, bmesh, mesh, vlp, vlc, fl, cp}, gg = GeoGraphics[FilterRules[{opts}, Options[GeoGraphics]]]; {geoModel, geoProjection} = {GeoModel, GeoProjection} /. AbsoluteOptions[gg, {GeoModel, GeoProjection}]; gp = FirstCase[gg, Polygon[vl_, ___, VertexTextureCoordinates -> _, ___] :> Polygon[DeveloperToPackedArray@vl],$Failed, \[Infinity]];
(* Area[gp] sometimes never ends: Area[bmesh] is immediate *)
bmesh = BoundaryDiscretizeRegion@gp;
mesh = TriangulateMesh[bmesh,
MaxCellMeasure -> Area@bmesh/Times @@ (OptionValue[PlotPoints] /.
{Automatic -> 250} /. {n :
Except[_List] :> {n}} /. {{n_} :> {n, n}})
];
vlp = GeoGridPosition[MeshCoordinates[mesh], geoProjection, geoModel];
vlc = Quiet[
GeoPositionXYZ[vlp], {CompiledFunction::cfne, Power::infy,
Infinity::indet, GeoPosition::ltrng}];
If[! MatrixQ[vlc[[1]], NumericQ],
(* For few projections, the way I used to get the projected domain is unreliable,
not all projected nodes can be unprojected, so it's necessary to process one at time (slow!) *)
vlc = Quiet[
Power::infy, Infinity::indet, GeoPosition::ltrng}];
With[{sel = FreeQ[#, GeoGridPosition] & /@ vlc},
vlp = MapAt[Pick[#, sel] &, vlp, {1}];
vlc = GeoPositionXYZ[
Pick[vlc, sel][[All, 1]] // DeveloperToPackedArray]
];
];
fl = fun @@@ vlc[[1]];
Evaluate@FilterRules[{opts}, Options[ListContourPlot]]];
GeoGraphics[cp[[1]], FilterRules[{opts}, Options[GeoGraphics]]]
]


The results are promising:

gl = Table[
GeoContourPlot[Norm@*field,
GeoRange -> All, GeoProjection -> geoProjection,
GeoGridLines -> Quantity[15, "AngularDegrees"],
ImageSize -> {200, 200}, PlotPoints -> 50, PlotRangePadding -> .2,
PlotLabel -> geoProjection,
ColorFunction -> (Opacity[0.6, ColorData["TemperatureMap"][#]] &),
ContourStyle -> Thin, Contours -> 8        ],
{geoProjection,
DeleteCases[GeoProjectionData[],
s_String /;
StringMatchQ[s,
StartOfString ~~ ("SPCS" | "UTMZone" | "UPS") ~~ __]]}
];
Grid[Partition[gl, 5, 5, {1, -1}, ""], Dividers -> All]


Some projection is still unsupported and I made choices that deserve furter investigations, but for me, I think these results are enough.

If someone is interested, in a future update to this self-answer, I can add some details on the other methods and on some choice made for this method.

• Nice results :) – Kuba May 12 '15 at 6:03
• @Kuba Thanks, but still unsatisfied by the unreliable way used to get the domain on the projected plane... – unlikely May 12 '15 at 8:01
• Polygons should be automatically cut on edge or antimeridian, this functionality is just Beta version! ps. take a look at Mercator and MillerCylindrical, don't you think that something is odd in area of inner disks? It looks like sometimes area of SouthAmerica taken is varying. – Kuba May 12 '15 at 8:09
• This is astonishing example of d3.geo library: jasondavies.com/maps/transition I tired to do the same, it is notpossible in MMA unless you fix GeoGraphics for Wolfram. So I decided I have no time for that :/ – Kuba May 12 '15 at 8:14