# This file shows how to read LGM2026 in Python.  The "h5py" and "numpy"
# libraries must be installed on your machine ("pip install h5py numpy") and
# your working directory must be the "software" directory.
#
# Note the casts from single to double precision before adding the offsets.
# This is necessary, as single precision does not have enough digits to store
# some of the datasets in their full accuracy (e.g., the gravitational
# potential).
#
# The program should print an output like this:
#
# At a 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('../grids/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 to "vz" that is associated with our "tile".  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 a 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()