# 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. – Daniel Lichtblau Oct 17 '17 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
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. – Bob Hanlon Oct 18 '17 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. – Carl Woll Oct 18 '17 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? – creidhne Oct 20 '17 at 4:58

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

Using @Carl's example:

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