# GeoListPlot and ListContourPlot

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

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?

• What code was used to produce these? Very hard to reverse engineer this from images. Oct 17, 2017 at 19:47

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]


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


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

Show[g, Epilog -> First @ plot]


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

• Using the option ContourShading -> None in definition of plot enables the map to show. Oct 18, 2017 at 0:36

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


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]


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]


The points on the map and contours now correspond.

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


• Try using GeoListPlot[GeoPosition[points], GeoProjection->"Albers"] to see why you can't simply use ListContourPlot with reversed coordinates in general. Oct 18, 2017 at 2:52
• @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? Oct 20, 2017 at 4:58

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@#, "Equirectangular"]&/@coords)[[All,1]];

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