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
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
 

#!/usr/bin/env python3

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

    # Import special modules ...
    try:
        import matplotlib
        matplotlib.rcParams.update(
            {
                       "backend" : "Agg",                                       # NOTE: See https://matplotlib.org/stable/gallery/user_interfaces/canvasagg.html
                    "figure.dpi" : 300,
                "figure.figsize" : (9.6, 7.2),                                  # NOTE: See https://github.com/Guymer/misc/blob/main/README.md#matplotlib-figure-sizes
                     "font.size" : 8,
            }
        )
        import matplotlib.pyplot
    except:
        raise Exception("\"matplotlib\" is not installed; run \"pip install --user matplotlib\"") from None
    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; you need to have the Python module from https://github.com/Guymer/PyGuymer3 located somewhere in your $PYTHONPATH") from None

    # Define constants ...
    bins = [
        "all10/a11g",
        "all10/b10g",
        "all10/c10g",
        "all10/d10g",
        "all10/e10g",
        "all10/f10g",
        "all10/g10g",
        "all10/h10g",
        "all10/i10g",
        "all10/j10g",
        "all10/k10g",
        "all10/l10g",
        "all10/m10g",
        "all10/n10g",
        "all10/o10g",
        "all10/p10g",
    ]
    nx = 43200                                                                  # [px]
    ny = 21600                                                                  # [px]
    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]
    ]

    # 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."

    # **************************************************************************
    # *                 CREATE PNG IMAGES FROM REMOTE SOURCES                  *
    # **************************************************************************

    # Start session ...
    with pyguymer3.start_session() as sess:
        # Define ZIP file name and download it if it is missing ...
        zfile = "all10g.zip"
        if not os.path.exists(zfile):
            print(f"Downloading \"{zfile}\" ...")
            if not pyguymer3.download_file(sess, "https://www.ngdc.noaa.gov/mgg/topo/DATATILES/elev/all10g.zip", zfile):
                raise Exception("download failed", "https://www.ngdc.noaa.gov/mgg/topo/DATATILES/elev/all10g.zip") from None

    # Loop over elevation pairs ...
    for minElev, maxElev in pairs:
        # Check inputs ...
        if minElev != 0:
            raise Exception("currently this script assumes that the lower elevation is 0") from None

        # Deduce stub ...
        stub = f"{minElev:04d}m-{maxElev:04d}m"

        # Deduce PNG name ...
        pfile = f"{stub}.png"

        # Add to JSON dictionary ...
        data[stub] = {}
        data[stub]["__comment__"] = "A 30-arc-second (1-km) gridded, quality-controlled global Digital Elevation Model (DEM)"
        data[stub]["__projection__"] = "PlateCarree"
        data[stub]["__source__"] = "https://www.ngdc.noaa.gov/mgg/topo/globe.html"
        data[stub][stub] = pfile

        # Skip if PNG already exists ...
        if os.path.exists(pfile):
            continue

        print(f"Creating \"{pfile}\" ...")

        # Make map ...
        arr = numpy.zeros((ny, nx), dtype = numpy.int16)                        # [m]

        # Load dataset ...
        with zipfile.ZipFile(zfile, "r") as fObj:
            # Initialize index ...
            iy = 0                                                              # [px]

            # Loop over y-axis ...
            for i in range(4):
                # Initialize index ...
                ix = 0                                                          # [px]

                # Loop over x-axis ...
                for j in range(4):
                    # Define tile size ...
                    if i in [0, 3]:
                        nrows = 4800                                            # [px]
                    else:
                        nrows = 6000                                            # [px]
                    ncols = 10800                                               # [px]

                    # Load tile ...
                    tile = numpy.frombuffer(
                        fObj.read(bins[j + i * 4]),
                        dtype = numpy.int16
                    ).reshape(nrows, ncols)                                     # [m]

                    # Fill map ...
                    arr[iy:iy + tile.shape[0], ix:ix + tile.shape[1]] = tile[:, :]  # [m]

                    # Increment index ...
                    ix += ncols                                                 # [px]

                # Increment index ...
                iy += nrows                                                     # [px]

        # Load colormap and make levels ...
        cm = matplotlib.pyplot.get_cmap("jet")
        levels = []
        for level in range(maxElev + 1):
            r, g, b, a = cm(float(level) / float(maxElev))
            levels.append([numpy.uint8(255.0 * r), numpy.uint8(255.0 * g), numpy.uint8(255.0 * b)])

        # Make image ...
        img = numpy.zeros((ny, nx, 3), dtype = numpy.uint8)

        # Loop over y-axis ...
        for iy in range(ny):
            # Loop over x-axis ...
            for ix in range(nx):
                # Determine colour level ...
                level = min(maxElev, max(0, arr[iy, ix]))

                # Set pixel ...
                img[iy, ix, 0] = levels[level][0]
                img[iy, ix, 1] = levels[level][1]
                img[iy, ix, 2] = levels[level][2]

        # Save PNG ...
        pyguymer3.image.save_array_as_PNG(img, pfile, ftype_req = 0)
        pyguymer3.image.optimize_image(pfile, strip = True)

    # **************************************************************************
    # *            CREATE DOWNSCALED PNG IMAGES FROM LOCAL SOURCES             *
    # **************************************************************************

    # Loop over elevation pairs ...
    for minElev, maxElev in pairs:
        # Check inputs ...
        if minElev != 0:
            raise Exception("currently this script assumes that the lower elevation is 0") from None

        # Deduce stub ...
        stub = f"{minElev:04d}m-{maxElev:04d}m"

        # Deduce PNG name and skip if it is missing ...
        pfile1 = f"{stub}.png"
        if not os.path.exists(pfile1):
            continue

        # Loop over downscaled sizes ...
        for width in [512, 1024, 2048, 4096, 8192]:
            # Deduce downscaled PNG file name ...
            pfile2 = f"{stub}{width:04d}px.png"

            # Add to JSON dictionary ...
            data[stub][f"{stub}{width:04d}px"] = pfile2

            # Skip if downscaled PNG already exists
            if os.path.exists(pfile2):
                continue

            print(f"Creating \"{pfile2}\" ...")

            # Convert PNG to downscaled PNG ...
            pyguymer3.image.image2png(pfile1, pfile2, screenHeight = width, screenWidth = width, strip = True)

    # 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
 

