Python Matplotlib Difference between two NetCDF datasets

Question

I am trying to map the difference between climate simulation data and observed data over a set geographical area.

To create the map of just the climate simulation I am using this code

import matplotlib.pyplot as plt
import iris
import iris.plot as iplt
import cartopy
from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter
import iris.analysis.cartography

def main():   
        #bring in all the models we need and give them a name
        CCCma = '/exports/csce/datastore/geos/users/s0XXXX/Climate_Modelling/AFR_44_tas/ERAINT/1979-2012/tas_AFR-44_ECMWF-ERAINT_evaluation_r1i1p1_CCCma-CanRCM4_r2_mon_198901-200912.nc'

        #Load exactly one cube from given file    
        CCCma = iris.load_cube(CCCma)    

        #we are only interested in the latitude and longitude relevant to Malawi   
        Malawi = iris.Constraint(grid_longitude=lambda v: 31 <= v <= 36.5, \
                             grid_latitude=lambda v: -18. <= v <= -8.) 

        CCCma = CCCma.extract(Malawi) 

        #time constraint to make all series the same
        iris.FUTURE.cell_datetime_objects = True
        t_constraint = iris.Constraint(time=lambda cell: 1989 <= cell.point.year <= 2008)

        CCCma = CCCma.extract(t_constraint)

        #Convert units to match, CORDEX data is in Kelvin but Observed data in Celsius, we would like to show all data in Celsius
        CCCma.convert_units('Celsius')

        #plot map with physical features
        cmap = plt.cm.afmhot_r    
        ax = plt.axes(projection=cartopy.crs.PlateCarree())
        ax.add_feature(cartopy.feature.COASTLINE)   
        ax.add_feature(cartopy.feature.BORDERS)
        ax.add_feature(cartopy.feature.LAKES, alpha=0.5)
        ax.add_feature(cartopy.feature.RIVERS)

        #set map boundary
        ax.set_extent([31, 36.5, -8,-18])

        #set axis tick marks
        ax.set_xticks([32, 34, 36])
        ax.set_yticks([-9, -11, -13, -15, -17])
        lon_formatter = LongitudeFormatter(zero_direction_label=True)
        lat_formatter = LatitudeFormatter()
        ax.xaxis.set_major_formatter(lon_formatter)
        ax.yaxis.set_major_formatter(lat_formatter)
        data = CCCma

        #take mean of data over all time
        plot = iplt.contourf(data.collapsed('time', iris.analysis.MEAN), \
                  cmap=cmap, levels=[15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31],\
                  extend='both')

        #add colour bar index
        plt.colorbar(plot)

        #give map a title
        plt.title('RCP4.5 Mean Temperature 1989-2008', fontsize=10)

        plt.show()

    if __name__ == '__main__':
        main()

