hainegroup / oceanspy

A Python package to facilitate ocean model data analysis and visualization.
https://oceanspy.readthedocs.io
MIT License
102 stars 32 forks source link

od.plot.vertical_section() doing something odd #265

Closed Mikejmnez closed 2 years ago

Mikejmnez commented 2 years ago

Description

I calculated a survey on the LLC4320 along A05 which worked fine. I then proceeded to plot its vertical section using oceanspy plotting functionality. This is how the section looked:

Screen Shot 2022-10-04 at 4 34 34 PM

The plot shows a spurious anomaly below 1000m uniformly along the section. I plotted a vertical section of Temperature around the mark 2000km and it didn't show anything wrong with the data.

To be sure after I stored the survey locally (created a checkpoint) I plotted it using xarray's basic plotting functionality as described below (od_survA05 is the same survey):

fig= plt.figure(figsize=(15, 6), facecolor='w')
od_survA05.dataset['Temp'].isel(time=0).plot(vmin=1, vmax=25, cmap='RdBu_r')
plt.show()

And the result was:

Screen Shot 2022-10-04 at 4 34 56 PM

I don't know what the issue is, but I have seen this behavior once before when plotting data along a single face (no transformation).

MaceKuailv commented 2 years ago

This is not a bug in plotting function. When oceanspy is doing the survey at time '2012-04-25T00', it also included that of '2012-04-25T01:00', and the data is apparently corrupted. image What we should fix is the data we stored.

ThomasHaine commented 2 years ago

Good sleuthing! @Mikejmnez can we write the data 2012-04-25T01:00 again (e.g., reconvert from binary to zarr, download again)?

Mikejmnez commented 2 years ago

Thanks @MaceKuailv I just confirmed this too. I was selecting, after the survey for a single snapshot,

.isel(time=0)

which does not show the temperature corrupted data.

If I do

.isel(time=1)

the corrupted data does shows up. Weird that when calculating the survey for a single snapshot of time, the end result is two-snapshots in time.

Here is the side-by-side comparison between the time=0 and time=1 indexes.

Screen Shot 2022-10-24 at 2 55 21 PM

But that also explains why I couldn't find the corrupted data when looking directly at the original dataset at time=0

Mikejmnez commented 2 years ago

Yes, we do have the binary file for that snapshot

Mikejmnez commented 2 years ago

I will close this issue then...

Mikejmnez commented 2 years ago

@ThomasHaine

I am re opening this issue just to share what I found by looking at the original (binary) file. I can see the same corrupted data.

temp_face2_from_binary

The file is here:

file = '/sciserver/filedb12-03/ocean/poseidon/llc4320_tests/10dayhourly/Theta.0000788112.data.shrunk'

This is the code I used to plot the variable:

import xarray as xr
import dask
from xmitgcm import llcreader
from fsspec.implementations.local import LocalFileSystem
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
plt.rcParams['figure.figsize']=12,8

fs = LocalFileSystem()

iters = [787968, 788112, 788256, 788400, 788544, 788688, 788832, 788976]

mask_path = '/sciserver/filedb01-01/ocean/poseidon/llc4320_masks.zarr/'
store_u2 = llcreader.BaseStore(fs, base_path='/sciserver/filedb12-03/ocean/poseidon/llc4320_tests/10dayhourly',
                               shrunk=True, mask_fs = fs, mask_path=mask_path, join_char='/')
ds_u2 = llcreader.LLC4320Model(store_u2).get_dataset(varnames=['Theta'], iter_start=iters[1], iter_stop=iters[1]+1)
ds_u2 = ds_u2.reset_coords()
ds_u2 = ds_u2.drop_vars([var for var in ds_u2.data_vars if var not in ['Theta']])

Temp = ds_u2['Theta'].isel(time=0, face=2, j = 3600)
fig= plt.figure(figsize=(12, 8), frameon=True, facecolor='w')
with dask.config.set(scheduler='single-threaded'):
    Temp.isel(i=slice(0, 2000, 2)).plot(vmin=-1, vmax=10, cmap='RdBu_r')

plt.savefig('temp_face2_from_binary.png')

I ran it in datascope, using one of the elephant nodes and within the Oceanography environment

Mikejmnez commented 2 years ago

When using xmitgcm version 0.5.2 it may require to specify the grid path. Here is the code that I use for that version of xmitgcm:

import xarray as xr
import dask
from xmitgcm import llcreader
from fsspec.implementations.local import LocalFileSystem
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
plt.rcParams['figure.figsize']=12,8

fs = LocalFileSystem()

iters = [787968, 788112, 788256, 788400, 788544, 788688, 788832, 788976]

mask_path = '/sciserver/filedb01-01/ocean/poseidon/llc4320_masks.zarr/'
grid_path = '/sciserver/filedb01-01/ocean/poseidon/llc4320_grid/'
store_u2 = llcreader.BaseStore(fs, base_path='/sciserver/filedb12-03/ocean/poseidon/llc4320_tests/10dayhourly',
                               shrunk=True, grid_path = grid_path, mask_fs = fs, mask_path=mask_path, join_char='/')
ds_u2 = llcreader.LLC4320Model(store_u2).get_dataset(varnames=['Theta'], iter_start=iters[1], iter_stop=iters[1]+1)
ds_u2 = ds_u2.reset_coords()
ds_u2 = ds_u2.drop_vars([var for var in ds_u2.data_vars if var not in ['Theta']])

Temp = ds_u2['Theta'].isel(time=0, face=2, j = 3600)
fig= plt.figure(figsize=(12, 8), frameon=True, facecolor='w')
with dask.config.set(scheduler='single-threaded'):
    Temp.isel(i=slice(0, 2000, 2)).plot(vmin=-1, vmax=10, cmap='RdBu_r')

plt.savefig('temp_face2_from_binary.png')
Mikejmnez commented 2 years ago

Closing as this is not an oceanspy issue