writing multidimensional data to netcdf for space-varying z coordinate, usign xArray

[Marco Miani]

Dear Group, 

my name is Marco and I am new here - hope this email finds you all well. 

I am having some hard time in trying to convert my model-generated data into netcdf, so to have it plotted in PRAVIEW (which would allow better visualization). My apologies for this lengthy message, but it's difficult to explain the complexity of the problem. 

Let me explain the problem: 

1. the original dataset is a matfile, generated by some numerical model, and contains the concentration of a tracer (salt), run for an academic test case in a small wave flume (2 m long, 30 cm deep). The model run is 1D (so no y dimension). 

2. All output variables, static (bathymetry, and horizontal coordinate) and  time-varying (remaining hydrodynamic variables), are condensed into that same matlfile. I wrote a python script to extract the data I want/need. 

3. Now, take a deep breath. The structure of the data, in my  specific case of salinity, is level-dependent, AND time-dependent. This means that if I run the model with, say, 130 vertical layers, I will have 130 variables in my matfile, containing the numerical value of the concentration. For instance: Salk100_0000_00. 

• Sal stands for salinity
• k100 stands for the number of layer
• 000000_000 is the time-step in HHMMSS_mms (thus here, initial condition)

Now this variable contains the numerical calculated salt concentration, and not where it is to be located in the computational domain. 

4. The vertical elevation of each layer interface, which is indeed time-varying,  is stored into another, separate, variable in the same file: zk100_0000_00

5. Unlike the vertical coordinate, the horizontal coordinate is always constant over time and space, i.e. it is fixed, it does not move. The vertical elevation of the layer interface, instead, may actually fluctuate, to follow the flow, or to adapt to the terrain (or bottom). 

6. In Python, I have no problem to lump everything together: to plot-contour data, I need data set and coordinate matrices to be consistent in size. 

• Therefore, Xp will be stacked (concatenated vertically) so to have the same size of Salk100_000_000. This is easy, I am just copying Xp and piling it up. 
• zk_100_0000_000 is varying. So, for that timestep, I need to loop over all k levels and stack them together. for k in range(0,100): new_zk=.....

No problem in doing that in python. I can plot the initial condition of my model and looks good, but when trying to convert my dataset into a netcdf with xarray, preserving the structure, is where things go ugly. 

The numerical model is run for a simple case over an up-sloping ramp. The settings are: Lx = 2m, Dz = 30cm, two different body of waters, fresh and salty, in initial unperturbed contact meet in the middle of the domain and mixing may be observed.

So I know that flow occurs on a sloping plane! And I can surely see that when I plot my results in python (when tapping into model-generated data), but when I save my data to netcdf, I loose that vertical structure and every thing looks perfectly horizontal. 

So, finally, my question: 

when prescribing X coordinate in my dataset, I impose 

qname= "Salk"

ds = xr.Dataset(

    {qname: (("z", "y","x"), vals)},

    coords={

        "x": Xp.squeeze(),

        "y": Y,

        "z": zk[:,0],    #----- I only take first vertical definition - erroneously applied everywhere  

    }

)

In essence, each column of z, should have its individual value. Here, instead, I take one value and I apply it all over the domain. 

I want to be able to include the whole vertical coordinate matrix zk into my file. 

Can I do that? 
Or should I include another dimension when preallocating data?  

I enclose two images in this email (one model generated, top; and one obtained in Paraview, after transforming my data, bottom). 

I am also sharing my so far written Python code for pre/post-post processing, as well as the model-generated data (zipped file attached to this email) 

My sincere gratitude to whoever might be willing and able to help with with some constructive comments, hints, inputs. 

Have a happy and sunny day 
Marco 

I 

[Ryan Abernathey]

Marco,

I'll try to give a brief answer to your question. Your problem with vertical coordinates is similar to the ROMS ocean model example in the xarray docs:

http://xarray.pydata.org/en/stable/examples/ROMS_ocean_model.html