{
    "0000m-1000m": {
        "0000m-1000m": "0000m-1000m.png", 
        "0000m-1000m0512px": "0000m-1000m0512px.png", 
        "0000m-1000m1024px": "0000m-1000m1024px.png", 
        "0000m-1000m2048px": "0000m-1000m2048px.png", 
        "0000m-1000m4096px": "0000m-1000m4096px.png", 
        "__comment__": "A 30-arc-second (1-km) gridded, quality-controlled global Digital Elevation Model (DEM)", 
        "__projection__": "PlateCarree", 
        "__source__": "https://www.ngdc.noaa.gov/mgg/topo/globe.html"
    }, 
    "0000m-2000m": {
        "0000m-2000m": "0000m-2000m.png", 
        "0000m-2000m0512px": "0000m-2000m0512px.png", 
        "0000m-2000m1024px": "0000m-2000m1024px.png", 
        "0000m-2000m2048px": "0000m-2000m2048px.png", 
        "0000m-2000m4096px": "0000m-2000m4096px.png", 
        "__comment__": "A 30-arc-second (1-km) gridded, quality-controlled global Digital Elevation Model (DEM)", 
        "__projection__": "PlateCarree", 
        "__source__": "https://www.ngdc.noaa.gov/mgg/topo/globe.html"
    }, 
    "0000m-3000m": {
        "0000m-3000m": "0000m-3000m.png", 
        "0000m-3000m0512px": "0000m-3000m0512px.png", 
        "0000m-3000m1024px": "0000m-3000m1024px.png", 
        "0000m-3000m2048px": "0000m-3000m2048px.png", 
        "0000m-3000m4096px": "0000m-3000m4096px.png", 
        "__comment__": "A 30-arc-second (1-km) gridded, quality-controlled global Digital Elevation Model (DEM)", 
        "__projection__": "PlateCarree", 
        "__source__": "https://www.ngdc.noaa.gov/mgg/topo/globe.html"
    }, 
    "0000m-4000m": {
        "0000m-4000m": "0000m-4000m.png", 
        "0000m-4000m0512px": "0000m-4000m0512px.png", 
        "0000m-4000m1024px": "0000m-4000m1024px.png", 
        "0000m-4000m2048px": "0000m-4000m2048px.png", 
        "0000m-4000m4096px": "0000m-4000m4096px.png", 
        "__comment__": "A 30-arc-second (1-km) gridded, quality-controlled global Digital Elevation Model (DEM)", 
        "__projection__": "PlateCarree", 
        "__source__": "https://www.ngdc.noaa.gov/mgg/topo/globe.html"
    }, 
    "0000m-5000m": {
        "0000m-5000m": "0000m-5000m.png", 
        "0000m-5000m0512px": "0000m-5000m0512px.png", 
        "0000m-5000m1024px": "0000m-5000m1024px.png", 
        "0000m-5000m2048px": "0000m-5000m2048px.png", 
        "0000m-5000m4096px": "0000m-5000m4096px.png", 
        "__comment__": "A 30-arc-second (1-km) gridded, quality-controlled global Digital Elevation Model (DEM)", 
        "__projection__": "PlateCarree", 
        "__source__": "https://www.ngdc.noaa.gov/mgg/topo/globe.html"
    }, 
    "0000m-6000m": {
        "0000m-6000m": "0000m-6000m.png", 
        "0000m-6000m0512px": "0000m-6000m0512px.png", 
        "0000m-6000m1024px": "0000m-6000m1024px.png", 
        "0000m-6000m2048px": "0000m-6000m2048px.png", 
        "0000m-6000m4096px": "0000m-6000m4096px.png", 
        "__comment__": "A 30-arc-second (1-km) gridded, quality-controlled global Digital Elevation Model (DEM)", 
        "__projection__": "PlateCarree", 
        "__source__": "https://www.ngdc.noaa.gov/mgg/topo/globe.html"
    }, 
    "0000m-7000m": {
        "0000m-7000m": "0000m-7000m.png", 
        "0000m-7000m0512px": "0000m-7000m0512px.png", 
        "0000m-7000m1024px": "0000m-7000m1024px.png", 
        "0000m-7000m2048px": "0000m-7000m2048px.png", 
        "0000m-7000m4096px": "0000m-7000m4096px.png", 
        "__comment__": "A 30-arc-second (1-km) gridded, quality-controlled global Digital Elevation Model (DEM)", 
        "__projection__": "PlateCarree", 
        "__source__": "https://www.ngdc.noaa.gov/mgg/topo/globe.html"
    }, 
    "0000m-8000m": {
        "0000m-8000m": "0000m-8000m.png", 
        "0000m-8000m0512px": "0000m-8000m0512px.png", 
        "0000m-8000m1024px": "0000m-8000m1024px.png", 
        "0000m-8000m2048px": "0000m-8000m2048px.png", 
        "0000m-8000m4096px": "0000m-8000m4096px.png", 
        "__comment__": "A 30-arc-second (1-km) gridded, quality-controlled global Digital Elevation Model (DEM)", 
        "__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. 512 px × 256 px (0.1 Mpx; 92.5 KiB)
  2. 1,024 px × 512 px (0.5 Mpx; 295.0 KiB)
  3. 2,048 px × 1,024 px (2.1 Mpx; 971.5 KiB)
  4. 4,096 px × 2,048 px (8.4 Mpx; 3.2 MiB)
  5. 8,192 px × 4,096 px (33.6 Mpx; 10.8 MiB)
  6. 43,200 px × 21,600 px (933.1 Mpx; 74.4 MiB)

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