Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
a script to create background images of elevation
#!/usr/bin/env python2
# -*- coding: utf-8 -*-
# Import modules ...
# NOTE: http://matplotlib.org/faq/howto_faq.html#matplotlib-in-a-web-application-server
import json
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot
import numpy
import os
import requests
import subprocess
import zipfile
# Import my modules ...
try:
import pyguymer
except:
raise Exception("you need to have the Python module from https://github.com/Guymer/PyGuymer located somewhere in your $PYTHONPATH")
# 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]
]
# 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 ...
sess = requests.Session()
sess.allow_redirects = True
sess.headers.update({"Accept" : "text/html,application/xhtml+xml,application/xml;q=0.9,*/*;q=0.8"})
sess.headers.update({"Accept-Language" : "en-GB,en;q=0.5"})
sess.headers.update({"Accept-Encoding" : "gzip, deflate"})
sess.headers.update({"DNT" : "1"})
sess.headers.update({"Upgrade-Insecure-Requests" : "1"})
sess.headers.update({"User-Agent" : "Mozilla/5.0 (Macintosh; Intel Mac OS X 10.14; rv:68.0) Gecko/20100101 Firefox/68.0"})
sess.max_redirects = 5
# Define ZIP file name and download it if it is missing ...
zfile = "all10g.zip"
if not os.path.exists(zfile):
print "Downloading \"{:s}\" ...".format(zfile)
if not pyguymer.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")
# Loop over elevation pairs ...
for pair in pairs:
# Check inputs ...
if pair[0] != 0:
raise Exception("currently this script assumes that the lower elevation is 0")
# Deduce stub ...
stub = "{:04d}m-{:04d}m".format(pair[0], pair[1])
# Deduce PNG name ...
pfile = "{:s}.png".format(stub)
# 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 "Creating \"{:s}\" ...".format(pfile)
# Make map ...
map = 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 xrange(4):
# Initialize index ...
ix = 0 # [px]
# Loop over x-axis ...
for j in xrange(4):
# Define tile size ...
if i == 0 or i == 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 ...
map[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 xrange(pair[1] + 1):
r, g, b, a = cm(float(level) / float(pair[1]))
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 xrange(ny):
# Loop over x-axis ...
for ix in xrange(nx):
# Determine colour level ...
level = min(pair[1], max(0, map[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 ...
pyguymer.save_array_as_PNG(img, pfile, ftype_req = 0)
subprocess.check_call(
["exiftool", "-overwrite_original", "-all=", pfile],
stderr = open(os.devnull, "wt"),
stdout = open(os.devnull, "wt")
)
subprocess.check_call(
["optipng", pfile],
stderr = open(os.devnull, "wt"),
stdout = open(os.devnull, "wt")
)
# End session ...
sess.close()
# ******************************************************************************
# * CREATE DOWNSCALED PNG IMAGES FROM LOCAL SOURCES *
# ******************************************************************************
# Loop over elevation pairs ...
for pair in pairs:
# Check inputs ...
if pair[0] != 0:
raise Exception("currently this script assumes that the lower elevation is 0")
# Deduce stub ...
stub = "{:04d}m-{:04d}m".format(pair[0], pair[1])
# Deduce PNG name and skip if it is missing ...
pfile1 = "{:s}.png".format(stub)
if not os.path.exists(pfile1):
continue
# Loop over downscaled sizes ...
for width in [512, 1024, 2048, 4096]:
# Deduce downscaled PNG file name ...
pfile2 = "{:s}{:04d}px.png".format(stub, width)
# Add to JSON dictionary ...
data[stub]["{:s}{:04d}px".format(stub, width)] = pfile2
# Skip if downscaled PNG already exists
if os.path.exists(pfile2):
continue
print "Creating \"{:s}\" ...".format(pfile2)
# Save downscaled PNG ...
subprocess.check_call(
["convert", pfile1, "-resize", "{:d}x".format(width), pfile2],
stderr = open(os.devnull, "wt"),
stdout = open(os.devnull, "wt")
)
subprocess.check_call(
["exiftool", "-overwrite_original", "-all=", pfile2],
stderr = open(os.devnull, "wt"),
stdout = open(os.devnull, "wt")
)
subprocess.check_call(
["optipng", pfile2],
stderr = open(os.devnull, "wt"),
stdout = open(os.devnull, "wt")
)
# Save JSON dictionary ...
open(
"images.json",
"wt",
).write(
json.dumps(
data,
indent = 4,
sort_keys = True
)
)
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.