NETCDF multi-arrays and integration with xarray

jlmaurer commented 2 years ago

I wasn't sure how doing a single file with multiple variables would work with your netcdf writer, but I noticed that pygmt has native support for xarray dataarrays, so I integrated xarray into the reading/writing of the multi-array for strain rate. At some point we should also include the interpolated velocities in the file, but saving that for a later PR. I think I got everything updated but if you find a bug let me know.

Description

Refactored code to use xarray Datasets for storing and writing out the strain rate variables
Updated plotting to use the dataarrays (surprisingly little to change as pygmt natively handles xarray DataArrays)
Also updated the max_shear, dilatation, and I2 calculation to be a bit simpler as numpy array operations are natively vectorized and I was getting some wild values

Motivation and Context

This allows for a single file to be written containing all the strain rate components and derived quantities on a single grid.
This also updates the derived quantities calculation

How Has This Been Tested?

Ran examples with gpsgridder and huang, and then ran the comparison script

jlmaurer commented 2 years ago

@kmaterna Not sure what xarray uses for registration, but since pygmt handles it natively it may not be an issue. I've used xarray a lot recently and it seems fairly robust, strong user base, etc. Also it talks back and forth with numpy natively.

kmaterna commented 2 years ago

This works on my computer too! I do like the multi-dimension netcdf file.

Tiny comments:

can I trouble you to change the doc-string on velocity_io.read_multiple_strain_files? You know the arguments better than me.
the dilatation and max shear values for a Delaunay example are actually different from before. Dilatation looks smaller and max shear looks larger. Does it make sense that they changed? It's possible there was an instability inherent in computing eigenvectors all along, eek!

After merging this PR, I'll do the following:

future: adjust the function "test_strain_functions.test_solid_body_rotation()" because one function was replaced
future: change the I2 max print statement to be "np.nanmax"
future: add long form pygmt arguments
future: document a line for extracting a layer into valid GRD file, respecting pixel node registration

jlmaurer commented 2 years ago

@kmaterna sounds good. I will update the doc-string. Re the values of dilatation and max shear, I did a simple test comparing the eigenvalue and direct method and noticed that the difference is exactly one-half (0.5), which I also notice you multiple the eigenvalue difference here: https://github.com/kmaterna/Strain_2D/blob/be7f5b2d1276834b150c272c4acaa3fbce82f03e/Strain_Tools/strain/strain_tensor_toolbox.py#L74

I double-checked the math for the direct method online and don't see a factor of two; however, I know there is this ambiguity between engineering shear strain and tensorial shear strain, so I'm not sure if we need the factor of 1/2 or not!

It's the exact same issue for dilatation except opposite. My version uses 1/2 the trace while yours is simply using the trace.

Not sure if there is a paper or textbook we could reference for standard conventions?

jlmaurer commented 2 years ago

Oh BTW, Github hint: if you 'ping' me using "@" I'll get a notification about it. :)

kmaterna / Strain_2D