I would try to represent your domain's geometry using non-dimension coordinates:

http://xarray.pydata.org/en/stable/data-structures.html#coordinates

I know this is not a complete solution, but hopefully this gives you some ideas to start from.

-Ryan

--
You received this message because you are subscribed to the Google Groups "xarray" group.
To unsubscribe from this group and stop receiving emails from it, send an email to xarray+un...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/xarray/2db0fdc1-990e-4545-b218-b5d32c74fcafn%40googlegroups.com.

[Marco Miani]

Hy Ryan, 

Thank you very much for your message. 

I had a look at the link, I confess I was unaware of that ROMS example.

Your suggestion is definitely a good start. Now that you raise this, do you know if for ROMS numbers of interfaces and layers is the same? In SWASH (which is I am trying to work on, right now) nrLayers = nrInterfaces-1

When preparing data for plotting, I find myself "cheating", by adding one fake (copied) extra layer to make dimensions consistent....

However, concerning the original problem, someone suggested the following strategy: 

1. foreach layer: read in quantity at that given layer, along with its belonging (spatial varying) vertical coordinate matrix

2. create a DataArray (and not a Dataset), where physical quantity and coordinate(s) are stored. 

3. At the end of the loop, merge all together and generate a DataArray.

By so doing, each Dataarray (i.e., each level) might have its own definition of z. 

I was doing something similar in Python then, after exporting it in netCDF, ingesting that into PARAVIEW - and the rest is history. 

Now, with the method described above, I want to see how far I can get. 

Have a great weekend everyone. 
Marco 

To view this discussion on the web visit https://groups.google.com/d/msgid/xarray/CAJAZx5A9JAFbLzFWRMue_ZwzGpot6V0sSHOsAFvP4A%2B7WgJQJw%40mail.gmail.com.

[Marco Miani]

Dear all, 

I am getting closer to the solution for my problem. 
Does anyone have any idea of how to generate, in python/x-array, the netcdf file (or a similar equivalent) to the one presented in this following example, shown on the official webpage? 

https://xesmf.readthedocs.io/en/v0.1.2/Curvilinear_grid.html

ds = xr.tutorial.open_dataset('rasm') ##---rasm.nc! 

If I manage to generate an equivalent empty dummy with the same structure, I can stuff it with the data I am interested in visualizing. 

Thanks! 

M.

--

[Deepak Cherian]

Take a look at xr.full_like, xr.zeros_like, xr.ones_like but it will be hard to change sizes.

The first example here: https://xarray.pydata.org/en/stable/data-structures.html#creating-a-dataset, seems like it creates what you want.

Deepak

To view this discussion on the web visit https://groups.google.com/d/msgid/xarray/CABjf4Fi8dV8wZ0f-05J7r_zwctnmtaPLLu39_YWgrM84x8YhHw%40mail.gmail.com.

[Marco Miani]

Hi Deepack, hi all

Thank you for your kind input. I have tried hard to get something out of it, but to no avail. 

I am trying to generate the same file (as in the example you had linked). Structure is the same, variable names are the same, metadata are all there. But still it doesn't seem to work, when visualizing in paraview. 

I am enclosing both files, mine and xarray's, maybe there is some kind sould out there that can help me spotting where the problem is? 

From xarray: 

xr.tutorial.open_dataset('rasm')

The one I have generated: 

test_from_Scratch.py (attached)

Resulting netcdf file generated with my python script is attached too...

All my gratitude and appreciation for any possible help you might be willing to give. 

cheers.

M

To view this discussion on the web visit https://groups.google.com/d/msgid/xarray/9a9b96c5-b92d-43ab-247c-2c8e353c59da%40gmail.com.

[Deepak Cherian]

On 10/21/20 8:51 AM, Marco Miani wrote:

But still it doesn't seem to work, when visualizing in paraview.

What does this mean? Paraview doesn't recognize the logical coordinates? You may need to specify .attrs["coordinates"] appropriately.

