bathymetry

README.md

A MATLAB translation of a function that reads OTPS GRID binary files.

read_bathy.py

#!/usr/bin/env python

import argparse
import struct
import numpy
from gaip import write_img
from eotools.geobox import GriddedGeoBox

def grd_in(fname):
    with open(fname, 'rb') as f:
        f.seek(4)
        n = struct.unpack('>i', f.read(4))[0]
        m = struct.unpack('>i', f.read(4))[0]
        lats = list(struct.unpack('>ff', f.read(8)))
        lons = list(struct.unpack('>ff', f.read(8)))
        dt = struct.unpack('>f', f.read(4))

        if (lons[0] < 0) & (lons[1] < 0) & (dt > 0):
            lons = [lon + 360 for lon in lons]

        ll_lims = [lons, lats]

        nob = struct.unpack('>i', f.read(4))[0]

        if nob == 0:
            f.seek(20, 1)
            iob = []
        else:
            f.seek(8, 1)
            bin_seq = 2 * nob 
            iob = struct.unpack('>{}i'.format(bin_seq), f.read(bin_seq * 4))
            iob = numpy.array(iob).reshape((2, nob))
            f.seek(8, 1)

        nm = n * m
        hz = struct.unpack('>{}f'.format(nm), f.read(nm * 4))
        hz = numpy.flipud(numpy.array(hz).reshape((m, n)))
        f.seek(8, 1)
        mz = struct.unpack('>{}i'.format(nm), f.read(nm * 4))
        mz = numpy.flipud(numpy.array(mz, dtype='int32').reshape((m, n)))

    return n, m, lats, lons, dt, ll_lims, nob, iob, hz, mz


if __name__ == '__main__':

    desc = 'Converts a OTPS GRID file to GeoTIFF'
    parser = argparse.ArgumentParser(description=desc)

    parser.add_argument('--gridfile', required=True,
                        help='The input GIRD file.')
    parser.add_argument('--outfname1', required=True,
                        help='The output filename')
    parser.add_argument('--outfname2', required=True,
                        help='The output filename')

    parsed_args = parser.parse_args()

    fname = parsed_args.gridfile
    outfname1 = parsed_args.outfname1
    outfname2 = parsed_args.outfname2

    x, y, _, _, _, extents, _, _, hz, mz = grd_in(fname)

    px = (extents[0][1] - extents[0][0]) / x
    py = (extents[1][1] - extents[1][0]) / y

    geobox = GriddedGeoBox(shape=(y, x), origin=(extents[0][0], extents[1][0]),
                           pixelsize=(px, py))
    write_img(hz, outfname1, fmt='GTiff', geobox=geobox)
    write_img(mz, outfname2, fmt='GTiff', geobox=geobox)