Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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

enter image description here

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]

enter image description here

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

share|improve this question

1 Answer 1

up vote 11 down vote accepted

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
]

geographics screenshot

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?

share|improve this answer

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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