Extract a point with emsarray

extract-point.py

import emsarray
from matplotlib import pyplot as plt
import numpy
import shapely

# Open a dataset
# This dataset starts with four dimensions: time, depth, and two surface dimensions
dataset = emsarray.open_dataset("https://dapds00.nci.org.au/thredds/dodsC/fx3/GBR1_H2p0_B3p2_Cfur_Dnrt.ncml")

# We can speed up future operations by limiting the timesteps
# and selecting only the variables we are interested in.
# This step is optional, and is here only to make the demonstration speedy.
dataset = dataset\
    .ems.select_variables(['temp', 'eta', 'botz'])\
    .isel(time=slice(-10, None))

# Extract a single point from the dataset
print("Selecting a single point...")
hervey_bay = shapely.Point(152.888677, -25.020682)
hervey_bay_dataset = dataset.ems.select_point(hervey_bay)

# Because we only selected a small subset of data,
# we can fetch it all now and speed up future computations.
# Don't do this on the whole aggregate endpoint
# unless you want to download many terabytes!
print("Fetching timeseries data for Hervey Bay...")
hervey_bay_dataset.load()
hervey_bay_dataset.to_netcdf("./hervey-bay.nc")

# This new point dataset now only has time and depth dimensions
print("Extracted point data for Hervey Bay:", hervey_bay_dataset)

# 'botz' is zero-dimensional, i.e. it is just a single value
print("Depth:", hervey_bay_dataset['botz'].values)

# 'eta' is one-dimensional, tracking the surface height of this point over time
eta = hervey_bay_dataset['eta']
print("Surface height mean:", numpy.mean(eta.values), "std:", numpy.std(eta.values))

# 'temp' is two-dimensional, as it also has depth.
temp = hervey_bay_dataset['temp']
# Hervey Bay is shallow. Most of the deep cells will be below the ocean floor.
# These cells will be filled with 'nan'.
# This will slice off all the data that is nan leaving only cells with values
deepest_index = numpy.where(numpy.isnan(temp.isel(time=0).values))[0][-1] + 1
temp = temp.isel(k=slice(deepest_index, None))
# Find the average temperature across time for each depth
mean_temp = temp.mean('time')
# Plot the mean temperature
fig = plt.figure()
plt.title("Mean temperature")
plt.ylabel("Depth (m)")
plt.xlabel("Temperature  (°C)")
plt.plot(mean_temp.values, mean_temp.coords['zc'].values)
plt.show()