maplibre / maplibre-style-spec

MapLibre Style Specification & Utilities
https://maplibre.org/maplibre-style-spec/
Other
92 stars 65 forks source link

Design Proposal: Contour Line Source from Raster DEM Tiles #583

Open msbarry opened 7 months ago

msbarry commented 7 months ago

Design Proposal: Contour Line Source from Raster DEM Tiles

Motivation

Give users a built-in way to render contour lines in maplibre from the same DEM tiles that are already used for terrain and hillshading, like this:

image

Proposed Change

Create a new contour source type in maplibre style spec that takes a raster-dem source as input and generates contour isolines as output that can be styled using line layers, for example:

sources: {
  dem: {
    type: "raster-dem",
    encoding: "terrarium",
    tiles: ["https://elevation-tiles-prod.s3.amazonaws.com/terrarium/{z}/{x}/{y}.png"],
    maxzoom: 13,
    tileSize: 256,
  },
  contours: {
    type: "contour",
    source: "dem",
    unit: "feet" | "meters" | number, // default=meters, for custom unit use length of the unit like unit: 1.8288 for fathoms
    // similar syntax to ["step", ["zoom"], ...] style expression 
    // to define contour interval by zoom level
    intervals: [
      200, // 200m interval <=z11
      12, 100, // 100m interval at z12 and z13
      14, 50, // 50m interval at z14
      15, 20 // 20m interval >= z15
    },
    // put a "major=false/true" tag on every Nth line by zoom so styles
    // can highlight major/minor lines differently
    majorMultiplier: [
      5, // every 5th line at < z14
      14, 4, // every 4th line at z14
      15, 5 // every 5th line for >= z15
    ],
    // minzoom inferred from raster-dem source and maxzoom determined automatically by maplibre
    // overzoom z10 tiles to generate z11 contour lines, z11 to make z12, etc...
    overzoom: 1,
  },
},

The generated isolines will have these attributes:

Layers can refer to the contours with source: contours but they can omit sourceLayer.

This offloads details about how to retrieve and parse DEM tiles to the DEM source definition, and gives style layers the flexibility to render any number of visible lines derived from that contour source.

I've already prototyped this in the maplibre-contour plugin which I'm using for contour lines on onthegomap.com. Here are some of the issues I had to work through to get these contours to look nice:

DEM "overzooming" (smoothing)

The contour lines look blocky when you zoom in much further past the maximum zoom for a DEM source, but they can look nice and smooth if you "overzoom" the DEM tiles by applying iterative bilinear interpolation before generating isolines. For example for onthegomap I use 512px z11 tiles, but overzoom the z11 tiles up to z15 so that the contour lines look smooth even at high zooms. This is why the proposal lets you specify a maxzoom that is higher than the maxzoom of the raster-dem source.

Also to generate smooth contour lines at the border between tiles, the algorithm needs to look at adjacent tiles. This means you need 9 DEM tiles to render a single contour tile. To mitigate this, the overzoom parameter lets you use overzoomed DEM tiles from a lower zoom level to generate contours at the current zoom level, for example overzoom=1 means use the top-left, top-right, botom-left, or bottom-right z10 tile to render a z11 contour line tile. This means you only need 4 DEM tiles to render a single contour tile:

image

Contour levels and units

The user needs to be able to choose what elevations to draw contour lines at, which changes by zoom level (rendering every contour would get too expensive at low zooms in hilly areas). They may also designate "major" and "minor" levels, for example generate thin contour lines every 200m but bold every 1000m. For now we will push this to layers that use the style, but in the future we can either add a major/minor designation to ticks, or pass-through the level and interval so styles can highlight every 5th or 10th line or something.

The unit attribute multiplies raw meter values by a certain amount to change the unit, for example unit=feet changes from meters to feet. When you click the distance indicator on onthegomap, it toggles between unit=meters and unit=feet. You could also set unit to a custom value for less common units like unit=1.8288 for fathoms.

Performance and Bundle Size

I've already implemented the smoothing logic and isoline generation in the maplibre-contour plugin so we would just need to bring that into maplibre-gl-js and port into the native projects. The overall plugin is 33kb (11kb gzipped) but most of that is replicating the web worker communication, cancelable message passing, and vector tile encoding that maplibre-gl-js already has. The actual smoothing+isoline business logic is only 3.7kb (1.6kb gzipped).

The isoline generation algorithm was derived from d3-contour but is much more efficient because it generates isolines in a single pass through the DEM tile instead of using a pass per contour level. For onthegomap users on a range of devices (mostly mobile phones) overzooming a 512px dem tile and generating isolines takes:

API Modifications

This should only change the style spec, but shouldn't require any js or native API changes, unless we wanted to expose the default contour layer, elevation or level key as constants?

Migration Plan and Compatibility

This is new functionality, so no migration is necessary.

Rejected Alternatives

Build a plugin for this

