nicrie / xmca Goto Github PK

View Code? Open in Web Editor NEW

50.0 2.0 15.0 7.06 MB

Maximum Covariance Analysis in Python

License: MIT License

Python 100.00%

pca mca xarray python numpy principal-component-analysis correlation rotation maximum complex

xmca's People

Contributors

Stargazers

Watchers

Forkers

xmca's Issues

improve docs with out-of-box MWE

basically use per-processed era5 data from the test folder

Does xMCA work on gridded data too?

Hi,
My dataset looks like:

xarray.Dataset
Dimensions:
lon: 171lat: 128time: 3285bnds: 2
Coordinates:
lon (lon) float64 -5.344e+05 -5.281e+05 ... 5.281e+05
lat (lat) float64 -1.959e+06 ... -1.166e+06
time (time) datetime64[ns] 2013-01-01 ... 2021-12-31

Lat and Lon in 'meters' (below)

When I applied EOFs, on this, it however, did not fail.
But, the result looked strange.
Please check the eof[0] pattern...

the eof[0] lokks like:

xarray.DataArray'EOFs'lon: 171mode: 6
array([[-1.76054379e-22, -3.28724969e-19, -1.11229007e-16,
        -1.27920154e-17, -2.97548278e-16,  1.83119557e-16],
       [-6.15561168e-21,  3.36886857e-18, -8.68632332e-17,
        -1.61311084e-16, -2.85747382e-17, -1.29165757e-16],
       [-1.66259177e-20,  1.44649751e-17, -4.54359525e-16,
        -2.18833121e-16, -3.33146482e-16,  2.95736462e-16],
       ...,
       [ 0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
        -0.00000000e+00, -0.00000000e+00,  0.00000000e+00],
       [ 0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
        -0.00000000e+00, -0.00000000e+00,  0.00000000e+00],
       [ 0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
        -0.00000000e+00, -0.00000000e+00,  0.00000000e+00]])
Coordinates:
lon
(lon)
float64
-5.344e+05 -5.281e+05 ... 5.281e+05
array([-534375., -528125., -521875., -515625., -509375., -503125., -496875.,
       -490625., -484375., -478125., -471875., -465625., -459375., -453125.,
       -446875., -440625., -434375., -428125., -421875., -415625., -409375.,
       -403125., -396875., -390625., -384375., -378125., -371875., -365625.,
       -359375., -353125., -346875., -340625., -334375., -328125., -321875.,
       -315625., -309375., -303125., -296875., -290625., -284375., -278125.,
       -271875., -265625., -259375., -253125., -246875., -240625., -234375.,
       -228125., -221875., -215625., -209375., -203125., -196875., -190625.,
       -184375., -178125., -171875., -165625., -159375., -153125., -146875.,
       -140625., -134375., -128125., -121875., -115625., -109375., -103125.,
        -96875.,  -90625.,  -84375.,  -78125.,  -71875.,  -65625.,  -59375.,
        -53125.,  -46875.,  -40625.,  -34375.,  -28125.,  -21875.,  -15625.,
         -9375.,   -3125.,    3125.,    9375.,   15625.,   21875.,   28125.,
         34375.,   40625.,   46875.,   53125.,   59375.,   65625.,   71875.,
         78125.,   84375.,   90625.,   96875.,  103125.,  109375.,  115625.,
        121875.,  128125.,  134375.,  140625.,  146875.,  153125.,  159375.,
        165625.,  171875.,  178125.,  184375.,  190625.,  196875.,  203125.,
        209375.,  215625.,  221875.,  228125.,  234375.,  240625.,  246875.,
        253125.,  259375.,  265625.,  271875.,  278125.,  284375.,  290625.,
        296875.,  303125.,  309375.,  315625.,  321875.,  328125.,  334375.,
        340625.,  346875.,  353125.,  359375.,  365625.,  371875.,  378125.,
        384375.,  390625.,  396875.,  403125.,  409375.,  415625.,  421875.,
        428125.,  434375.,  440625.,  446875.,  453125.,  459375.,  465625.,
        471875.,  478125.,  484375.,  490625.,  496875.,  503125.,  509375.,
        515625.,  521875.,  528125.])
lat () float64 -1.959e+06
standard_name : projection_y_coordinate
long_name : y coordinate of projection
units : metre
axis : Y
array( -1959375.)
mode (mode) int32 1 2 3 4 5 6
array([1, 2, 3, 4, 5, 6])
Indexes:
lon PandasIndex
PandasIndex(Float64Index([-534375.0, -528125.0, -521875.0, -515625.0, -509375.0, -503125.0,
              -496875.0, -490625.0, -484375.0, -478125.0,              ...
               471875.0,  478125.0,  484375.0,  490625.0,  496875.0,  503125.0,
               509375.0,  515625.0,  521875.0,  528125.0],
             dtype='float64', name='lon', length=171))
mode PandasIndex
PandasIndex(Int64Index([1, 2, 3, 4, 5, 6], dtype='int64', name='mode'))

