MCA errors

[Murk89]

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)

[nicrie]

The error is due to an unfortunate limitation of this package that the spatial coordinates have to be named lon and lat.

Currently, your main spatial coordinates are called west_east and south_north. You should try renaming them to lon and lat. Since lon and lat already exist in your Dataset you may have to rename these as well (or even drop them, not sure).
For example, before doing the analysis:

var = var.rename({'lat': 'lat_old', 'lon': 'lon_old'})
var = var.rename({'south_north': 'lat', 'west_east': 'lon'})

If you're still running into problems try to boil down your DataArrays to a form where they only have three dimensions called time, lat and lon.

---

Note: In this package, lon and lat are unfortunately hard coded. In xeofs (which will also feature MCA soon), you will have full flexibility to choose any given coordinate combination within your DataArrays. So perhaps watch out for coming releases :)

[Murk89]

Hi,
Many thanks for getting back. So after renaming the coords and dims as suggested, the code runs fine. However I am afraid that there is some issue with the output. May be due to renaming the south_north and west_east as lat and lon, the plotting method doesn't correctly pick up WRF's projection?
The first figure is the output I get from mca.plot, though my WRF domain is as shown in the second figure.

[nicrie]

I think you're right. The plot function is a convenient way of quickly inspecting the modes, but it seems that for your use case it is not general enough and indeed struggles to pick up the spatial coordinates/projection correctly.

In any case, you can easily plot the modes yourselves, e.g.:

from cartopy.crs import NearsidePerspective, PlateCarree

# ... your analysis

pcs = mca.pcs()                           
eofs = mca.eofs()  

proj = NearsidePerspective(central_longitude=15, central_latitude=65)
mode = 1

fig = plt.figure(figsize=(14, 7))
ax_eof_left = fig.add_subplot(221, projection=proj)
ax_eof_right = fig.add_subplot(222, projection=proj)
ax_pc_left = fig.add_subplot(223)
ax_pc_right = fig.add_subplot(224)

eofs['left'].sel(mode=mode).plot(ax=ax_eof_left, transform=PlateCarree())
eofs['right'].sel(mode=mode).plot(ax=ax_eof_right, transform=PlateCarree())
pcs['left'].sel(mode=mode).plot(ax=ax_pc_left)
pcs['right'].sel(mode=mode).plot(ax=ax_pc_right)

Let me know if it works.

[Murk89]

Hi Niclas,

Thanks again for helping out.
So I tried to get the projection right by using wrf-python along with xmca. Is a bit lengthy (being utterly honest), but I am new to coding so please forgive my clumsiness!

Here is the final script;

import xarray as xr
import matplotlib.pyplot as plt
from xmca.xarray import xMCA  # use with xr.DataArray
import netCDF4
from netCDF4 import Dataset
import wrf                   ##wrf-python package 
from wrf import (getvar, get_cartopy, cartopy_xlim,cartopy_ylim, latlon_coords)

var=xr.open_dataset("F:\\era5_2000_2020_vars_salem.nc") ##this is a concatenated WRF dataset. Concatenation was done using the salem package.
t2=var.T2C
snow=var.SNOWH

t2_new = t2.rename({'lat': 'lat_old', 'lon': 'lon_old'}) ##for resolving the first error mentioned in this issue
t2_new=t2_new.rename({'south_north': 'lat', 'west_east': 'lon'}) 

snow_new=snow.rename({'lat': 'lat_old', 'lon': 'lon_old'})
snow_new=snow_new.rename({'south_north': 'lat', 'west_east': 'lon'})

##MCA two fields
mca = xMCA(t2_new, snow_new)  # MCA of field A and B
mca.solve(complexify=False)   # True for complex MCA

eigenvalues = mca.singular_values()     # singular vales
pcs = mca.pcs()                         # expansion coefficient (PCs)
eofs = mca.eofs()   

##getting the WRF projection right for plotting

data=Dataset("F:\\rcp4.5\\D02-2090-2100\\wrfout_d02_2090-11-01_00_00_00.nc.nc") 
##doesn't get the lat,lon from the dataset concatenatd using salem hence have to read in an original wrfout file

sst = getvar(data, "SST") ##select a random variable from the WRF dataset 'data'
lats, lons = latlon_coords(sst)
cart_proj = get_cartopy(sst) ##gets the correct WRF projection
proj=cart_proj
 
##plotting the PCs
mode = 2

fig = plt.figure(figsize=(14, 7))
ax_eof_left = fig.add_subplot(221, projection=proj)
ax_eof_right = fig.add_subplot(222, projection=proj)
ax_pc_left = fig.add_subplot(223)
ax_pc_right = fig.add_subplot(224)

eofs['left'].sel(mode=mode).plot(ax=ax_eof_left, )
eofs['right'].sel(mode=mode).plot(ax=ax_eof_right, )
pcs['left'].sel(mode=mode).plot(ax=ax_pc_left)
pcs['right'].sel(mode=mode).plot(ax=ax_pc_right)    

[nicrie]

Looks valid to me, no? Perhaps you want to add coastlines to identify your region of interest:

# analysis

fig = plt.figure(figsize=(14, 7))
# define axes ...

# add coastlines
ax_eof_left.coastlines()
ax_eof_right.coastlines()

# rest of plotting
# ...
eofs['left'].sel(mode=mode).plot(ax=ax_eof_left, )

Do you think your issue with plotting is solved by this?

[Murk89]

Hi Niclas,

I think the plotting is correct for now now. Though I might work on the plotting further with the salem package, because it produces better maps with WRF grid (as shown in the first set of figures I had provided). Will update here if I can get salem to produce them.
Thanks.

[gkb999]

Hi team, the package xMCA is simply amazing I should say. It keeps things really straightforward.
However, there is limited information in API references regarding projection. Although I referred to previous comments and modified my plot, still I need to tune it.
Any help/support would be greatly appreciated

What I got from :
from cartopy.crs import Orthographic # for different map projections

map projections for "left" and "right" field

projections = {
'left': Orthographic(-180, -90),
'right': Orthographic(-180, -90),
}
pkwargs = {
"figsize" : (8, 3),
"orientation" : 'horizontal',
'cmap_eof' : 'RdBu_r', # colormap amplitude
"projection" : projections, }
is:

Also, any projection I am changing, be it, orthographic, stereo, or SouthPolarStereo - the graph doesn't change. What could be the reason?

What I want - a bit clean this way:

[nicrie]

Thank you for your kind words! The convenience function plot unfortunately sometimes shows some undesirable/irregular behaviour, especially when using more fancy projections than PlateCarree. I haven't had (and probably won't have) the time to fix this issue anytime soon, but you can refer to the snippet above to extract the EOFs/PCs and create your own plot using matplotlib and cartopy (or really any other library you like :) ).

from cartopy.crs import NearsidePerspective, PlateCarree

# ... your analysis

pcs = mca.pcs()                           
eofs = mca.eofs()  

proj = NearsidePerspective(central_longitude=15, central_latitude=65)
mode = 1

fig = plt.figure(figsize=(14, 7))
ax_eof_left = fig.add_subplot(221, projection=proj)
ax_eof_right = fig.add_subplot(222, projection=proj)
ax_pc_left = fig.add_subplot(223)
ax_pc_right = fig.add_subplot(224)

eofs['left'].sel(mode=mode).plot(ax=ax_eof_left, transform=PlateCarree())
eofs['right'].sel(mode=mode).plot(ax=ax_eof_right, transform=PlateCarree())
pcs['left'].sel(mode=mode).plot(ax=ax_pc_left)
pcs['right'].sel(mode=mode).plot(ax=ax_pc_right)

[gkb999]

Wonderful @nicrie !!! This is what I was exactly looking for.

Thanks alot.

[gkb999]

Thank you for your kind words! The convenience function plot unfortunately sometimes shows some undesirable/irregular behaviour, especially when using more fancy projections than PlateCarree. I haven't had (and probably won't have) the time to fix this issue anytime soon, but you can refer to the snippet above to extract the EOFs/PCs and create your own plot using matplotlib and cartopy (or really any other library you like :) ).

from cartopy.crs import NearsidePerspective, PlateCarree

# ... your analysis

pcs = mca.pcs()                           
eofs = mca.eofs()  

proj = NearsidePerspective(central_longitude=15, central_latitude=65)
mode = 1

fig = plt.figure(figsize=(14, 7))
ax_eof_left = fig.add_subplot(221, projection=proj)
ax_eof_right = fig.add_subplot(222, projection=proj)
ax_pc_left = fig.add_subplot(223)
ax_pc_right = fig.add_subplot(224)

eofs['left'].sel(mode=mode).plot(ax=ax_eof_left, transform=PlateCarree())
eofs['right'].sel(mode=mode).plot(ax=ax_eof_right, transform=PlateCarree())
pcs['left'].sel(mode=mode).plot(ax=ax_pc_left)
pcs['right'].sel(mode=mode).plot(ax=ax_pc_right)

@nicrie Just to get it to your notice,
the same code doesn't work well for complex plotting (I mean to say, when set complex to True).

The following is the error:
Plotting requires coordinates to be numeric, boolean, or dates of type numpy.datetime64, datetime.datetime, cftime.datetime or pandas.Interval. Received data of type complex128 instead.

A simple plot function does work, but again, it doesn't look quite clean.
Any idea to sort this out?

[nicrie]

@GIRIJA-KALYANI this is probably not a bug but - as the error message implies - xarray does not know how to plot complex data. Since you set complexify=True, the output of mca.eofs() and mca.pcs() are the complex EOFs and PCs. So you either have to use the EOFs/PCs to real or imaginary part only (e.g. eofs['left'].real) or you plot the amplitude and phase (which is by definition real). To get the amplitude use mca.spatial_amplitudes() and mca.spatial_phase().