I maintain the maplibre-contour plugin which already lets you do this by using the addProtocol integration, but it has a few downsides:

  1. It's an extra step to install: rendering contour lines is a common use-case that users should be able to do by default
  2. 90% of the plugin is duplicating things that maplibre already does like spawning a web worker, communicating with it using cancelable messages, and decoding DEM tiles. The actual code for computing the contours is a small fraction of the overall plugin.
  3. It has to do a wasteful extra step of encoding the result as vector tile bytes only for maplibre to decode immediately after in its own web workers (see https://github.com/onthegomap/maplibre-contour/blob/main/architecture.png)
  4. It can't make use of other registered maplibre request interceptors or protocols like DEM tiles served out of a pmtiles archive
  5. It doesn't work in maplibre-native

Pre-render contour lines

You can render contour line vector tiles ahead of time and serve those for the planet, this will save some browser CPU cycles but rendering them on the fly from DEM tiles has a few advantages:

  1. There are a lot of parameters you can tweak when generating contour lines from elevation data like units, thresholds, and smoothing parameters. Pre-generated contour vector tiles require 100+gb of storage for each variation you want to generate and host. Generating them on-the-fly in the browser gives infinite control over the variations you can use on a map from the same source of raw elevation data that maplibre uses to render terrain and hillshade.
  2. You're likely already downloading DEM tiles for hillshading and terrain, so this eliminates the extra bandwidth used to download those vector tiles.

Implement as a new layer type

We could implement this as a new layer type instead of a source type, but that would tightly couple the display parameters to the logic for how contour lines are generated, and potentially require us to generate the contours in multiple passes over the source DEM data. It seems cleaner to generate contour lines so you can generate as many layers as you want from them afterwards.

Take a DEM tile source URL as input

A contour layer could take as input tiles: ["server.com/{z}/{x}/{y}.png"], but there are a lot of different knobs to tune for how these are interpreted, so by depending on a DEM source we re-use the DEM source control all of those parameters.

voncannon commented 7 months ago

Naming is hard. So adding more terminology to the docket...

I think traditionally the thicker contour lines are called index contours:

  1. page 9 https://pubs.usgs.gov/circ/1955/0368/report.pdf
  2. https://www3.nd.edu/~cneal/planetearth/Lab-SurfaceHydrology/TopoMaps.html

Following the example below, one advantage here is both intervals and indexIntervals track the same unit as specified .. rather than one being a unit and one a interval/multiplier/skip/etc, respectively.

N/multiplier/skip can be computed per zoom level under the covers: multiplier = index_interval / interval

sources: {
  dem: {
    type: "raster-dem",
    encoding: "terrarium",
    tiles: ["https://elevation-tiles-prod.s3.amazonaws.com/terrarium/{z}/{x}/{y}.png"],
    maxzoom: 13,
    tileSize: 256,
  },
  contours: {
    type: "contour",
    source: "dem",
    unit: "feet" | "meters" | number, // default=meters
    // similar syntax to ["step", ["zoom"], ...] style expression 
    // to define contour interval by zoom level
    intervals: [
      200, // 200m interval <=z11
      12, 100, // 100m interval at z12 and z13
      14, 50, // 50m interval at z14
      15, 20 // 20m interval >= z15
    },
    indexIntervals: [
      5, // every 1000m at < z14
      14, 200, every 200m at z14
      15, 100 // every 100m for >= z15
    ],
    minzoom: 11,
    // overzoom z13 DEM tiles up to z16 to generate smooth high-zoom contour lines:
    maxzoom: 16,
    // overzoom z10 tiles to generate z11 contour lines, z11 to make z12, etc...
    overzoom: 1,
  },
},
HarelM commented 7 months ago

If the indexInterval (the word index is far from something intuitive to me BTW, but it's not important for this comment) isn't a full step of the interval things become tricky. What happen if the interval is 100m and you want to place an index line in 250m? Is it valid? If not, do you get a console warning? The idea behind the multiplier is that effectively, the above can't happen. It forces integer multiplier. I also believe the common case would be the 5th line, or something similar, for all zoom levels which doesn't require to rewrite the intervals again, but simply write “[5]“. I have thought about what you are proposing before posting my suggestion, but given the above (which can be solve obviously, not saying they are a blocker) I believe what I suggested would be better.

voncannon commented 7 months ago

Great thoughts, I agree, that could make a mess of things where as the multiplier directly in the options wouldn't allow the user to get in a pickle as easy.

With that said, the index contour terminology still goes. But either way for me, major/index. I just think traditionally index is what has been used.

msbarry commented 7 months ago

I did find that there is one place where a source definition already uses a style expression for the geojson clusterProperties field: https://github.com/maplibre/maplibre-style-spec/blob/69300a8f5d23b2c5f4040793b9c737e746e5bcde/src/validate/validate_source.ts#L64-L78

Think we should continue with the abbreviated step-zoom syntax? Or turn intervals and major/index into expressions? We still might run into issues because it's not operating in the context of an individual feature...

HarelM commented 7 months ago

I believe it makes more sense to have an expression there if this is not the first time an expression is used, the documentation will need to provide a good example and describe what properties are available in this context (probably only zoom, if I needed to guess).

If we encounter serious issues with the implementation we can reduce the complexity of the spec, but I think it makes sense to try.

zstadler commented 7 months ago

I does seem that the term "Index" is used by English speakers, at least in the US. My efforts to find the UK term were unsuccessful... Perhaps a index-multiplier keyword would make things clear for most people.

It uses kebab-case, similar to other keywords in the style spec, such as icon-text-fit-padding.

msbarry commented 6 months ago

What do people think of omitting maxzoom from the contour source definition? Getting into the gl-js implementation, there are 4 stages of the data transformation:

  1. raw DEM tiles
  2. "overzoomed" DEM tiles with smoothed values ready to generate contours from
  3. vector tile with contours generated from overzoomed DEM tile
  4. overzoomed vector tile

So the overzooming to get from a z12 raster DEM tile to a z16, 17, 18+ vector tile could happen in either step 2 or 4. If you set the source maxzoom parameter correctly, the results will look nearly identical, but if you set them wrong they could look bad.

I'd propose that we let maplibre manage how much happens in 2 vs. 4 in order to optimize performance and visual appearnce without a user needing to know how the internals work. Is it weird for a source not to have a maxzoom though?

msbarry commented 6 months ago

And for unit I was initially thinking of it as a multiplier (like 3.28084 for feet) but wondering if it makes more sense for it to be the length of the actual unit (0.3048 for feet, 1.8288 for fathoms)?

msbarry commented 6 months ago

Making progress on the gl-js implementation... I've got an initial rough proof of concept working:

image

Still need to figure out source tile caching and handling border pixels but I have found that using expressions in the source definition for intervals and major/index works well (both on the evaluation side, and validating that it only uses zoom as an input)

HarelM commented 6 months ago

I think it makes sense to relay on the definitions of the referenced source to simplify stuff. This is the first source that references another source, so it's the first of it's kind.

Light already uses expression in the position property, so I believe using expression here for the interval would make sense as well.

zstadler commented 6 months ago

What do people think of omitting maxzoom from the contour source definition?

Perhaps there is no need for minzoom, maxzoom and overzoom in a contour source definition. Let the contour source implementation do its best given the minzoom and maxzoom of the Raster DEM source it is based on, and let each style layers define its minzoom and maxzoom for what needs to be shown in the map.

No harm will be done if the style layer has no minzoom and the contour source cannot perform under-zoom. Similar to hillshade style layers, it will just not be shown.

msbarry commented 6 months ago

Perhaps there is no need for minzoom, maxzoom and overzoom in a contour source definition. Let the contour source implementation do its best given the minzoom and maxzoom of the Raster DEM source it is based on, and let each style layers define its minzoom and maxzoom for what needs to be shown in the map.

That sounds good to omit minzoom/maxzoom - "overzoom" is a bit different though - it uses lower zoom tile to generate higher zoom contours so that 1) the results look smoother and 2) there's less data to fetch to cover the viewport, for example:

