3
$\begingroup$

I'm doing a research and I have to add to a GeoListPlot a ListContourPlot, but it doesn't work with simple Show and I have no idea how to make it working

GeoListPlot image ListContourPlot plot

As you can see these points doesn't match well... but they're plotted from the same data. Does anyone know how to make it working?

$\endgroup$
  • 3
    $\begingroup$ What code was used to produce these? Very hard to reverse engineer this from images. $\endgroup$ – Daniel Lichtblau Oct 17 '17 at 19:47
4
$\begingroup$

GeoListPlot uses GeoProjection/GeoModel when constructing the output. This means the data points do not match the latitude-longitude coordinates as shown in a equirectangular projection, which is effectively what ListContourPlot uses. One possibility to fix this is to extract the projection coordinates from GeoListPlot. Here is some sample data:

SeedRandom[1]
coords = Thread[{RandomReal[{35,45}, 10], RandomReal[{-90,-80}, 10]}];
coordValues = RandomInteger[10, 10];

Here is a method to get a list of rules to convert the latitude-longitude coordinates coords into the projection coordinates used by GeoListPlot:

g = GeoListPlot[Tooltip[GeoPosition[#], #]& /@ coords]
rules = Cases[g, Tooltip[Inset[_, geo_], latlong_List] :> latlong->geo, Infinity]

enter image description here

{{35.6574, -81.7484} -> {-81.7484, 38.2106}, {36.1142, -82.5134} -> {-82.5134, 38.7744}, {36.878, -87.5251} -> {-87.5251, 39.7246}, {37.3115, -84.2194} -> {-84.2194, 40.2681}, {37.4136, -80.2283} -> {-80.2283, 40.3965}, {38.9601, -87.0713} -> {-87.0713, 42.3641}, {40.4225, -80.7472} -> {-80.7472, 44.2647}, {42.0047, -87.9195} -> {-87.9195, 46.3682}, {42.8953, -85.7715} -> {-85.7715, 47.5752}, {43.1739, -87.8817} -> {-87.8817, 47.9564}}

Now, we construct a dataset using the geo coordinates and the coordValues above:

dataset = Join[coords /. rules, List/@coordValues, 2];

Here is a ListContourPlot of dataset:

plot = ListContourPlot[dataset, Mesh->All]

enter image description here

Finally, we can add the ListContourPlot output to the GeoGraphics object produced by GeoListPlot:

Show[g, Epilog -> First @ plot]

enter image description here

Notice how the ListContourPlot output matches up with the coordinates in the GeoListPlot output.

| improve this answer | |
$\endgroup$
  • $\begingroup$ Using the option ContourShading -> None in definition of plot enables the map to show. $\endgroup$ – Bob Hanlon Oct 18 '17 at 0:36
3
$\begingroup$

The problem lies in how GeoListPlot and ContourListPlot interpret their coordinates. GeoListPlot displays latitude on the vertical axis and longitude on the horizontal axis. ContourListPlot displays X values on the horizontal axis and Y values on the vertical axis. The order of the horizontal and vertical coordinates are reversed.

Here's a list of points, where each point is latitude, longitude, and elevation in meters.

points = {
  {37.80241, -122.4432, 7.},
  {37.77555, -122.4279, 46.},
  {37.80316, -122.4271, 26.},
  {37.80600, -122.4375, 5.},
  {37.78155, -122.4362, 44.},
  {37.78379, -122.4573, 63.},
  {37.77206, -122.4403, 76.},
  {37.76710, -122.4212, 9.},
  {37.76705, -122.4160, 4.},
  {37.79626, -122.4434, 47.}
};

We can plot the points on a map with:

map = GeoListPlot[GeoPosition[points]]

GeoListPlot

Here's a contour plot of these points. This demonstrates the problem. Notice that negative longitude values are plotted on the vertical axis, and the latitude values are on the horizontal axis.

ListContourPlot[points, Mesh -> All]

ListContourPlot with reversed axes

Solve the problem by simply reversing the X and Y points for ListContourPlot.

xypoints = #[[{2, 1, 3}]] & /@ points;
contours = ListContourPlot[xypoints, Mesh -> All, PlotLegends -> Automatic]

ListContourPlot with corrected axes

The points on the map and contours now correspond.

Show[map, ListContourPlot[xypoints, Mesh -> All, ContourShading -> None]]

enter image description here

| improve this answer | |
$\endgroup$
  • 1
    $\begingroup$ Try using GeoListPlot[GeoPosition[points], GeoProjection->"Albers"] to see why you can't simply use ListContourPlot with reversed coordinates in general. $\endgroup$ – Carl Woll Oct 18 '17 at 2:52
  • $\begingroup$ @Carl Yes, the swapping-coordinates method only works for the Mercator projection, but it has an advantage in that the ListContourPlot axes show the correct latitude and longitude without the need to find corrected values for DataRange. For example, the southernmost tip of Lake Michigan is at 42 N, and not nearly 46 N as the transformation plots imply. Did the question ask for a generalized coordinate transform, or rather, simply a way to match a GeoListPlot to a ListContourPlot? $\endgroup$ – creidhne Oct 20 '17 at 4:58
3
$\begingroup$

You can use a combination of GeoPosition and GeoGridPosition to get the projection of coords.

Using @Carl's example:

SeedRandom[1]
coords = Thread[{RandomReal[{35,45}, 10], RandomReal[{-90,-80}, 10]}];
coordValues = RandomInteger[10, 10];

Create projected coordinates for coords and combine with coordValues:

procejtedcoords = (GeoGridPosition[GeoPosition@#, "Mercator"]&/@coords)[[All,1]];

lcpdata = Join[procejtedcoords, List /@ coordValues, 2];
Show[g, ListContourPlot[lcpdata, Mesh -> All, BaseStyle->Opacity[.5]]]

enter image description here

| improve this answer | |
$\endgroup$

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.