[Marco Miani]

Hi Deepak. 

It means that PARAVIEW displays it on an x,y plane, instead of a x,z plane. And that is not what I was expecting, since I (think I) have prescribed it the right way. Structure, names & number of variables are all consistent with the netcdf file of the example. I feel I am very close, but not quite there yet. 

When doing a ncdump, I get: 

ncdump -h from_scratch.nc 

netcdf from_scratch {

dimensions:

time = 100 ;

y = 12 ;

x = 500 ;

variables:

double salinity(time, y, x) ;

salinity:_FillValue = NaN ;

salinity:coordinates = "yc xc" ;

int64 time(time) ;

time:units = "days since 2014-09-06 00:00:00" ;

time:calendar = "proleptic_gregorian" ;

double xc(y, x) ;

xc:_FillValue = NaN ;

xc:axis = "X" ;

xc:long_name = "x-coordinate in Cartesian system" ;

xc:units = "m" ;

double yc(y, x) ;

yc:_FillValue = NaN ;

yc:axis = "Z" ;

yc:long_name = "sigma" ;

yc:positive = "down" ;

yc:units = "m" ;

// global attributes:

:convention = "CF-1.4" ;

 

--

You received this message because you are subscribed to the Google Groups "xarray" group.
To unsubscribe from this group and stop receiving emails from it, send an email to xarray+un...@googlegroups.com.

To view this discussion on the web visit https://groups.google.com/d/msgid/xarray/9aed84b5-5627-f4d6-1d82-c74bc4830deb%40gmail.com.

[Bane Sullivan]

Hi there! Ryan forwarded this to me as someone familiar with ParaView (and being from Kitware). 

From a quick glance, it seems the NetCDF file has the coordinates listed as X and Y. In order for ParaView to load it in on the XZ plane, the Y coordinates array likely needs to be switched to Z. I assume the `axis` property is just a label and not used for interpretation. Take this with a grain of salt though as I am not too familiar with the specifics of the NetCDF format.

To have ParaView show the data on the correct XZ plane (as is, without altering the data), you can do a representation rotation (see attached screenshot). To do this, search for “transformation” on the properties of the loaded dataset and apply a 90-degree rotation about the X-axis.

Hope this helps!

Bane Sullivan

[Marco Miani]

Dear Bane, 

Thank you so much for your kind answer. 

The rotation is something I had in mind too. Until I won't find  a more "consistent" solution, I will go for that and see. 

I have downloaded an empty netcdf template file from NOAA repository, a series of vertical profiles with the correct axis definition: file is attached. I will try to mimic that structure and see how far I can get.

Again, thanks for your precious help - appreciated! 
Have a great day ahead. 
Marco  

To view this discussion on the web visit https://groups.google.com/d/msgid/xarray/4bbadf70-4a8a-406b-99aa-2c539cd120afn%40googlegroups.com.

[Bane Sullivan]

Hi Marco, after looking at this a bit more and talking with some others more familiar with the NetCDF format and VTK's (and ParaView's by extension) reader for the format, it appears there could be a bug with the file reader, and your file is actually labeled correctly using the CF convention. I'll continue looking into this and potentially open an issue in VTK if there is a problem. I will check back.

In the meantime, I'd suggest only applying the rotation during rendering and leave the file as-is.

- Bane

[Marco Miani]

Hi Bane. 

Thanks for your feedback and for the effort you are putting into this. 

Ok, I am applying a mere rotation so far and will do until the bug is resolved. 

Have a good weekend. 

Ciao, 
M.

To view this discussion on the web visit https://groups.google.com/d/msgid/xarray/d88ec5de-6ac4-4daf-98a8-8107db3cee03n%40googlegroups.com.

[Bane Sullivan]

For the records, I posted an issue in VTK about this: https://gitlab.kitware.com/vtk/vtk/-/issues/18035

I'm not sure when it may be addressed.