How can I amend this to take the difference between the two datasets? I tried this code:

        import matplotlib.pyplot as plt
        import iris
        import iris.plot as iplt
        import cartopy
        from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter
        import iris.analysis.cartography


        #this file is split into parts as follows:
            #PART 1: load and format CORDEX models
            #PART 2: load and format observed data
            #PART 3: format data
            #PART 4: plot data

        def main():
            #PART 1: CORDEX MODELS
            #bring in all the models we need and give them a name
            CCCma = '/exports/csce/datastore/geos/users/s0XXXX/Climate_Modelling/AFR_44_tas/ERAINT/1979-2012/tas_AFR-44_ECMWF-ERAINT_evaluation_r1i1p1_CCCma-CanRCM4_r2_mon_198901-200912.nc'

            #Load exactly one cube from given file    
            CCCma = iris.load_cube(CCCma)    


            #we are only interested in the latitude and longitude relevant to Malawi   
            Malawi = iris.Constraint(grid_longitude=lambda v: 31 <= v <= 36.5, \
                                 grid_latitude=lambda v: -18. <= v <= -8.) 

            CCCma = CCCma.extract(Malawi) 


            #time constraint to make all series the same
            iris.FUTURE.cell_datetime_objects = True
            t_constraint = iris.Constraint(time=lambda cell: 1989 <= cell.point.year <= 2008)

            CCCma = CCCma.extract(t_constraint)


            #PART 2: OBSERVED DATA
            #bring in all the files we need and give them a name
            CRU= '/exports/csce/datastore/geos/users/s0XXXX/Climate_Modelling/Actual_Data/cru_ts4.00.1901.2015.tmp.dat.nc'

            #Load exactly one cube from given file
            CRU = iris.load_cube(CRU, 'near-surface temperature')

            #we are only interested in the latitude and longitude relevant to Malawi 
            Malawi = iris.Constraint(longitude=lambda v: 32.5 <= v <= 36., \
                                 latitude=lambda v: -17. <= v <= -9.) 
            CRU = CRU.extract(Malawi)


            #time constraint to make all series the same
            iris.FUTURE.cell_datetime_objects = True
            t_constraint = iris.Constraint(time=lambda cell: 1989 <= cell.point.year <= 2008)

            CRU = CRU.extract(t_constraint)


            #PART 3: FORMAT DATA

            #Convert units to match
            CCCma.convert_units('Celsius')
            CRU.convert_units('Celsius')

            #Take difference between two datasets
            Bias_CCCma = CCCma-CRU

            #PART 4: PLOT MAP
            #plot map with physical features
            cmap = plt.cm.afmhot_r    
            ax = plt.axes(projection=cartopy.crs.PlateCarree())
            ax.add_feature(cartopy.feature.COASTLINE)   
            ax.add_feature(cartopy.feature.BORDERS)
            ax.add_feature(cartopy.feature.LAKES, alpha=0.5)
            ax.add_feature(cartopy.feature.RIVERS)

            #set map boundary
            ax.set_extent([31, 36.5, -8,-18])

            #set axis tick marks
            ax.set_xticks([32, 34, 36])
            ax.set_yticks([-9, -11, -13, -15, -17])
            lon_formatter = LongitudeFormatter(zero_direction_label=True)
            lat_formatter = LatitudeFormatter()
            ax.xaxis.set_major_formatter(lon_formatter)
            ax.yaxis.set_major_formatter(lat_formatter)
            data = Bias_CCCma

            #take mean of data over all time
            plot = iplt.contourf(data.collapsed('time', iris.analysis.MEAN), \
                      cmap=cmap, levels=[15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31],\
                      extend='both')

            #add colour bar index
            plt.colorbar(plot)

            #give map a title
            plt.title('RCP4.5 Mean Temperature 1989-2008', fontsize=10)

            plt.show()

        if __name__ == '__main__':
            main()

However this gives me the following error:

ValueError: This operation cannot be performed as there are differing coordinates (grid_latitude, grid_longitude, time) remaining which cannot be ignored.

I was pretty sure this wasn't going to be so simple, but I'm not sure how to fix it. Any ideas? TIA!

Answer 1

My guess is that CCCma and CRU are on different grids, so when you try to subtract them you get an error. You probably need to interpolate them to the same grid first (otherwise, how would iris know which grid you want the result to lie on?).

Answer 2

Iris is very strict about matching up the cube coordinates for binary operations and there is an open issue discussing whether and how to make it more flexible ready for version 2. In the meantime, if your cubes are the same shape and you don't mind loading the data, you could just do

Bias_CCCma = CCCma - CRU.data

If your cubes are different shapes (i.e. the models are on different grids, as Jeremy suggested) or you don't want to load the data, there are a few things to look at:

