Created
November 13, 2023 23:30
-
-
Save mx-moth/82b56b7b6452bae5d6b7e34048acee01 to your computer and use it in GitHub Desktop.
Extract a point with emsarray
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
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() |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment