iris.analysis.stats.pearsonr metadata error

[EnergyResearcher]

Hello,
I am trying to calculate correlation between two cubes. I have regridded them so that they are on the same grid but I get an error which seems to be about metadata. I can't understand from the error what I should be fixing in metadata to get the function to work.
The code I use is cc = iris.analysis.stats.pearsonr(cube1, cube2) which throws the following error:

ValueError: Insufficient matching coordinate metadata to resolve cubes, cannot map dimension (0,) of the RHS cube 'unknown' to the LHS cube 'unknown'.

I tried to search through the source code, the error seems to be in the module resolve.py but the code there is a bit hard to navigate for me. Could anyone point me to how I could fix the error?

EDIT:
cube1 print:

significant_height_of_wind_and_swell_waves / (unknown) (time: 6552; latitude: 361; longitude: 720)
    Dimension coordinates:
        time                                                x               -               -
        latitude                                            -               x               -
        longitude                                           -               -               x
    Scalar coordinates:
        forecast_period                                0 hours
        originating_centre                             European Centre for Medium Range Weather Forecasts

cube2 print:

eastward_wind / (m s-1)             (time: 29224; latitude: 361; longitude: 720)
    Dimension coordinates:
        time                             x                -               -
        latitude                         -                x               -
        longitude                        -                -               x
    Scalar coordinates:
        height                      10.0 m
    Cell methods:
        mean                        area
        point                       time
    Attributes:
        Conventions                 CF-1.7 CMIP-6.2
        activity_id                 CMIP
        branch_method               standard
        branch_time_in_child        0.0
        branch_time_in_parent       0.0
        cmor_version                3.5.0
        comment                     Eastward component of the near-surface (usually, 10 meters)  wind
        contact                     cmip6-mpi-esm@dkrz.de
        data_specs_version          01.00.30
        experiment                  all-forcing simulation of the recent past
        experiment_id               historical
        external_variables          areacella
        forcing_index               1
        frequency                   3hrPt
        further_info_url            https://furtherinfo.es-doc.org/CMIP6.MPI-M.MPI-ESM1-2-HR.historical.no...
        grid                        gn
        grid_label                  gn
        initialization_index        1
        institution                 Max Planck Institute for Meteorology, Hamburg 20146, Germany
        institution_id              MPI-M
        license                     CMIP6 model data produced by MPI-M is licensed under a Creative Commons...
        mip_era                     CMIP6
        nominal_resolution          100 km
        parent_activity_id          CMIP
        parent_experiment_id        piControl
        parent_mip_era              CMIP6
        parent_source_id            MPI-ESM1-2-HR
        parent_time_units           days since 1850-1-1 00:00:00
        parent_variant_label        r1i1p1f1
        physics_index               1
        product                     model-output
        project_id                  CMIP6
        realization_index           1
        realm                       atmos
        references                  'MPI-ESM: Mauritsen, T. et al. (2019), Developments in the MPI‐M Earth...
        source                      'MPI-ESM1.2-HR (2017): \naerosol: none, prescribed MACv2-SP\natmos: ECHAM6.3...
        source_id                   MPI-ESM1-2-HR
        source_type                 AOGCM
        sub_experiment              none
        sub_experiment_id           none
        table_id                    3hr
        table_info                  Creation Date:(09 May 2019) MD5:e6ef8ececc8f338646ebfb3aeed36bfc
        title                       MPI-ESM1-2-HR output prepared for CMIP6
        variable_id                 uas
        variant_label               r1i1p1f1

cube1 coords print:

[DimCoord(array([78888, 78889, 78890, ..., 85437, 85438, 85439]), standard_name='time', units=Unit('hours since 1970-01-01 00:00:00', calendar='gregorian')),
DimCoord(array([ 90. ,  89.5,  ......., -88.5, -89. , -89.5, -90. ]), bounds=array([[ 90.  ,  89.75],
       [ 89.75,  89.25],
       ....
       [-89.25, -89.75],
       [-89.75, -90.  ]]), standard_name='latitude', units=Unit('degrees'), coord_system=GeogCS(6371229.0)), 
DimCoord(array([  0. ,   0.5,358.5, 359. , 359.5]), bounds=array([[-2.5000e-01,  2.5000e-01],
       [ 2.5000e-01,  7.5000e-01],
       [ 7.5000e-01,  1.2500e+00],
       ...,
       [ 3.5825e+02,  3.5875e+02],
       [ 3.5875e+02,  3.5925e+02],
       [ 3.5925e+02,  3.5975e+02]]), standard_name='longitude', units=Unit('degrees'), coord_system=GeogCS(6371229.0), circular=True), 
DimCoord(array([0], dtype=int32), standard_name='forecast_period', units=Unit('hours')), 
AuxCoord(array(['European Centre for Medium Range Weather Forecasts'], dtype='<U50'), standard_name=None, units=Unit('no_unit'), long_name='originating_centre')]