overzoom=0 image overzoom=1 image overzoom=2 image overzoom=3 image

In practice I've found an "overzoom" of 1 works well with maplibre - an overzoom of 0 has too many kinks in the lines for mapilbre's label placement to find suitable spots to put the height labels (although there are probably improvements we could make to label placement to automatically "de-kink" spots like this to improve line label placement)

Maybe this is a concept we'd want to enable across other sources? Or maybe we it should be called something like "smoothing"

zstadler commented 6 months ago

It seems like some peaks loose elevation when smoothing/overzoom is applied:

image

The size and shape of the 7,000' contour line around this peak also changes significantly as smoothing is increased.

In practice I've found an "overzoom" of 1 works well with maplibre - an overzoom of 0 has too many kinks in the lines for mapilbre's label placement to find suitable spots to put the height labels (although there are probably improvements we could make to label placement to automatically "de-kink" spots like this to improve line label placement)

A similar issue exists for rivers and potentially other natural features. I mitigated that by increasing text-max-angle and text-letter-spacing.

If the elevation labels are the reason for implementing smoothing/overzooming, then I suggest addressing the issue in the symbol layer of the labels, rather than decreasing the accuracy of the contour lines.

P.S., I would love to see "smoothing" available for line labels in general. Perhaps this is a topic for a different discussion.

zstadler commented 6 months ago

Some maps, including OpenTopoMap and later also Israel Hiking Map, help the reader to identify the "up" direction by orienting the labels according the direction of the slope:

I wonder how this label would behave when "text-keep-upright": false is set in the contour labels layer definition.

image

If the label keeps this orientation, please reverse the direction of the contour lines. The reversal would only affect label layers that use "text-keep-upright": false and enable achieving the desired orientation.