Background Images Of Elevation

metadata

This post is a follow up to my previous post about Replacing Cartopy’s Background Image. As discussed in the previous post, Cartopy’s default background image can be replaced with whatever you want: you just need to create a folder of PNG images with a JSON database file that describes them. To use this folder within your scripts you just need to set the environment variable CARTOPY_USER_BACKGROUNDS either by explicitly setting it in the shell before you call Python or by overriding it within Python by running something like os.environ["CARTOPY_USER_BACKGROUNDS"] = "/path/to/background/images/folder".

Recently I wanted to create some maps using elevation data and I came across the Global Land One-km Base Elevation (GLOBE) project, which is a 30-arc-second (1-km) gridded, quality-controlled global Digital Elevation Model (DEM). Their dataset is provided as a ZIP file containing a bunch of binary signed 16-bit integer tiles. To create images of the world I needed to load up all the tiles into RAM simultaneously (and then optionally clip them) to save them as a single PNG image.

The script below downloads the ZIP file (if it is missing) and then loops over pairs of elevations. For each elevation pair, such as “0m and 3000m”, it creates an array of the entire world and populates it with the elevation data from the tiles in the ZIP file. This array is then clipped to be between 0m and 3000m before being saved as a PNG image. Finally, the script creates smaller PNG images to allow quicker drawing before it creates the JSON database file that describes the whole folder.

  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
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
 

#!/usr/bin/env python3

