From ba78eb9cec2b90c85d77fbac46a329a5c23c3b4c Mon Sep 17 00:00:00 2001 From: Juan Nunez-Iglesias Date: Wed, 27 Nov 2024 16:24:37 +1100 Subject: [PATCH 1/6] Add xarray atmospheric data example --- examples/xarray-climate.py | 99 ++++++++++++++++++++++++++++++++++++++ 1 file changed, 99 insertions(+) create mode 100644 examples/xarray-climate.py diff --git a/examples/xarray-climate.py b/examples/xarray-climate.py new file mode 100644 index 00000000000..3cf96eccab2 --- /dev/null +++ b/examples/xarray-climate.py @@ -0,0 +1,99 @@ +"""Viewing thredds wr45/ops_aps3 data in napari. + +Download bash commands: + +# download model prediction data +curl -O https://thredds.nci.org.au/thredds/fileServer/wr45/ops_aps3/access-g/1/20241104/0000/fc/ml/air_temp.nc +curl -O https://thredds.nci.org.au/thredds/fileServer/wr45/ops_aps3/access-g/1/20241104/0000/fc/ml/spec_hum.nc + +# download corresponding 10 days' worth of measurements +mkdir an && cd an # use 'an' folder for single time points +for day in 04 05 06 07 08 09 10 11 12 13; do + for hour in 00 06 12 18; do + curl https://thredds.nci.org.au/thredds/fileServer/wr45/ops_aps3/access-g/1/202411${day}/${hour}00/an/ml/spec_hum.nc -o ${day}-${hour}-spec_hum.nc + curl https://thredds.nci.org.au/thredds/fileServer/wr45/ops_aps3/access-g/1/202411${day}/${hour}00/an/ml/air_temp.nc -o ${day}-${hour}-air_temp.nc + done +done +""" + +from glob import glob +from pathlib import Path + +import numpy as np +import xarray as xr + +import napari + +root_dir = Path('/Users/jni/data/thredds/20241104') + +def get_scale_translate(dataset, array_name, invert_lat=False): + """Get the translate/offset and scale parameters for an xarray dataset. + + This code assumes that the dataset is regularly spaced. You should + interpolate your data if it is sampled at irregular spaces. + + Parameters + ---------- + dataset : xr.Dataset + The dataset containing the array to be displayed. + array_name : str + The name of the xarray DataArray within `dataset` to be displayed in + napari. + invert_lat : bool + Whether to invert the latitude values. + + napari's axes follow zyx convention, with y pointing down. + + latitude decreases going down conventionally (putting North at the top + and having southern latitudes be negative), so here we multiply by -1 + to display the globe conventionally. Unfortunately, this means the + coordinates displayed by the viewer on hover will show northern + latitudes as negative. The true fix is to add a transformation + between the world/scene coordinates and the canvas in the napari code + base, but that might take some time. + """ + array = getattr(dataset, array_name) + if array is None: + raise ValueError(f'{dataset} has no array with name {array_name}') + dims = [getattr(dataset, dim) for dim in array.dims] + translate = [float(d[0]) for d in dims] + scale = [float(d[1] - d[0]) for d in dims] + if invert_lat: + lat_pos, lat_name = next( + (i, dim) + for i, dim in enumerate(array.dims) + if dim.startswith('lat') + ) + scale[lat_pos] *= -1 + translate[lat_pos] = float(getattr(dataset, lat_name)[-1]) + return {'scale': scale, 'translate': translate} + + +an = xr.open_mfdataset(sorted(glob(str(root_dir / 'an/*spec_hum.nc')))) +ds = xr.open_dataset(root_dir / 'spec_hum.nc', chunks={'time': 1}) +start, stop, step = [ + np.array(elem)[()] + for elem in (ds.time[0], ds.time[-1], ds.time[1] - ds.time[0]) + ] +ds_reg = ds.interp( + coords={'time': np.arange(start, stop, step)}, + method='nearest', + ) + +viewer = napari.Viewer() +model = viewer.add_image( + ds_reg.spec_hum, + name='model', + **get_scale_translate(ds_reg, 'spec_hum', invert_lat=True), + colormap='magenta', + ) +gt = viewer.add_image( + an.spec_hum, + name='measurement', + **get_scale_translate(an, 'spec_hum', invert_lat=True), + colormap='green', + blending='additive', + ) + +if __name__ == '__main__': + napari.run() From ebe47a774df9b4879ab63cb63a7d5ac6d0051337 Mon Sep 17 00:00:00 2001 From: Juan Nunez-Iglesias Date: Thu, 28 Nov 2024 18:12:59 +1100 Subject: [PATCH 2/6] Add raw data viewer; add axis labels --- examples/xarray-climate.py | 47 ++++++++++++++++++++++++++------------ 1 file changed, 32 insertions(+), 15 deletions(-) diff --git a/examples/xarray-climate.py b/examples/xarray-climate.py index 3cf96eccab2..8000829ef78 100644 --- a/examples/xarray-climate.py +++ b/examples/xarray-climate.py @@ -69,8 +69,21 @@ def get_scale_translate(dataset, array_name, invert_lat=False): return {'scale': scale, 'translate': translate} -an = xr.open_mfdataset(sorted(glob(str(root_dir / 'an/*spec_hum.nc')))) +# open the model dataset ds = xr.open_dataset(root_dir / 'spec_hum.nc', chunks={'time': 1}) + +# Show the raw (not resampled) model data +viewer, model_layer = napari.imshow( + ds.spec_hum, + name='model', + **get_scale_translate(ds, 'spec_hum', invert_lat=True), +) +viewer.dims.axis_labels = ds.spec_hum.dims + +# open the measurement data +an = xr.open_mfdataset(sorted(glob(str(root_dir / 'an/*spec_hum.nc')))) + +# resample the model data to have even spacing in time start, stop, step = [ np.array(elem)[()] for elem in (ds.time[0], ds.time[-1], ds.time[1] - ds.time[0]) @@ -79,21 +92,25 @@ def get_scale_translate(dataset, array_name, invert_lat=False): coords={'time': np.arange(start, stop, step)}, method='nearest', ) +an = xr.open_mfdataset(sorted(glob(str(root_dir / 'an/*spec_hum.nc')))) -viewer = napari.Viewer() -model = viewer.add_image( - ds_reg.spec_hum, - name='model', - **get_scale_translate(ds_reg, 'spec_hum', invert_lat=True), - colormap='magenta', - ) -gt = viewer.add_image( - an.spec_hum, - name='measurement', - **get_scale_translate(an, 'spec_hum', invert_lat=True), - colormap='green', - blending='additive', - ) +# show the resampled model data overlaid on the measurement data +# note: this has performance issues due to the live resampling. +viewer2 = napari.Viewer() +model = viewer2.add_image( + ds_reg.spec_hum, + name='model', + **get_scale_translate(ds_reg, 'spec_hum', invert_lat=True), + colormap='magenta', +) +gt = viewer2.add_image( + an.spec_hum, + name='measurement', + **get_scale_translate(an, 'spec_hum', invert_lat=True), + colormap='green', + blending='additive', +) +viewer2.dims.axis_labels = ds_reg.spec_hum.dims if __name__ == '__main__': napari.run() From 93ec6cdb0e9c5a7abd31a9efc47371c4fbdbb4c8 Mon Sep 17 00:00:00 2001 From: Juan Nunez-Iglesias Date: Fri, 29 Nov 2024 00:37:00 +1100 Subject: [PATCH 3/6] =?UTF-8?q?don't=20try=20to=20open=20the=20same=20data?= =?UTF-8?q?set=20twice=20=E2=80=94=20this=20crashes=20xarray?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- examples/xarray-climate.py | 1 - 1 file changed, 1 deletion(-) diff --git a/examples/xarray-climate.py b/examples/xarray-climate.py index 8000829ef78..946f7ea2774 100644 --- a/examples/xarray-climate.py +++ b/examples/xarray-climate.py @@ -92,7 +92,6 @@ def get_scale_translate(dataset, array_name, invert_lat=False): coords={'time': np.arange(start, stop, step)}, method='nearest', ) -an = xr.open_mfdataset(sorted(glob(str(root_dir / 'an/*spec_hum.nc')))) # show the resampled model data overlaid on the measurement data # note: this has performance issues due to the live resampling. From f49504e688d36857b910ac8e649aec5972101dac Mon Sep 17 00:00:00 2001 From: Juan Nunez-Iglesias Date: Fri, 29 Nov 2024 19:17:58 +1100 Subject: [PATCH 4/6] formatting and use private camera to filp y --- examples/xarray-climate.py | 42 +++++++++++++++++++++++--------------- 1 file changed, 26 insertions(+), 16 deletions(-) diff --git a/examples/xarray-climate.py b/examples/xarray-climate.py index 946f7ea2774..1b5aab22e53 100644 --- a/examples/xarray-climate.py +++ b/examples/xarray-climate.py @@ -50,7 +50,11 @@ def get_scale_translate(dataset, array_name, invert_lat=False): coordinates displayed by the viewer on hover will show northern latitudes as negative. The true fix is to add a transformation between the world/scene coordinates and the canvas in the napari code - base, but that might take some time. + base, but that might take some time. An alternative for now is to use + private APIs to find the VisPy camera and change its orientation, like + so: + + viewer.window._qt_window._qt_viewer.canvas.camera._2D_camera.flip = (0, 0, 0) """ array = getattr(dataset, array_name) if array is None: @@ -74,11 +78,14 @@ def get_scale_translate(dataset, array_name, invert_lat=False): # Show the raw (not resampled) model data viewer, model_layer = napari.imshow( - ds.spec_hum, - name='model', - **get_scale_translate(ds, 'spec_hum', invert_lat=True), -) + ds.spec_hum, + name='model', + **get_scale_translate(ds, 'spec_hum'), + ) viewer.dims.axis_labels = ds.spec_hum.dims +# currently no private API to flip the camera, so increasing y is up, +# so we use these private attributes. +viewer.window._qt_window._qt_viewer.canvas.camera._2D_camera.flip = (0, 0, 0) # open the measurement data an = xr.open_mfdataset(sorted(glob(str(root_dir / 'an/*spec_hum.nc')))) @@ -97,19 +104,22 @@ def get_scale_translate(dataset, array_name, invert_lat=False): # note: this has performance issues due to the live resampling. viewer2 = napari.Viewer() model = viewer2.add_image( - ds_reg.spec_hum, - name='model', - **get_scale_translate(ds_reg, 'spec_hum', invert_lat=True), - colormap='magenta', -) + ds_reg.spec_hum, + name='model', + **get_scale_translate(ds_reg, 'spec_hum'), + colormap='magenta', + ) gt = viewer2.add_image( - an.spec_hum, - name='measurement', - **get_scale_translate(an, 'spec_hum', invert_lat=True), - colormap='green', - blending='additive', -) + an.spec_hum, + name='measurement', + **get_scale_translate(an, 'spec_hum'), + colormap='green', + blending='additive', + ) viewer2.dims.axis_labels = ds_reg.spec_hum.dims +# currently no private API to flip the camera, so increasing y is up, +# so we use these private attributes. +viewer2.window._qt_window._qt_viewer.canvas.camera._2D_camera.flip = (0, 0, 0) if __name__ == '__main__': napari.run() From 7d3396e518b601227cfbd4f3111a725182deea28 Mon Sep 17 00:00:00 2001 From: Juan Nunez-Iglesias Date: Mon, 9 Dec 2024 16:32:05 +1100 Subject: [PATCH 5/6] Improve chunking based on h5 chunks --- examples/xarray-climate.py | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/examples/xarray-climate.py b/examples/xarray-climate.py index 1b5aab22e53..747c02111d6 100644 --- a/examples/xarray-climate.py +++ b/examples/xarray-climate.py @@ -74,7 +74,10 @@ def get_scale_translate(dataset, array_name, invert_lat=False): # open the model dataset -ds = xr.open_dataset(root_dir / 'spec_hum.nc', chunks={'time': 1}) +ds = xr.open_dataset( + root_dir / 'spec_hum.nc', + chunks={'time': 1, 'theta_lvl': 1}, + ) # Show the raw (not resampled) model data viewer, model_layer = napari.imshow( From 09889cd11f44c8bc0a411eaffcc96c265f793236 Mon Sep 17 00:00:00 2001 From: Juan Nunez-Iglesias Date: Mon, 9 Dec 2024 16:32:26 +1100 Subject: [PATCH 6/6] =?UTF-8?q?Use=20assume=5Fsorted=20in=20interp=20?= =?UTF-8?q?=E2=80=94=20helps=20speed?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- examples/xarray-climate.py | 1 + 1 file changed, 1 insertion(+) diff --git a/examples/xarray-climate.py b/examples/xarray-climate.py index 747c02111d6..88a0dfc1b34 100644 --- a/examples/xarray-climate.py +++ b/examples/xarray-climate.py @@ -101,6 +101,7 @@ def get_scale_translate(dataset, array_name, invert_lat=False): ds_reg = ds.interp( coords={'time': np.arange(start, stop, step)}, method='nearest', + assume_sorted=True, ) # show the resampled model data overlaid on the measurement data