How can I use multiple overlaid map layers (e.g. basemap, labels) in GeoGraphics?

Question

We can use GeoServer to use map tiles from an external service. There are many compatible services:

• Which external services are compatible with GeoServer?

These services often provide unlabelled background maps, as well as various layers to superimpose on a base map (e.g. labels only, borders only, etc.).

How can we use multiple layers, possibly loaded from multiple different servers, in a way that integrates well with the GeoGraphics functionality and also work with DynamicGeoGraphics?

---

For example, using the definitions from the above linked answer,

base = GeoGraphics[{Entity["Country", "Hungary"]}, 
  GeoServer -> stamen["watercolor"]]

labels = GeoGraphics[{Entity["Country", "Hungary"]}, 
  GeoServer -> carto["light_only_labels"]]

We could layer these together manually using one of multiple possible ways, e.g.

Overlay[{base, labels}]

But this is not a good solution because it requires multiple GeoGraphics calls that are independent, and does not work with functions like DynamicGeoGraphics.

Answer 1

It's possible to superimpose multiple layers within the same GeoGraphics if the layers are "geoprimitives" for which you can define different geostyles (including opacity and the choice of geoserver). Here is a simple example:

GeoGraphics[{{GeoStyling["StreetMap", Opacity[1], 
    GeoServer -> "http://tile.stamen.com/toner/`1`/`2`/`3`.png"], 
   Polygon[Entity["Country", "Hungary"]]}, {GeoStyling["StreetMap", 
    Opacity[0.2], 
    GeoServer -> "http://a.tile.openstreetmap.org/`1`/`2`/`3`.png"
    ], Polygon[Entity["Country", "Hungary"]]}}, 
 GeoBackground -> GeoStyling[Opacity[0.]], ImageSize -> 500]

(Instead of Polygon you can use other geographic primitives like in this post.)

The difficulty of course is to set the right opacities for the layers (when you have to deal with just more than labels like in your trivial example).

You can try the previous example with DynamicGeoGraphics, it does not seem to work, zooming is not possible, the layers are not dynamically updated. And you'll notice that you cannot get rid of the GeoBackground anymore which actually is the only one to update ! Clearly, DynamicGeoGraphics need some more dev ...

---

Here's a sample application to show MODIS snow cover data in China:

GeoGraphics[{
  EdgeForm[Black], 
  GeoStyling["StreetMap", 
      GeoServer -> {"https://map1.vis.earthdata.nasa.gov/wmts-webmerc/MODIS_Terra_Snow_Cover/default//GoogleMapsCompatible_Level8/`1`/`3`/`2`.png", 
                   "TileDataType" -> "PNG", "ZoomRange" -> {1, 8}}
  ], 
  Polygon[Entity["Country", "China"]]
 }, 
 GeoBackground -> GeoStyling["ReliefMap"], ImageSize -> Large
]