Python-kode som les inn bin-filer

# -*- coding: utf-8 -*-
'''
Under er eit døme på å lese inn filer av det binære .bin-formatet som er brukt
for separasjonsmodellane mellom referansenivå i kystsona og til havs, publisert
av Kartverket. Dette er same format som er brukt for høgdereferansemodellane
for NN2000 og NN1954.

Det er definert ein funksjon for å lese inn filene til ein dictionary,
deretter vert funksjonen brukt for å lese inn fila
"ChartDatum_above_Ellipsoid_EUREF89_v2021a.bin", og resultatet plotta.

I koden under er det antatt at den .bin-fila ligg i same mappe som skriptet.

Lisensiert under Creative Commons BY-SA 3.0
https://creativecommons.org/licenses/by-sa/3.0/
'''

import struct
import numpy as np
import matplotlib.pyplot as plt


def readbinfile(filename):
    '''Returnerer dictionary med data frå .bin-fil.

    Parametre
    ---------
    filename : string
        Filnamn eller fullstendig sti til .bin-fil som skal lesast inn

    Returnerer
    ----------
    D : dictionary
        Oppføringer i dict er
         - lat : Nx1 numpy array med breiddegradsverdier
         - lon : Mx1 numpy array med lengdegradsverdier
         - data : NxM numpy masked array med modellverdiar
         - dlon : avstand mellom gridpunkt i aust-vest-retning
         - dlat : avstand mellom gridpunkt i nord-sør-retning
         - ellipsoide : viss UTM-koordinatar, gjev aktuell ellipsoide
         - UTM-sone : viss UTM-koordinatar, gjev UTM-sonen

    Notat
    -----
    Formatet i .bin-filene som HREF og separasjonsmodellane kjem i er eit
    Gravsoft-binærformat. Verdiane er på eit regulært grid, med radene langs
    linjer av breiddegrad, og kolonnene langs linjer av lengdegrad.
    Fila er delt i blokker på 64 byte, den fyrste blokken inneheld metadata.

    Resten av fila er modellverdiar som "single precision" flyttal, det vil
    seie 4 byte for kvart tal.

    Alle tal i same 64 byte-blokk er på same breiddegrad, og verdiane går
    frå vest mot aust. Fyrste rad er den nordlegaste. Det vil seie, den fyrste
    verdien er i det nordvestre hjørnet, den siste i det søraustre hjørnet.

    Innhald i metadata-blokken, fyrste 64 byte:
    - Byte 1-4: Icode - 4 byte integer, variabel som skal ha verdi 777.
        Viss ikkje er det feil med big/little endian.
    - Byte 5-12: Rnmin - 8 byte float/double, gir minste breiddegrad for grid
    - Byte 13-20: Rnmax - 8 byte float/double, gir største breiddegrad for grid
    - Byte 21-28: Remin - 8 byte float/double, gir minste lengdegrad for grid
    - Byte 29-36: Remin - 8 byte float/double, gir største lengdegrad for grid
    - Byte 37-44: deltaN - 8 byte float/double, gir gridstørrelse i
        breiddegradsretning
    - Byte 45-52: deltaE - 8 byte float/double, gir gridstørrelse i
        lengdegradsretning
    - Byte 53-56: Iutm - 4 byte integer. Viss geografisk grid lik 0
    - Byte 57-60: Iell - 4 byte integer. Lik 0 for geografisk grid.
        For UTM-grid, kode for å definere ellipsoide.
    - Byte 61-64: Izone - 4 byte integer. Angir UTM-sone,
        viss geografisk grid er den lik 0.

    Griddet i fila vil då ha N = (Rnmax - Rnmin)/deltaN + 1 rader og
    M = (Remax - Remin)/deltaE + 1 kolonner.

    OBS: viss M ikkje går opp i 16 i talet på verdiar i kvar blokk, vil slutten av
    den siste blokken på ei rad vere fylt med null. T.d. om M = 2801, vil det
    vere mod(2801, 16) = 1 reelle datapunkt i den siste blokken på ei rad,
    og dermed får ein 16 - 1 = 15 nullverdiar på slutten av kvar rad.

    Metode for lesing tatt frå  https://stackoverflow.com/a/8711061
    '''

    # les inn heile fila gitt ved stien
    with open(filename, mode='rb') as file:
        fileContent = file.read()

    # pakk ut fyrste fire byte til int, og sjekk om lik 777
    Icode = struct.unpack('i', fileContent[:4])
    if Icode[0] != 777:
        raise('Icode er ikkje 777, feil endian')

    # pakk ut neste 48 byte som double precision floats
    Rnmin, Rnmax, Remin, Remax, dn, de = struct.unpack('dddddd',
                                                       fileContent[4:4+8*6])

    # pakk ut siste del av metadatablokk som tre ints
    Iutm, Iell, Izone = struct.unpack('iii', fileContent[4+8*6:64])

    # dict med ellipsoider
    ellipsoids = {
        0: 'LatLon',
        1: 'WGS84/NAD83',
        2: 'Hayford/ED50',
        3: 'Clarke/NAD27',
        }

    modEll = ellipsoids[Iell] if Iell in ellipsoids.keys() else 'Undefined'

    # har ikkje testa med UTM-grid, so gi ut advarsel viss det er tilfelle
    if Iutm != 0:
        Warning('UTM-koordinatar i fil -- veit ikkje om dette fungerer')

    # pakk ut resten av fila som single precision float, og legg i numpy array
    data = struct.unpack('f' * ((len(fileContent) - 64) // 4),
                         fileContent[64:])
    data = np.array(data)

    # rekn ut storleik på grid
    Nlon = round((Remax-Remin)/de) + 1
    Nlat = round((Rnmax-Rnmin)/dn) + 1

    # lag dict med lat/lon vektorar og anna metadata
    # merk lat-vektoren går frå nord til sør for å samsvare med data-matrisa
    D = {'lat': np.linspace(Rnmax, Rnmin, Nlat),
         'lon': np.linspace(Remin, Remax, Nlon),
         'dlat': dn,
         'dlon': de,
         'ellipsoide': modEll,
         'UTM-sone': Izone}

    N = data.size

    if (N / Nlat) % 1 != 0:
        raise('Talet på datapunkt går ikkje opp i talet med breiddegradpunkt.')

    # kor mange datapunkt det er per rad, inkludert nullar
    Nlon_with_zeros = N//Nlat

    # (differansen mellom matrisestørrelse med og utan null-verdiar) delt på
    # (differansen mellom radstørrelser med og utan null-verdiar)
    # burde verte lik talet på rader
    if (N - Nlon*Nlat) / (Nlon_with_zeros - Nlon) != Nlat:
        raise('Feil tal på datapunkt i forhold til lat-lon-griddet.')

    # reshape data til array, og fjern nullverdiane på slutten av kvar rad
    data = data.reshape((int(Nlat), int(Nlon_with_zeros)))[:, :Nlon]

    # lag masked_array og legg til i dict
    # maske basert på størrelsen til verdiane
    D['data'] = np.ma.masked_array(data, mask=np.abs(data) > 500)

    return D


# definer namn på fila som skal lesast inn, og bruk funksjonen over
filnamn = 'ChartDatum_above_Ellipsoid_EUREF89_v2021a.bin'
modeldata = readbinfile(filnamn)

# lag matriser med lengde- og breiddegradsverdiar
lon, lat = np.meshgrid(modeldata['lon'], modeldata['lat'])

# plott resultatet
fig, ax = plt.subplots()
g = ax.pcolormesh(lon, lat, modeldata['data'])
ax.set_xlabel('Lengdegrad')
ax.set_ylabel('Breiddegrad')
ax.set_title('Sjøkartnull over ellipsoiden')
cb = fig.colorbar(g)
cb.set_label('Meter')
plt.show()