generate.py

#!/usr/bin/env python3

import sys
import os
import re
import math
from contextlib import closing

scale = 3.0
splitter = re.compile(r'\s+')

if math.fmod(360.0, scale) != 0:
    raise RuntimeError('non-exact scale')
if math.fmod(180.0, scale) != 0:
    raise RuntimeError('non-exact scale')

bandsize = int(360.0 / scale + 1)
bandcount = int(180.0 / scale + 1)
tablesize = bandsize * bandcount
table = [None] * tablesize

print('{} degree grid, {} table entries'.format(scale, tablesize), file=sys.stderr)
print('Generating INPUT.DAT', file=sys.stderr)
with closing(open('INPUT.DAT', 'w')) as f:
    for lat in range(bandcount):
        for lon in range(bandsize):
            print(lat * scale - 90.0, lon * scale - 180.0, file=f)

print('Generating OUTPUT.DAT', file=sys.stderr)
os.system('./interpolate')

with closing(open('OUTPUT.DAT', 'r')) as f:
    for line in f:
        line = line.strip()
        lat, lon, geoid = splitter.split(line.rstrip('\r\n'))
        lat, lon, geoid = float(lat), float(lon), float(geoid)

        lat_i = int((lat + 90) / scale)
        lon_i = int((lon + 180) / scale)
        index = (lat_i * bandsize) + lon_i
        table[index] = geoid

for i in range(tablesize):
    if table[i] is None:
        raise RuntimeError('missing values at index {}'.format(i))

print('Generating geoid.rs', file=sys.stderr)
with closing(open('geoid.rs', 'w')) as f:
    prog = '''
// Generated by generate.py, DO NOT EDIT!

/// Grid size in degrees.
pub(crate) const SCALE: f32 = {scale};

/// How many entries each band has.
pub(crate) const BAND_SIZE: usize = {bandsize};

/// How many bands there are in total.
pub(crate) const BAND_COUNT: usize = {bandcount};

/// An array of geoid undulation data with BAND_SIZE * BAND_COUNT entries.
///
/// Data is expressed in meters.
#[allow(clippy::approx_constant)]
pub(crate) const DATA: &[f32] = &{table};
'''.format(scale=scale,
           bandsize=bandsize,
           bandcount=bandcount,
           table=table)
    print(prog, file=f)

print('Formatting geoid.rs', file=sys.stderr)
os.system('rustfmt geoid.rs')