20
$\begingroup$

Some brief, only slightly important background. I am doing a research project using the data from NASA's GRACE mission. I wrote a short Perl script to take two data files and find the change in groundwater between the two dates. This gave me a set of 64,800 3D coordinates (One for every degree latitude and longitude on the Earth's surface). Using Mathematica, I created a ListDensityPlot to visualize the changes in groundwater. As you can see from the code below, the way I deal with clipping is pretty clumsy and doesn't look very good on the map. Otherwise, I am pretty happy with this plot. It pretty well shows everything I want it to. Most of the code courtesy of @Mr.Wizard.

den = ListDensityPlot[jpl200313,ColorFunction ->(ColorData["ThermometerColors"][1 - #] &),
    ClippingStyle -> {RGBColor[0.5, 0.02, 0.03],RGBColor[0, 0.01, 0.56]},
    PlotLegends ->BarLegend[Automatic,LegendMarkerSize -> 180,LegendFunction -> "Frame",
    LegendMargins -> 5,LegendLabel -> "Water Level Change (cm)"],PlotRange -> {-20, 20}];
prim = First@Cases[den, Graphics[a_, ___] :> a, {0, -1}, 1];
geo = GeoGraphics[{Opacity[0.6], prim},GeoBackground -> GeoStyling["StreetMapNoLabels"], 
   ImageSize -> 1000];
geo~Legended~den[[2]]

This shows global change in groundwater (in centimeters),from July 2003 to July 2013

The final piece that I would like to figure out is how to narrow down to specific countries while keeping the legend. Eventually I will build a table or possibly an animate function of several maps of the same country with time being the manipulatable variable. These pictures are from code courtesy of @FJRA.

southamerica =ListDensityPlot[jpl200313, AspectRatio -> 1/2, Frame->None, 
   PlotRangePadding -> 0, PlotRange -> {-20, 20},ColorFunction ->
    (ColorData["ThermometerColors"][1 - #] &)];
img1 = Rasterize[southamerica, "Image", RasterSize -> 360];
img2 = SetAlphaChannel[img1, .8];
geoplot = GeoGraphics[{GeoStyling[{"GeoImage", img2},GeoRange -> {{-90, 90}, {-180, 180}}],
    Polygon[EntityClass["Country", "SouthAmerica"]]},GeoBackground ->
    GeoStyling["StreetMapNoLabels"],GeoZoomLevel -> 3,GeoProjection -> "Equirectangular"]

The code for the picture of India is identical except for the name and the Entity function. South America changes in groundwater India changes in groundwater

Anyway, my big question at this point is whether or not the functionality of looking at individual countries can be combined with the readability of the first plot where I can add legends, titles labels etc. Thanks again!

$\endgroup$
  • $\begingroup$ What projection are you using for the plot above? $\endgroup$ – Mr.Wizard May 28 '15 at 0:32
  • $\begingroup$ That projection is completely square. That is the latitude and longitude lines make a perfect grid. $\endgroup$ – disc otter May 28 '15 at 0:34
  • $\begingroup$ I think based on that a simple method will work. Please attempt to apply my method and see if it aligns correctly. $\endgroup$ – Mr.Wizard May 28 '15 at 0:45
  • $\begingroup$ You're welcome. You may want to hold off on the Accept however as there are probably prettier ways to do it than simply setting Opacity[0.5]. $\endgroup$ – Mr.Wizard May 28 '15 at 0:54
  • 2
21
$\begingroup$

Please see the Utility function section for a concise summary.

An arbitrary density plot for the example:

den = DensityPlot[Sin[x] Sin[y], {x, -180, 180}, {y, -90, 90}]

enter image description here:

Extract the graphics primitives from the density plot:

prim = First @ Cases[den, Graphics[a_, ___] :> a, {0, -1}, 1];

Plot them directly with GeoGraphics while setting the desired GeoStyling for the GeoBackground:

GeoGraphics[
  {Opacity[0.8], prim},
  GeoBackground -> GeoStyling["ReliefMap"]
]

enter image description here

With GeoStyling["ContourMap"]:

enter image description here

ImageSize proves to be important; with "StreetMapNoLabels" an and and ImageSize of 512 or less no country borders are shown; 513 or greater and they are:

GeoGraphics[
  {Opacity[0.6], prim},
  GeoBackground -> GeoStyling["StreetMapNoLabels"], 
  ImageSize -> 600
]

enter image description here


Projections

To enable arbitrary projections we need to convert the plain coordinates in in the DensityPlot primitives to GeoPosition coordinates. prim as extracted above is a GraphicsComplex object which we can convert with:

prim2 = MapAt[GeoPosition @* Map[Reverse], prim, 1];

Now:

GeoGraphics[
  {Opacity[0.7], prim2},
  GeoBackground -> GeoStyling["StreetMapNoLabels"], 
  ImageSize -> 700,
  GeoProjection -> "Albers"
]

enter image description here

Legends

Including the legend from the original DensityPlot may be done like this:

den = DensityPlot[Sin[x] Sin[y], {x, -180, 180}, {y, -90, 90},
        PlotLegends -> Automatic];

prim = First @ Cases[den, Graphics[a_, ___] :> a, {0, -1}, 1];

geo = GeoGraphics[{Opacity[0.6], prim},
        GeoBackground -> GeoStyling["StreetMapNoLabels"], 
        ImageSize -> 600];

geo ~Legended~ den[[2]]

enter image description here


Utility function

The methods above may be combined into a single utility function.

toGeoGraphics[
  Shortest[opac : _?NumericQ : 0.6],
  opts : OptionsPattern[GeoGraphics]
][in_] :=
  With[{trans = If[MatchQ[OptionValue[GeoProjection], Automatic | "Equirectangular"], {},
     gc_GraphicsComplex :> MapAt[GeoPosition@*Map[Reverse], gc, 1]]},
    in /. Graphics[prim_, ___] :>
      GeoGraphics[{Opacity @ opac, prim /. trans}, opts, Options @ toGeoGraphics]
  ]

Define any default options that you want:

Options[toGeoGraphics] =
  {GeoBackground -> GeoStyling["StreetMapNoLabels"], 
   ImageSize -> 600};

Now use it like this:

DensityPlot[Sin[x] Sin[y], {x, -180, 180}, {y, -90, 90},
 PlotLegends -> Automatic] // 
   toGeoGraphics[GeoProjection -> "Mollweide"]

enter image description here

The first parameter of toGeoGraphics is the opacity; the remainder are any options you wish to pass to GeoGraphics, overriding defaults.

big = DensityPlot[Sin[x/38] Sin[y/25], {x, -180, 180}, {y, -90, 90}, 
  ColorFunction -> "CMYKColors", PlotPoints -> 100, MeshFunctions -> {#3 &, #3 &}, 
  Mesh -> {Range[-1, 1, 0.4], Range[-0.8, 0.8, 0.4]}, MeshStyle -> {Black, Dashed}, 
  PlotLegends -> Automatic];

big // toGeoGraphics[0.4, GeoProjection -> "Albers"]

enter image description here

$\endgroup$
  • $\begingroup$ @discotter Did you respect the First @ den part of my proposal? You need to extract the primitives from the Graphics expression or you will get that message. If you still get it I will take another look at this. $\endgroup$ – Mr.Wizard May 28 '15 at 1:32
  • $\begingroup$ No, I'm sorry. It did end up working. Apparently the PlotLegend piece of my code brought up the graphics error. Once I removed the PlotLegend command, I was able to combine them. Are there any projections I can use that don't have the continent and ocean names on them? $\endgroup$ – disc otter May 28 '15 at 1:56
  • $\begingroup$ @discotter I didn't consider legends; indeed a plot with a legend will have a format of Legended[Graphics[ . . . ], . . . ] and you would need to strip that too, e.g. First @ First @ den or den[[1, 1]]. Perhaps more robust would be FirstCase[den, Graphics[a_, ___] :> a, {}, {0, -1}]. Good question about the labels; let me check and get back to you. $\endgroup$ – Mr.Wizard May 28 '15 at 2:07
  • $\begingroup$ @discotter I realized that I was over-complicating things with Epilog. Please see the update. $\endgroup$ – Mr.Wizard May 28 '15 at 2:16
  • 2
    $\begingroup$ A word of caution: The legend does not account for the opacity of the density/contour plot, which might cause difficulties connecting the colors to their corresponding values. (+1 of course) $\endgroup$ – shrx May 29 '15 at 21:39
10
$\begingroup$

I would prefer to use the "GeoImage" styling, because you can use other projections when using it.

Let's say you have data for the whole world in a matrix:

data = Table[
  Sin[x Degree] Sin[y Degree], {y, -90, 90}, {x, -180, 180}]

Then you use ListDensityPlot:

den1 = ListDensityPlot[data, AspectRatio -> 1/2, Frame -> None, 
  PlotRangePadding -> 0];

And convert it to an image, and add some transparency level:

img1 = Rasterize[den1, "Image", RasterSize -> 360];
img2 = SetAlphaChannel[img1, .5]

enter image description here

Then you can use GeoStyling with "GeoImage":

GeoGraphics[{GeoStyling[{"GeoImage", img2}, 
   GeoRange -> {{-90, 90}, {-180, 180}}], 
  FilledCurve[
   GeoPath[{{-90, -180}, {90, -180}, {90, 0}, {90, 180}, {-90, 
      180}, {-90, 0}}, "Rhumb"]]}, GeoRange -> "World", 
 GeoBackground -> GeoStyling["StreetMapNoLabels"], GeoZoomLevel -> 2]

enter image description here

As I said, the good thing is that you can project it (but be careful with the GeoPath, needs to be modified for some projections that are not defined in the poles):

GeoGraphics[{GeoStyling[{"GeoImage", img2}, 
   GeoRange -> {{-90, 90}, {-180, 180}}], 
  FilledCurve[
   GeoPath[{{-86, -180}, {86, -180}, {86, 0}, {86, 180}, {-86, 
      180}, {-86, 0}}, "Rhumb"]]}, GeoRange -> "World", 
 GeoBackground -> GeoStyling["StreetMapNoLabels"], GeoZoomLevel -> 2, 
 GeoProjection -> "Mercator"]

enter image description here

Or clip only an specific area:

GeoGraphics[{GeoStyling[{"GeoImage", img2}, 
   GeoRange -> {{-90, 90}, {-180, 180}}], 
  Polygon[EntityClass["Country", "SouthAmerica"]]}, 
 GeoBackground -> GeoStyling["StreetMapNoLabels"], GeoZoomLevel -> 3, 
 GeoProjection -> "LambertAzimuthal"]

enter image description here

$\endgroup$
  • 1
    $\begingroup$ I added projections to my method as well. +1 for making me step up my game. :-) $\endgroup$ – Mr.Wizard May 29 '15 at 15:16
  • $\begingroup$ @FJRA This is great! However, I am still encountering one problem. Is there any way to keep a Plot Legend on the side of the map? For this data, it is pretty important that the viewer be able to see the units of what the density plot is showing. Is there a way to keep the plot legend separate and not combine it as part of the image? Thanks! $\endgroup$ – disc otter May 29 '15 at 18:08
  • $\begingroup$ @discotter I updated my answer for legends. $\endgroup$ – Mr.Wizard May 29 '15 at 18:29

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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