[gkb999]

Thanks a ton for the suggestion @nicrie . I'm trying to plot complex EOFs, and only eofs['left'].real works. Using pca.spatial_amplitudes() or pca.spatial_phase() - doesn't seem to work at least me. Maybe I'm going wrong somewhere.
The error I get: 'xMCA' object has no attribute 'spatial_amplitudes'

When I try pca.plot(mode=1) - it works

What could be the error?

Also, when I try using varimax rotation
I get the following error.
.

Also, when I look into the data of complex EOFs, I see nan+nanj values,

any insights on how the code handles nan values please?

Thanks in advance

[nicrie]

I'm sorry you encounter these errors, let me go through them one by one:

1. When performing complex EOF, (in contrast to MCA) you can retrieve only one field (called left). Therefore, eofs['left'].real and eofs['left'].imag should work while eofs['right'].real and eofs['right'].imag will raise an error.
2. You can retrieve the spatial amplitude and phase through mca.spatial_amplitude() and mca.spatial_phase() - in my previous answer there was a typing, it should be amplitude (singular) in contrast to amplitudes. For reference, you can check the API for all possible methods. However, I'm surprised that mca.spatial_phase() fails - can you copy the error message please?
3. The correct usage for rotation is just pca.rotate(), no assignment needed. Actually, the method does not return anything, the rotation is done internally. When assigning pca = pca.rotate(), the variable pca will just be None. Therefore, any further analysis will fail because you have basically overwritten your pca class. So just write the following:

pca.rotate(n_rot=10, power=1)
expvar_rot = pca.explained_variance()
....

4. All the NaNs you see when inspecting the eofs are probably right and just represent "invalid" grid points such as the part of Antarctica that has no measurements. If there were only NaNs in the result pca.plot(1) would only show a blank figure without data. Since you actually get something quite colorful (and beautiful :)), I wouldn't worry much about the shown NaNs.

Does that answer your questions?

[gkb999]

I'm sorry you encounter these errors, let me go through them one by one:

1. When performing complex EOF, (in contrast to MCA) you can retrieve only one field (called left). Therefore, eofs['left'].real and eofs['left'].imag should work while eofs['right'].real and eofs['right'].imag will raise an error.
2. You can retrieve the spatial amplitude and phase through mca.spatial_amplitude() and mca.spatial_phase() - in my previous answer there was a typing, it should be amplitude (singular) in contrast to amplitudes. For reference, you can check the API for all possible methods. However, I'm surprised that mca.spatial_phase() fails - can you copy the error message please?
3. The correct usage for rotation is just pca.rotate(), no assignment needed. Actually, the method does not return anything, the rotation is done internally. When assigning pca = pca.rotate(), the variable pca will just be None. Therefore, any further analysis will fail because you have basically overwritten your pca class. So just write the following:

pca.rotate(n_rot=10, power=1)
expvar_rot = pca.explained_variance()
....

4. All the NaNs you see when inspecting the eofs are probably right and just represent "invalid" grid points such as the part of Antarctica that has no measurements. If there were only NaNs in the result pca.plot(1) would only show a blank figure without data. Since you actually get something quite colorful (and beautiful :)), I wouldn't worry much about the shown NaNs.

Does that answer your questions?

I'm sorry for responding late.
But yes, most of my questions are answered. Although I'm aware of calling only the 'left' field for EOF analysis (As it deals with a single field), I wasn't aware of eofs['left'].imag for imaginary mode when I want to see complex imaginary. Thanks for that.
2) mca.spatial_amplitude()andmca.spatial_phase() work well now, but a bit working around to look at the plots.
3) So, I don't have to assign anything . Like --> pca.rotate() ?? (did I get it well? No need to use pca.rotate(n_rot=10, power=1)??)
4) True, I don't have to worry much about the NaNs anyway :)

Plus, sometimes I get SVD did not converge as runtime error, this mostly happens when I set complex to True initially, then change again to False , just to notice the differences

I also want to see correlation maps (I have seen this from API that the module gives this option, how do I derive one and visualize??
Also, I'm interested in plotting heterogeneous_patterns and homogenous_patterns - any light on that?

Thanks a ton in advance

[nicrie]

If you want to rotate, just call

pca.rotate(<your arguments here>)

The following will not work

pca = pca.rotate(<your arguments here>)

In any case, the arguments are still needed!

The plotting of homogeneous/heterogeneous and or correlation maps works exactly as for the spatial amplitudes. The objects are xr.DataArray. So if you manage to plot the spatial amplitudes/phases just replace these with the quantity you are interested in. Check the API which quantities you can retrieve.

---

As a side note: please consider switching over to xeofs as it is starting to incorporate xmca. For now, you can do e.g. EOF, MCA and rotation with xeofs. You'll also find some examples how to do it, so getting started should be easy.