Here, the lat() is empty.
Why is it?
Should I reproject the data, as xMCA may not support this format of data?

Any assistance on this would be helpful.

Thanks in advance.

cartopy dependency is too restrictive

The very restrictive cartopy dependency makes it tricky to install into an existing conda environment

cartopy==0.18.0

I can see it was changed at coecms@896e0b5

When I changed it back to >=18.0 it installed fine using pip install --user and imported fine with cartopy-0.19.0.post1.

I ran the tests like so

python -m unittest discover -v -s tests/

but three of the tests didn't pass

test_save_load_cplx (integration.test_integration_xarray.TestIntegration) ... ERROR    
test_save_load_rot (integration.test_integration_xarray.TestIntegration) ... ERROR                                       
test_save_load_std (integration.test_integration_xarray.TestIntegration) ... ERROR

Some other error messages:

Error: Rotation process did not converge!

======================================================================                                                   
ERROR: test_save_load_cplx (integration.test_integration_xarray.TestIntegration)                                         
----------------------------------------------------------------------                                                   
Traceback (most recent call last):
  File "/home/502/aph502/.local/lib/python3.9/site-packages/parameterized/parameterized.py", line 533, in standalone_func
    return func(*(a + p.args), **p.kwargs)
  File "/home/502/aph502/code/python/xmca/tests/integration/test_integration_xarray.py", line 148, in test_save_load     
    rmtree(join(getcwd(), 'tests/integration/temp/'))
  File "/g/data3/hh5/public/apps/miniconda3/envs/analysis3-22.01/lib/python3.9/shutil.py", line 731, in rmtree           
    onerror(os.rmdir, path, sys.exc_info())
  File "/g/data3/hh5/public/apps/miniconda3/envs/analysis3-22.01/lib/python3.9/shutil.py", line 729, in rmtree           
    os.rmdir(path)
OSError: [Errno 39] Directory not empty: '/home/502/aph502/code/python/xmca/tests/integration/temp/'                     

======================================================================                                                   
ERROR: test_save_load_rot (integration.test_integration_xarray.TestIntegration)                                          
----------------------------------------------------------------------                                                   
Traceback (most recent call last):
  File "/home/502/aph502/.local/lib/python3.9/site-packages/parameterized/parameterized.py", line 533, in standalone_func
    return func(*(a + p.args), **p.kwargs)
  File "/home/502/aph502/code/python/xmca/tests/integration/test_integration_xarray.py", line 148, in test_save_load     
    rmtree(join(getcwd(), 'tests/integration/temp/'))
  File "/g/data3/hh5/public/apps/miniconda3/envs/analysis3-22.01/lib/python3.9/shutil.py", line 731, in rmtree           
    onerror(os.rmdir, path, sys.exc_info())
  File "/g/data3/hh5/public/apps/miniconda3/envs/analysis3-22.01/lib/python3.9/shutil.py", line 729, in rmtree           
    os.rmdir(path)
OSError: [Errno 39] Directory not empty: '/home/502/aph502/code/python/xmca/tests/integration/temp/'                     

======================================================================                                                   
ERROR: test_save_load_std (integration.test_integration_xarray.TestIntegration)                                          
----------------------------------------------------------------------                                                   
Traceback (most recent call last):
  File "/home/502/aph502/.local/lib/python3.9/site-packages/parameterized/parameterized.py", line 533, in standalone_func
    return func(*(a + p.args), **p.kwargs)
  File "/home/502/aph502/code/python/xmca/tests/integration/test_integration_xarray.py", line 148, in test_save_load     
    rmtree(join(getcwd(), 'tests/integration/temp/'))
  File "/g/data3/hh5/public/apps/miniconda3/envs/analysis3-22.01/lib/python3.9/shutil.py", line 731, in rmtree           
    onerror(os.rmdir, path, sys.exc_info())
  File "/g/data3/hh5/public/apps/miniconda3/envs/analysis3-22.01/lib/python3.9/shutil.py", line 729, in rmtree           
    os.rmdir(path)
OSError: [Errno 39] Directory not empty: '/home/502/aph502/code/python/xmca/tests/integration/temp/'                     

----------------------------------------------------------------------

but it wasn't obvious that this was an issue with the version of cartopy.

method to get the PC spectrum

thanks to @chjones2 for bringing this to my attention

currently only accepts 'lon' and 'lat' as label arguments

increase flexibility by allowing different labels for 'longitude' and 'latitude'

save_analysis currently does not save cos lat weighting

Just stumbled over this:

When saving a model via xmca.xarray.save_analysis and cosine latitude weighting was applied (apply_coslat), the current implementation does not invoke xmca.xarray.apply_coslat when the model gets loaded via xmca.xarray.load_analysis thus creating false PCs.

I hope to provide a fix to this soon.

May cause errors after using mca. normalize()

Hello, first of all, thank you for your contribution. However, I found a bug that may cause SVD to fail to calculate (due to the existence of NaN value). The details are as follows:

