# This script shows how to read LGM2026 grid files "./lgm2026-grd-*.h5" in
# Python.
#
# To run the script, you must have installed the "h5py" and "numpy" libraries
# on your machine ("pip install h5py numpy") and your working directory must be
# the directory, where the "h5" files are stored.
#
# Note the casts from single to double precision before adding the offsets.
# This is necessary, as single precision does not have enough bits to store
# some of the datasets in full precision (e.g., the gravitational potential).
#
# The program should print an output like this:
#
# At the point with
#
#         latitude: "29.99609375 degree"
#         longitude: "0.00390625 degree"
#         spherical radius: "1735438.1171875 metre"
#
# the value of vz is
#
#         "-162745.94203948975 mGal".



import h5py
import numpy as np

# Open the stream for some tile
tile = h5py.File('./lgm2026-grd-N00E000.h5', 'r')

# Dataset to read: the z-element of the gravitational vector at the topography
dataset = '/gravitational-field/topography/vz'

# Read the dataset (note the cast to double precision)
vz = tile[dataset][()].astype(np.float64)

# Add the offset associated with the "tile" to "vz".  Note that we are casting
# the offset to double precision before adding it to "vz", which has already
# been cast to double precision.  This is the correct way to work with LGM2026
# data.
vz += tile[dataset].attrs['offset'].astype(np.float64)

# Get the unit of "vz"
vz_unit = tile[dataset].attrs['unit']

# Spherical latitudes.  No casting to double precision is needed here, as
# latitudes are stored in double precision.
dataset   = '/position/topography/latitude'
lat       = tile[dataset][()]
lat      += tile[dataset].attrs['offset']
lat_unit  = tile[dataset].attrs['unit']

# Spherical longitudes.  No casting to double precision is needed here, as
# longitudes are stored in double precision.
dataset   = '/position/topography/longitude'
lon       = tile[dataset][()]
lon      += tile[dataset].attrs['offset']
lon_unit  = tile[dataset].attrs['unit']

# Spherical radii
dataset = '/position/topography/radius'
r       = tile[dataset][()].astype(np.float64)
r      += tile[dataset].attrs['offset'].astype(np.float64)
r_unit  = tile[dataset].attrs['unit']

print(f'At the point with\n\n\tlatitude: \"{lat[0]} {lat_unit}\"'
      f'\n\tlongitude: \"{lon[0]} {lon_unit}\"'
      f'\n\tspherical radius: \"{r[0, 0]} {r_unit}\"'
      f'\n\nthe value of vz is\n\n\t\"{vz[0, 0]} {vz_unit}\".')

tile.close()