1. If the grids are different then you will need to regrid one of the cubes to match the other.
2. For the subtraction operation, the grid coordinate names need to match up. If you are confident that grid_latitude and grid_longitude mean the same as latitude and longitude, you can rename the grid coordinates on one of your cubes. You will also need to ensure the other coordinate metadata matches (e.g. var_name is often an issue).
3. The time coordinate coming up in your error message is almost certainly due to the unit mismatch you identified in your previous question. I think this issue should go away if you reorder your code to do the time averaging first and then take the difference (the binary operations don't care so much about scalar coordinates).

Answer 3

Thank you all for your answers. In the end I needed to regrid the data first, as @RuthC suggested.

So the code changed to look like this:

    import matplotlib.pyplot as plt
    import matplotlib.cm as mpl_cm
    import numpy as np
    from cf_units import Unit
    import iris
    import iris.plot as iplt
    import cartopy
    from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter
    import iris.analysis.cartography
    import iris.coord_categorisation as iriscc

    #this file is split into parts as follows:
        #PART 1: load and format CORDEX models
        #PART 2: load and format observed data
        #PART 3: format data
        #PART 4: plot data

def main():
    iris.FUTURE.netcdf_promote=True    

    #PART 1: CORDEX MODELS
    #bring in all the models we need and give them a name
    CCCma = '/exports/csce/datastore/geos/users/s0899345/Climate_Modelling/AFR_44_tasmax/ERAINT/1979-2012/tasmax_AFR-44_ECMWF-ERAINT_evaluation_r1i1p1_CCCma-CanRCM4_r2_mon_198901-200912.nc'

    #Load exactly one cube from given file    
    CCCma = iris.load_cube(CCCma)     


    #remove flat latitude and longitude and only use grid latitude and grid longitude to make consistent with the observed data, also make sure all of the longitudes are monotonic 
    lats = iris.coords.DimCoord(CORDEX.coord('latitude').points[:,0], \
                            standard_name='latitude', units='degrees')
    lons = CORDEX.coord('longitude').points[0]
    for i in range(len(lons)):
         if lons[i]>100.:
             lons[i] = lons[i]-360.
    lons = iris.coords.DimCoord(lons, \
                            standard_name='longitude', units='degrees')

    CORDEX.remove_coord('latitude')
    CORDEX.remove_coord('longitude')
    CORDEX.remove_coord('grid_latitude')
    CORDEX.remove_coord('grid_longitude')
    CORDEX.add_dim_coord(lats, 1)
    CORDEX.add_dim_coord(lons, 2)

    #PART 2: OBSERVED DATA
    #bring in all the files we need and give them a name
    CRU= '/exports/csce/datastore/geos/users/s0XXXX/Climate_Modelling/Actual_Data/cru_ts4.00.1901.2015.tmp.dat.nc'

    #Load exactly one cube from given file
    CRU = iris.load_cube(CRU, 'near-surface temperature')

    #PART 3: FORMAT DATA
    #Regrid observed data onto rotated pole grid
    CRU = CRU.regrid(CORDEX, iris.analysis.Linear())

    #we are only interested in the latitude and longitude relevant to Malawi 
    Malawi = iris.Constraint(longitude=lambda v: 32.5 <= v <= 36.5, \
                     latitude=lambda v: -17. <= v <= -9.) 

    CORDEX = CORDEX.extract(Malawi) 
    CRU = CRU.extract(Malawi)

    #time constraint to make all series the same
    iris.FUTURE.cell_datetime_objects = True
    t_constraint = iris.Constraint(time=lambda cell: 1990<= cell.point.year <= 2008)

    CORDEX = CORDEX.extract(t_constraint)
    CRU = CRU.extract(t_constraint)

    #Convert units to match
    CORDEX.convert_units('Celsius')
    CRU.unit = Unit('Celsius') # This fixes CRU which is in 'Degrees Celsius' to read 'Celsius'

    #add years to data
    iriscc.add_year(CORDEX, 'time')
    iriscc.add_year(CRU, 'time')

    #We are interested in plotting the data for the average of the time period.
    CORDEX = CORDEX.collapsed('time', iris.analysis.MEAN)
    CRU = CRU.collapsed(['time'], iris.analysis.MEAN)

    #Take difference between two datasets
    Bias = CRU-CORDEX

    #PART 4: PLOT MAP
    #load color palette
    colourA = mpl_cm.get_cmap('brewer_YlOrRd_09')

    #plot map with physical features 
    ax = plt.axes(projection=cartopy.crs.PlateCarree())
    ax.add_feature(cartopy.feature.COASTLINE)   
    ax.add_feature(cartopy.feature.BORDERS)
    ax.add_feature(cartopy.feature.LAKES, alpha=0.5)
    ax.add_feature(cartopy.feature.RIVERS)

    #set map boundary
    ax.set_extent([32.5, 36., -9, -17]) 

    #set axis tick marks
    ax.set_xticks([33, 34, 35]) 
    ax.set_yticks([-10, -12, -14, -16]) 
    lon_formatter = LongitudeFormatter(zero_direction_label=True)
    lat_formatter = LatitudeFormatter()
    ax.xaxis.set_major_formatter(lon_formatter)
    ax.yaxis.set_major_formatter(lat_formatter)

    #plot data and set colour range
    plot = iplt.contourf(CORDEX, cmap=colourA, levels=np.arange(13,32,1), extend='both')

    #add colour bar index and a label
    plt.colorbar(plot, label='Celsius')

    #give map a title
    plt.title('Tasmax 1990-2008 - CanRCM4_ERAINT ', fontsize=10)

    #save the image of the graph and include full legend
    plt.savefig('ERAINT_CCCma_Tasmax_MAP_Annual', bbox_inches='tight')

    plt.show()

    if __name__ == '__main__':
        main()