# Geographic background containing both outlines and contours: how?

I would like to create a map which will show lakes and rivers together with elevation contours and contour labels. What is the simplest way to do this with GeoGraphics in version 10?

The solution I have (below) looks much too convoluted.

I need to plot this point on the map:

pos = WeatherData["KP60", "Coordinates"]


This is good, except the background is not what I need:

g = GeoGraphics[{PointSize[Large], Darker@Red, Point[pos]}]


Since I can't seem to mix two GeoStyles, I tried overlaying a polygon that contains the contours. This lead to the following:

{range, projection} = {GeoRange, GeoProjection} /. Options[g, {GeoRange, GeoProjection}]

{a, b, d, c} = Tuples[Reverse@range]

GeoGraphics[{{GeoStyling["ContourMap", ContourShading -> None,  ContourLabels -> True], Polygon[{a, b, c, d}]}, {PointSize@Large, Darker@Red, Point[pos]}}, GeoRange -> range]


However, now the projection has changed, and adding GeoProjection -> projection causes the polygon not to show. Is there a fix for this?

-

The problem seems to be related to mixing geopositions and raw long/lat pairs. pos is expressed as a GeoPosition and the polygon is expressed as raw pairs. The graphics are better behaved if all coordinates are expessed in the same fashion.

Option one: use geopositions throughout

Change the assignment to {a, b, d, c} from this:

{a, b, d, c} = Tuples[Reverse@range]
(* {{-110.652, 44.3787}, {-110.652, 44.7093}, {-110.19, 44.3787}, {-110.19, 44.7093}} *)


to this:

{a, d, b, c} = GeoPosition /@ Tuples[range]
(* {GeoPosition[{44.3787, -110.652}], GeoPosition[{44.3787, -110.19}],
GeoPosition[{44.7093, -110.652}], GeoPosition[{44.7093, -110.19}]}*)


Option two: use raw long/lat pairs throughout

Change the assignment to pos from this:

pos = WeatherData["KP60", "Coordinates"]
(* GeoPosition[{44.544, -110.421}] *)


to this:

pos = Reverse @ First @ WeatherData["KP60", "Coordinates"]
(* {-110.421, 44.544} *)


Why do mixed coordinates behave like this?

I don't know, that remains to be seen... :)

Interestingly, if the Point uses a raw pair and Polygon uses geopositions (option 3: the reverse situation from the original code), then the result is well-behaved. Perhaps there is a bug in Polygon?

If we are using the mixed coordinates from the original code, the following expression suggests that there is some kind of projection mismatch:

GeoGraphics[
{ { GeoStyling["ContourMap", ContourShading -> None,  ContourLabels -> True]
, Polygon[{a, b, c, d}]
}
, {PointSize@Large, Orange, Point[pos]}
}
, GeoProjection -> projection
]


The GeoRange option has been deleted so that we can see the entire range. The GeoProjection option has been added. Note how the contoured polygon is being plotted in a different location despite having overlapping coordinates. Is the wrong projection being used?

-