17
$\begingroup$

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]

Mathematica graphics

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"...

Improve this question
$\endgroup$
2
  • $\begingroup$ "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? $\endgroup$
    – J. M.'s lack of A.I. ♦
    May 9, 2015 at 22:07
  • $\begingroup$ 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. $\endgroup$
    – unlikely
    May 9, 2015 at 22:19

2 Answers 2

Reset to default
23
$\begingroup$

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.

enter image description here

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]]

enter image description here

Improve this answer
$\endgroup$
5
  • $\begingroup$ 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... $\endgroup$
    – unlikely
    May 10, 2015 at 18:58
  • 1
    $\begingroup$ 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. $\endgroup$
    – J. M.'s lack of A.I. ♦
    May 10, 2015 at 19:48
  • $\begingroup$ @J. M. I see, well, I don't know any good solution if we don't have those transformations. $\endgroup$
    – Kuba ♦
    May 10, 2015 at 19:55
  • 1
    $\begingroup$ 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 $\endgroup$
    – Kuba ♦
    May 10, 2015 at 19:58
  • $\begingroup$ @J. M. Thanks for translating $\endgroup$
    – unlikely
    May 10, 2015 at 21:43
7
$\begingroup$

Animation

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[Developer`ToPackedArray@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[
    GeoPositionXYZ /@ Thread@vlp, {CompiledFunction::cfne, 
     Power::infy, Infinity::indet, GeoPosition::ltrng}];
  With[{sel = FreeQ[#, GeoGridPosition] & /@ vlc},
   vlp = MapAt[Pick[#, sel] &, vlp, {1}];
   vlc = GeoPositionXYZ[
     Pick[vlc, sel][[All, 1]] // Developer`ToPackedArray]
   ];
  ];
 fl = fun @@@ vlc[[1]];
 cp = ListContourPlot[MapThread[Append, {vlp[[1]], fl}], 
   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"], 
    GeoRangePadding -> [email protected],
    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]

Image

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.

Improve this answer
$\endgroup$
4
  • $\begingroup$ Nice results :) $\endgroup$
    – Kuba ♦
    May 12, 2015 at 6:03
  • $\begingroup$ @Kuba Thanks, but still unsatisfied by the unreliable way used to get the domain on the projected plane... $\endgroup$
    – unlikely
    May 12, 2015 at 8:01
  • $\begingroup$ 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. $\endgroup$
    – Kuba ♦
    May 12, 2015 at 8:09
  • $\begingroup$ 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 :/ $\endgroup$
    – Kuba ♦
    May 12, 2015 at 8:14

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.