when MCA class is initialized, it will be called get_nan_cols and remove_ nan_ cols to remove the NaN value, but if you call mca.normalize() at this time. New NaN values can appear and be brought into SVD calculation. This is because if the value of each time step of a grid point in the input array is the same (For example, it never rains in a place), the value obtained after standardization is NaN, which causes SVD unable to solve problem.

error when loading and plotting PCA model

note to myself; add test suite for this case

SVD did not converge

Hi Niclas,

The XMCA worked fine when I used it directly on my raw data.
As it produced results of what I was expecting.
However, I tried to use processed data (like anomalies and detrend)
it gives the following error - SVG didn't converge

Does XMCA only accept raw data or is something wrong with my i/p?
This is how my data looks:

Even the previously worked data was in a similar format.
What could be the issue?

multivariate EOF analysis / MCA

add this feature in next release

note: this will be probably be a major change since it requires to rewrite the internal structure of the package and therefore will break backwards version compatibility

installation error

Hello Niclas,

I installed xmca but when I import the library I get this error:
libmkl_rt.so: cannot open shared object file: No such file or directory

any ideas or suggestions would be greatly appreciated.

Cheers,

Charles.

MCA errors

Hello,
I am trying to run the MCA with two variables, which are a climate model, WRF's output.

I get this error right after the bit:

mca.plot(mode=1, **pkwargs) :

ValueError: coordinate lon has dimensions ('south_north', 'west_east'), but these are not a subset of the DataArray dimensions ['lat', 'lon', 'mode']

Would really appreciate any help with this error. Many thanks.

# Load packages and data:
import xarray as xr
import matplotlib.pyplot as plt
import cartopy

var=xr.open_dataset("F:\\era5_2000_2020_vars_salem.nc")

t2=var.T2C
snow=var.SNOWH

#The variables, e.g., t2 is  structured as follows: 
t2:
<xarray.DataArray 'T2C' (time: 3512, south_north: 111, west_east: 114)>
Coordinates:
    lat          (south_north, west_east) float32 ...
    lon          (south_north, west_east) float32 ...
    xtime        (time) datetime64[ns] ...
  * time         (time) datetime64[ns] 2000-11-01 2000-11-02 ... 2020-04-29
  * west_east    (west_east) float64 -2.766e+05 -2.666e+05 ... 8.534e+05
  * south_north  (south_north) float64 -1.353e+05 -1.253e+05 ... 9.647e+05
Attributes:
    FieldType:    104
    MemoryOrder:  XY 
    description:  2m Temperature
    units:        C
    stagger:      
    pyproj_srs:   +proj=lcc +lat_0=64 +lon_0=10 +lat_1=64 +lat_2=68 +x_0=0 +y...
    coordinates:  XLONG XLAT XTIME

mca = xMCA(t2, snow)                  # MCA of field A and B
mca.solve(complexify=False)            # True for complex MCA


eigenvalues = mca.singular_values()     
pcs = mca.pcs()                           
eofs = mca.eofs()   

mca.set_field_names('t2','snow')
pkwargs = {'orientation' : 'vertical'}
mca.plot(mode=1, **pkwargs)

Feedback

I'm posting my feedback on behalf of Alex vM, who asked me to have a look at this package.
The feedback is subjective and you may disagree with some if not all suggestions below.

In general, the readme is well written. Clear and concise. I'd add a figure(s) though. For instance, after you call mca.plot() in your example.
Currently, this package is for developers or at least those who have python experience. Well, maybe covariance estimation requires knowledge of python and programming skills, but I recommend making sphinx documentation (and publishing it on readthedocs) to bring users. You've already written docs for each if not all functions and classes. The only step left is to wrap it in sphinx with default parameters and paths.
To give credits to your package and show that it's well maintained, I recommend adding badges (readthedocs, travis build and test coverage). Use CircleCI. Here is a config example (you need only build-pip job since it's the easiest).
tools folder must be in the xmca folder.
setup.py install requires must be read from requirements.txt file (example) and not hard-coded.
GPL-3 license is too restrictive. Consider BSD-3 or even MIT.
Each test that involves randomness must start with a numpy.random.seed. Currently, you're setting the seed globally. It's not a good idea because the test results depend on the test order, which, of course, should not happen.

Good luck!

Best,
Danylo

Why are there different outputs of nicrie / xmca compared with outputs of Yefee's xMCA?

Hi, nicrie.
Recently, I perform a svd analysis using nicrie / xmca code, the both shapes of left and right field matrixes were (time, lat, lon), and the value of time axis is (2000, 2001, ..., 2022).
I can get outputs of mca.eofs and mca.pcs, but I wonder that outputs were opposite compared with outputs of mca.expansionCoefs() of Yefee's xMCA.
Finally, I found that the same outputs can be obtained when time axis was reorded in inverse order, that is, (2022,2021,...,2000).
I want to know, which one is right? Is it a bug in code?
Thanks a lot.
the outputs of mca.pcs() using nicrie / xmca code were showed as

the outputs of mca.expansionCoefs() using Yefee's xMCA code were showed as