# Use the proper idiom in the main module ...
# NOTE: See https://docs.python.org/3.12/library/multiprocessing.html#the-spawn-and-forkserver-start-methods
if __name__ == "__main__":
    # Import standard modules ...
    import argparse
    import json
    import os

    # Import special modules ...
    try:
        import numpy
    except:
        raise Exception("\"numpy\" is not installed; run \"pip install --user numpy\"") from None

    # Import my modules ...
    try:
        import pyguymer3
        import pyguymer3.image
    except:
        raise Exception("\"pyguymer3\" is not installed; run \"pip install --user PyGuymer3\"") from None

    # **************************************************************************

    # Create argument parser and parse the arguments ...
    parser = argparse.ArgumentParser(
           allow_abbrev = False,
            description = "Convert BIN files to PNG images.",
        formatter_class = argparse.ArgumentDefaultsHelpFormatter,
    )
    parser.add_argument(
        "--debug",
        action = "store_true",
          help = "print debug messages",
    )
    parser.add_argument(
        "--maximum-size",
        default = 250.0,
           dest = "maxSize",
           help = "the maximum size of image to make a PNG for (in mega-pixels)",
           type = float,
    )
    args = parser.parse_args()

    # **************************************************************************

    # NOTE: This script has been significantly reworked since the original was
    #       published. This script now uses the downscaled binary output created
    #       by my VGD project. You will have to run all the steps, up to and
    #       including, "downscale" before you can run this script. See:
    #         * https://github.com/Guymer/vgd

    # **************************************************************************

    # Define pairs of elevations ...
    pairs = [
        (0, 1000),                                                              # [m], [m]
        (0, 2000),                                                              # [m], [m]
        (0, 3000),                                                              # [m], [m]
        (0, 4000),                                                              # [m], [m]
        (0, 5000),                                                              # [m], [m]
        (0, 6000),                                                              # [m], [m]
        (0, 7000),                                                              # [m], [m]
        (0, 8000),                                                              # [m], [m]
        (0, 9000),                                                              # [m], [m]
    ]

    # Load colour tables and create short-hand ...
    with open(f"{pyguymer3.__path__[0]}/data/json/colourTables.json", "rt", encoding = "utf-8") as fObj:
        colourTables = json.load(fObj)
    turbo = numpy.array(colourTables["turbo"]).astype(numpy.uint8)

    # Create JSON dictionary ...
    data = {}
    data["__comment__"] = "JSON file specifying the image to use for a given type/name and resolution. Read in by cartopy.mpl.geoaxes.read_user_background_images."
    for minElev, maxElev in pairs:
        data[f"{minElev:04d}m-{maxElev:04d}m"] = {}
        data[f"{minElev:04d}m-{maxElev:04d}m"]["__comment__"] = f"A 30-arc-second (1-km) gridded, quality-controlled global Digital Elevation Model (DEM) scaled from {minElev:d} m to {maxElev:d} m"
        data[f"{minElev:04d}m-{maxElev:04d}m"]["__projection__"] = "PlateCarree"
        data[f"{minElev:04d}m-{maxElev:04d}m"]["__source__"] = "https://www.ngdc.noaa.gov/mgg/topo/globe.html"

    # **************************************************************************

    # Define the size of the dataset ...
    nx = 43200                                                                  # [px]
    ny = 21600                                                                  # [px]

    # Loop over scales ...
    for iscale in range(6):
        # Determine scale ...
        scale = pow(2, iscale)                                                  # [km]

        # Determine BIN file name ...
        bName = f"scale={scale:02d}km.bin"

        # Skip this scale if the BIN does not exist ...
        if not os.path.exists(bName):
            print(f"Skipping \"{bName}\" (BIN is missing).")
            continue

        # Find out how many mega-pixels there are at this scale ...
        mega = float((nx // scale) * (ny // scale)) / 1.0e6                     # [Mpx]

        # Skip this scale if the PNGs would be too big ...
        if mega > args.maxSize:
            print(f"Skipping \"{bName}\" (the PNGs would be {mega:.1f} Mpx).")
            continue

        # Loop over elevation pairs ...
        for minElev, maxElev in pairs:
            # Check inputs ...
            assert minElev == 0, "currently this script assumes that the lower elevation is 0 m"

            # Determine PNG file name ...
            pName = f"scale={scale:02d}km_minElev={minElev:04d}m_maxElev={maxElev:04d}m.png"

            # Add to JSON dictionary ...
            data[f"{minElev:04d}m-{maxElev:04d}m"][f"{scale:02d}km"] = pName

            # Skip this elevation pair if the PNG already exists ...
            if os.path.exists(pName):
                print(f"Skipping \"{pName}\" (PNG already exists).")
                continue

            print(f"Making \"{pName}\" ...")

            # Load dataset as 64-bit floats ...
            elev = numpy.fromfile(
                bName,
                dtype = numpy.int16,
            ).astype(numpy.float64).reshape(ny // scale, nx // scale, 1)        # [m]

            # Scale data from 0 to 255, mapping it from 0m to 2000m ...
            elev = 255.0 * elev / float(maxElev)
            numpy.place(elev, elev <   0.0,   0.0)
            numpy.place(elev, elev > 255.0, 255.0)
            elev = elev.astype(numpy.uint8)

            # Make PNG ...
            src = pyguymer3.image.makePng(
                elev,
                calcAdaptive = True,
                 calcAverage = True,
                    calcNone = True,
                   calcPaeth = True,
                     calcSub = True,
                      calcUp = True,
                     choices = "all",
                       debug = args.debug,
                         dpi = None,
                      levels = [9,],
                   memLevels = [9,],
                     modTime = None,
                    palUint8 = turbo,
                  strategies = None,
                      wbitss = [15,],
            )
            with open(pName, "wb") as fObj:
                fObj.write(src)

    # **************************************************************************

    # Save JSON dictionary ...
    with open("images.json", "wt", encoding = "utf-8") as fObj:
        json.dump(
            data,
            fObj,
            ensure_ascii = False,
                  indent = 4,
               sort_keys = True,
        )
You may also download “background-images-of-elevation-process.py” directly or view “background-images-of-elevation-process.py” on GitHub Gist (you may need to manually checkout the “main” branch).

Be careful using the script as it takes a long time and uses a lot of RAM to generate the initial PNG images before they are downsized. The GLOBE dataset in the ZIP file is 43,200 px × 21,600 px (which equates to a 889.9 mega-pixel image!). As a raw 2D 16-bit array that is 1.73 GiB of RAM; as a raw 3-channel image that is 2.61 GiB of RAM. On my Core i7 NAS it took around 1 hour 35 minutes to generate each initial PNG image - given the RAM usage (at least 4.34 GiB) I decided to not parallelise this script using multiprocessing. The JSON database file that is created is below.

 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
 

{
    "0000m-1000m": {
        "02km": "scale=02km_minElev=0000m_maxElev=1000m.png",
        "04km": "scale=04km_minElev=0000m_maxElev=1000m.png",
        "08km": "scale=08km_minElev=0000m_maxElev=1000m.png",
        "16km": "scale=16km_minElev=0000m_maxElev=1000m.png",
        "32km": "scale=32km_minElev=0000m_maxElev=1000m.png",
        "__comment__": "A 30-arc-second (1-km) gridded, quality-controlled global Digital Elevation Model (DEM) scaled from 0 m to 1000 m",
        "__projection__": "PlateCarree",
        "__source__": "https://www.ngdc.noaa.gov/mgg/topo/globe.html"
    },
    "0000m-2000m": {
        "02km": "scale=02km_minElev=0000m_maxElev=2000m.png",
        "04km": "scale=04km_minElev=0000m_maxElev=2000m.png",
        "08km": "scale=08km_minElev=0000m_maxElev=2000m.png",
        "16km": "scale=16km_minElev=0000m_maxElev=2000m.png",
        "32km": "scale=32km_minElev=0000m_maxElev=2000m.png",
        "__comment__": "A 30-arc-second (1-km) gridded, quality-controlled global Digital Elevation Model (DEM) scaled from 0 m to 2000 m",
        "__projection__": "PlateCarree",
        "__source__": "https://www.ngdc.noaa.gov/mgg/topo/globe.html"
    },
    "0000m-3000m": {
        "02km": "scale=02km_minElev=0000m_maxElev=3000m.png",
        "04km": "scale=04km_minElev=0000m_maxElev=3000m.png",
        "08km": "scale=08km_minElev=0000m_maxElev=3000m.png",
        "16km": "scale=16km_minElev=0000m_maxElev=3000m.png",
        "32km": "scale=32km_minElev=0000m_maxElev=3000m.png",
        "__comment__": "A 30-arc-second (1-km) gridded, quality-controlled global Digital Elevation Model (DEM) scaled from 0 m to 3000 m",
        "__projection__": "PlateCarree",
        "__source__": "https://www.ngdc.noaa.gov/mgg/topo/globe.html"
    },
    "0000m-4000m": {
        "02km": "scale=02km_minElev=0000m_maxElev=4000m.png",
        "04km": "scale=04km_minElev=0000m_maxElev=4000m.png",
        "08km": "scale=08km_minElev=0000m_maxElev=4000m.png",
        "16km": "scale=16km_minElev=0000m_maxElev=4000m.png",
        "32km": "scale=32km_minElev=0000m_maxElev=4000m.png",
        "__comment__": "A 30-arc-second (1-km) gridded, quality-controlled global Digital Elevation Model (DEM) scaled from 0 m to 4000 m",
        "__projection__": "PlateCarree",
        "__source__": "https://www.ngdc.noaa.gov/mgg/topo/globe.html"
    },
    "0000m-5000m": {
        "02km": "scale=02km_minElev=0000m_maxElev=5000m.png",
        "04km": "scale=04km_minElev=0000m_maxElev=5000m.png",
        "08km": "scale=08km_minElev=0000m_maxElev=5000m.png",
        "16km": "scale=16km_minElev=0000m_maxElev=5000m.png",
        "32km": "scale=32km_minElev=0000m_maxElev=5000m.png",
        "__comment__": "A 30-arc-second (1-km) gridded, quality-controlled global Digital Elevation Model (DEM) scaled from 0 m to 5000 m",
        "__projection__": "PlateCarree",
        "__source__": "https://www.ngdc.noaa.gov/mgg/topo/globe.html"
    },
    "0000m-6000m": {
        "02km": "scale=02km_minElev=0000m_maxElev=6000m.png",
        "04km": "scale=04km_minElev=0000m_maxElev=6000m.png",
        "08km": "scale=08km_minElev=0000m_maxElev=6000m.png",
        "16km": "scale=16km_minElev=0000m_maxElev=6000m.png",
        "32km": "scale=32km_minElev=0000m_maxElev=6000m.png",
        "__comment__": "A 30-arc-second (1-km) gridded, quality-controlled global Digital Elevation Model (DEM) scaled from 0 m to 6000 m",
        "__projection__": "PlateCarree",
        "__source__": "https://www.ngdc.noaa.gov/mgg/topo/globe.html"
    },
    "0000m-7000m": {
        "02km": "scale=02km_minElev=0000m_maxElev=7000m.png",
        "04km": "scale=04km_minElev=0000m_maxElev=7000m.png",
        "08km": "scale=08km_minElev=0000m_maxElev=7000m.png",
        "16km": "scale=16km_minElev=0000m_maxElev=7000m.png",
        "32km": "scale=32km_minElev=0000m_maxElev=7000m.png",
        "__comment__": "A 30-arc-second (1-km) gridded, quality-controlled global Digital Elevation Model (DEM) scaled from 0 m to 7000 m",
        "__projection__": "PlateCarree",
        "__source__": "https://www.ngdc.noaa.gov/mgg/topo/globe.html"
    },
    "0000m-8000m": {
        "02km": "scale=02km_minElev=0000m_maxElev=8000m.png",
        "04km": "scale=04km_minElev=0000m_maxElev=8000m.png",
        "08km": "scale=08km_minElev=0000m_maxElev=8000m.png",
        "16km": "scale=16km_minElev=0000m_maxElev=8000m.png",
        "32km": "scale=32km_minElev=0000m_maxElev=8000m.png",
        "__comment__": "A 30-arc-second (1-km) gridded, quality-controlled global Digital Elevation Model (DEM) scaled from 0 m to 8000 m",
        "__projection__": "PlateCarree",
        "__source__": "https://www.ngdc.noaa.gov/mgg/topo/globe.html"
    },
    "0000m-9000m": {
        "02km": "scale=02km_minElev=0000m_maxElev=9000m.png",
        "04km": "scale=04km_minElev=0000m_maxElev=9000m.png",
        "08km": "scale=08km_minElev=0000m_maxElev=9000m.png",
        "16km": "scale=16km_minElev=0000m_maxElev=9000m.png",
        "32km": "scale=32km_minElev=0000m_maxElev=9000m.png",
        "__comment__": "A 30-arc-second (1-km) gridded, quality-controlled global Digital Elevation Model (DEM) scaled from 0 m to 9000 m",
        "__projection__": "PlateCarree",
        "__source__": "https://www.ngdc.noaa.gov/mgg/topo/globe.html"
    },
    "__comment__": "JSON file specifying the image to use for a given type/name and resolution. Read in by cartopy.mpl.geoaxes.read_user_background_images."
}
You may also download “background-images-of-elevation-images.json” directly or view “background-images-of-elevation-images.json” on GitHub Gist (you may need to manually checkout the “main” branch).

The PNG file that it creates for the elevation pair “0m and 3000m”, is shown below.

Download:
  1. 1,350 px × 675 px (0.9 Mpx; 158.8 KiB)
  2. 2,700 px × 1,350 px (3.6 Mpx; 547.8 KiB)
  3. 5,400 px × 2,700 px (14.6 Mpx; 1.8 MiB)
  4. 10,800 px × 5,400 px (58.3 Mpx; 6.3 MiB)
  5. 21,600 px × 10,800 px (233.3 Mpx; 21.5 MiB)

The resulting folder of PNG images and JSON database file is entirely compatible with my function pyguymer3.geo.add_map_background() - enjoy!