Hypsometric area fun

Step 0: Get yourself a DEM

Imgur

If you don’t have one, there are plenty of online resources.

Step 1: Recode with GRASS

Imgur

We want polygons out of this thing, but if we let every pixel have a different value our polygon output will be a furious ball of nothing. I used GRASS r.recode, through QGIS, with a recode text file created with this bit of Python:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
minele = 200
maxele = 900
interval = 10

file = open("reclass.txt", "w")

x = minele

while x <= maxele:
    file.write(str(x) + ':' + str(x + interval) +
               ':' + str(round(x * 0.3048, 1)) + '\n')
    x = x + interval

file.close()

which output this:

1
2
3
4
5
6
7
8
200:210:61.0
210:220:64.0
220:230:67.1
230:240:70.1
240:250:73.2
250:260:76.2
260:270:79.2
...

Step 2: Majority Filter with SAGA

Imgur

r.recode helped clean up our DEM for polygon creation, but we can do better. SAGA’s Majority Filter takes areas too tiny to be useful and changes their value to the largest adjacent value. This will run for many hours. For our data (3.1’ pixels), I gave it a threshold of 0 and a kernel radius of 10, but adjust for your data as necessary. SAGA can be run via the SAGA GUI, the command line, or through QGIS. If running through the SAGA GUI or command line, note it creates a SAGA grid, which you’ll need to convert back to a TIFF.

Step 3: Hypsometric Areas with GRASS

Imgur

We’re ready to make polygons. I used GRASS r.to.vect to convert to vector polygons, making sure to check the Smooth corners of area features option for non-pixel-art polygon boundaries.

Step 4: Vectors to GeoJSON

In QGIS, right-click your polys and Save As to GeoJSON, making sure to convert to WGS84 if necessary.

Step 5: GeoJSON to Mapbox Vector Tiles

Tippecanoe to the rescue. Here’s the command I used, YMMV:

1
tippecanoe -o contours.mbtiles -P -y elev -E elev:mean -Z11 -z14 --drop-densest-as-needed contours.geojson

This created a 32MB MBTiles file. From a 2,300MB GeoJSON file. Gotta love vector tiles.

Step 6: Visual

Imgur

I used a GeoJSON polygon (in this case a county boundary) to be the water. The base tiles were served by mbtiles-server. Here be code:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
<!DOCTYPE html>
<html lang="en">

<head>
  <meta charset="UTF-8">
  <meta name="viewport" content="width=device-width, initial-scale=1.0">
  <meta http-equiv="X-UA-Compatible" content="ie=edge">
  <title>GL Test</title>

  <script src='https://api.tiles.mapbox.com/mapbox-gl-js/v0.45.0/mapbox-gl.js'></script>
  <link href='https://api.tiles.mapbox.com/mapbox-gl-js/v0.45.0/mapbox-gl.css' rel='stylesheet' />

  <style>
    * {
      box-sizing: border-box;
    }

    body {
      margin: 0;
      padding: 0;
    }

    #map {
      position: absolute;
      top: 0;
      bottom: 0;
      width: 100%;
    }

    #flood {
      position: absolute;
      bottom: 10px;
      width: 99%;
    }
  </style>
</head>

<body>
  <div id="map"></div>
  <input id="flood" type="range" min="80" max="250" step="0.3" value="80" oninput="setFloodHeight(this.value)" />


  <script>
    const mapStyle = {
      version: 8,
      sources: {
        "contours": {
          "type": "vector",
          "tiles": [
            "http://localhost:3000/contours/{z}/{x}/{y}.pbf"
          ],
          "maxzoom": 14,
          "minzoom": 12
        },
        "flood": {
          "type": "geojson",
          "data": "./county.geojson"
        }
      },
      layers: [{
          "id": "background",
          "type": "background",
          "paint": {
            "background-color": "rgb(239,239,239)"
          }
        },

        {
          "id": "hypsometric",
          "type": "fill-extrusion",
          "source": "contours",
          "source-layer": "contours",
          "minzoom": 11,
          "paint": {
            'fill-extrusion-color': {
              'property': 'elev',
              'stops': [
                [67, '#440154'],
                [87, '#472878'],
                [107, '#3e4a89'],
                [127, '#31688e'],
                [148, '#25838e'],
                [168, '#1e9e89'],
                [188, '#35b779'],
                [208, '#6cce59'],
                [228, '#b5de2c'],
                [248, '#fde725']
              ]
            },
            'fill-extrusion-height': ["get", "elev"],
            'fill-extrusion-opacity': 1,
            'fill-extrusion-base': 0
          }
        },

        {
          "id": "flood",
          "type": "fill-extrusion",
          "source": "flood",
          "minzoom": 11,
          "paint": {
            'fill-extrusion-color': "blue",
            'fill-extrusion-height': 80,
            'fill-extrusion-opacity': 0.6,
            'fill-extrusion-base': 0
          }
        }


      ]
    };

    const map = new mapboxgl.Map({
      container: 'map',
      style: mapStyle,
      attributionControl: false,
      center: [-80.84, 35.26],
      zoom: 12,
      minZoom: 12
    });

    function setFloodHeight(height) {
      map.setPaintProperty('flood', 'fill-extrusion-height', Number(height));
    }
  </script>
</body>

</html>

Resources