cube2 coords print:

[DimCoord(array([45655.125, 45655.25 , 45655.375, ..., 49307.75 , 49307.875,
       49308.   ]), standard_name='time', units=Unit('days since 1850-1-1 00:00:00', calendar='proleptic_gregorian'), long_name='time', var_name='time'),
DimCoord(array([ 90. ,  89.5,  89. ,  88.5,  88. ,  87.5,  87. ,  86.5,  86. ,
        85.5,  85. ,  84.5,  84. ,  83.5,  83. ,  82.5,  82. ,  81.5, .....  -90. ]), bounds=array([[ 90.  ,  89.75],
       [ 89.75,  89.25],
       [ 89.25,  88.75], 
....
[-89.25, -89.75],
       [-89.75, -90.  ]]), standard_name='latitude', units=Unit('degrees'), coord_system=GeogCS(6371229.0)),
DimCoord(array([  0. ,   0.5,   1. ,   1.5,   2. ,   2.5,   3. ,   3.5,   4. , .... 359. , 359.5]), bounds=array([[-2.5000e-01,  2.5000e-01],
       [ 2.5000e-01,  7.5000e-01],
       [ 7.5000e-01,  1.2500e+00],
       ...,
       [ 3.5825e+02,  3.5875e+02],
       [ 3.5875e+02,  3.5925e+02],
       [ 3.5925e+02,  3.5975e+02]]), standard_name='longitude', units=Unit('degrees'), coord_system=GeogCS(6371229.0), circular=True), 
DimCoord(array([10.]), standard_name='height', units=Unit('m'), long_name='height', var_name='height', attributes={'positive': 'up'})]

[pp-mo]

Obviously, can't say a whole lot without more details.
Could we at least have a cube printout ?

However, it sounds like it can't identify which dimension in cube A matches that in cube B.
So, at a guess, it could because there are no coordinates at all attached to some dimensions. Or, possibly, some dims coordinates have identical identities, i.e. same names and no distinguishing attributes.
( in which case, it might not be clear from a cube printout : might have to print coords as well )

[EnergyResearcher]

Sorry, I've updated the post with the printouts. (The last two odd dimensions in cube1 appeared when I read the grib files with iris. Reading netCDF files didn't result in these dimensions but coordinates were not assigned x and y axes so I couldn't regrid it.)

[rcomer]

I assume you do some subsetting so that the time dimensions are the same length. The problem resolve is flagging is that your time coordinates have different units. You could fix this with something like

import cf_units

# Make sure both cubes have the same calendar, otherwise convert_units doesn't work.
new_unit = cf_units.Unit(cube2.coord("time").units.origin, calendar="gregorian")
cube2.coord("time").units = new_unit

cube2.convert_units(cube1.coord("time"))

However, you will still get an error if the time coordinate points and bounds don't match, so it may be easier to replace the whole coordinate:

cube2.remove_coord("time")
cube2.add_dim_coord(cube1.coord("time"), 0)

[EnergyResearcher]

Thank you! I the latter worked but with the warning: Collapsing a non-contiguous coordinate. Metadata may not be fully descriptive for 'time'.. I was expecting to get a single value, I converted the scalar cube to xarray DataArray but I still don't see the value. When I try to show the numpy arrays behind the xr.DataArray the kernel crashes. Here is the code:

corr = stats.pearsonr(cube1, cube2)
corr_xr = xr.DataArray.from_iris(corr)
corr_xr.values

[rcomer]

I don't know anything about xarray but you can get your value out with corr.data.

The warning about non-contiguous coordinates always happens if your coordinates don't have bounds, and I always happily ignore it.

[EnergyResearcher]

I tried corr.data before converting to xarray, corr.data also crashes the kernel

[pp-mo]

corr.data also crashes the kernel

OK that is concerning. I guess you mean it aborts the whole Python session ?

It might be because the whole result (or calculating it) is simply too big for memory
I think the calculation is fully lazy.
If so, maybe it is possible to calculate only parts of it ?
However, from the dimensions in your printouts, it doesn't seem that big. As @rcomer says, it looks as if you must be getting two equal-length time series first, somehow.

Another possibility is that there might be an unexpected rounding error leading to square root of a negative result - but I don't think that would actually crash the calculation, let alone kill the session.

If it is a memory issue, it might also be due to problems with Dask, if the chunking in the two source cubes does not align well.

It's hard to dig for possible problems without the original data, so I'm really just guessing again ...
Is there any chance you can give us this whole actual problem in a reproducible form ?

[EnergyResearcher]

Thank you for the reply. Yes, it aborts the whole python session.

I think while getting the example ready, I managed to isolate the problem. Before calculating the correlation, I was regridding one of the cubes. corr.data crashes only when cube2 was regridded to cube1 grid but doesn't crash if cube1 is regridded to cube2 grid. I still haven't figured out why saving regridded cube crashes (doesn't matter which cube is regridded).

The cubes are here.
Here is the reproducible code:

import iris
from iris.analysis import stats
import xarray as xr
import numpy as np
import time
import numpy.ma as ma
import iris.plot as iplt

# download the cubes from: https://drive.google.com/drive/folders/1VTUx_SfcPku3EWIdxt2i2mFvidtIK-5P?usp=sharing
cube1 = iris.load('cube1.nc')[0]
cube2 = iris.load('cube2-1.nc')[0]

# check the coordinates of both cubes for regridding
cubes = [cube1, cube2]
for cube in cubes:
    info = {'bounds': [], 'cs': []}
    for coord in cube.coords():
        if cube.coord(coord).has_bounds():
            info['bounds'].append(coord.name())
        if coord.coord_system is not None:
            info['cs'].append(coord.coord_system)
        
    print('Coordinates of cube "{}" which have bounds: {}'.format(cube.name(0),info['bounds']))
    print('Coordinates of cube "{}" which have CS: {}'.format(cube.name(0),info['cs']))

def arrange_cs(cube, cs=True, bounds=True):
    new_cs = iris.coord_systems.GeogCS(iris.fileformats.pp.EARTH_RADIUS)
    if cs:
        cube.coord(axis='x').coord_system = new_cs
        cube.coord(axis='y').coord_system = new_cs
    if bounds:
        cube.coord('latitude').guess_bounds()
        cube.coord('longitude').guess_bounds()
    
    return cube

cube1 = arrange_cs(cube1)

# corr.data doesn't crash if cube1 is regridded with cube2 as a target grid
regridded = cube1.regrid(cube2,iris.analysis.AreaWeighted())

# corr.data DOES crash if cube 2 is regridded with cube1 as a target grid
#regridded = cube2.regrid(cube1,iris.analysis.AreaWeighted())

# this crashes the kernel with either regridding order
#iris.save(regridded, 'cube2-2.nc') 

def rearrange_time(cube1, cube2):
    
    # limit cube2 to the timeframe of cube1
    cube1_xr = xr.DataArray.from_iris(cube1)
    cube2_xr = xr.DataArray.from_iris(cube2)
    cube2_xr = cube2_xr.sel(time=slice(cube1_xr.time.min(), cube1_xr.time.max()))
    cube1 = cube1_xr.to_iris()
    cube2 = cube2_xr.to_iris()
    
    # make units uniform
    if not cube1.coord('time').units == cube2.coord('time').units:
        cube1.remove_coord('time')
        cube1.add_dim_coord(cube2.coord('time'),0)
    
    return cube1, cube2

# if cube1 was regridded:
cube1, cube2 = rearrange_time(regridded, cube2)
#if cube1 was regridded:
#cube1, cube2 = rearrange_time(cube1,regridded)

corr = stats.pearsonr(cube1,cube2)

corr.data

When I calculated the correlation I realised that the uas cube has values over land (and wave height has NaNs there) which were included in correlation calculation. I tried to mask the land values which worked when looking at the plots. However, when I calculate the correlation with the mask, the result is the same as without the mask. The relevant code is below:

cube1_xr = xr.DataArray.from_iris(cube1)
np_mask = cube1_xr.isnull().to_numpy()
cube1_masked = cube1.copy()
cube2_masked = cube2.copy()
cube1_masked.data = ma.array(cube1.data, mask=np_mask)
cube2_masked.data = ma.array(cube2.data,mask=np_mask)

iplt.pcolormesh(cube1_masked[0]) # significant wave height

iplt.pcolormesh(cube2_masked[0]) # uas

corr_masked = stats.pearsonr(cube1_masked,cube2_masked, common_mask=True)
corr = stats.pearsonr(cube1,cube2)
corr == corr_masked # returns True

[pp-mo]

Archiving "answered" Q+As