Python Matplotlib 两个 NetCDF 数据集之间的差异

Question

我正在尝试绘制特定地理区域内气候模拟数据和观测数据之间的差异。

要创建气候模拟地图，我正在使用此代码

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()

如何修改它以获取两个数据集之间的差异？我试过这段代码：

        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()

但是这给了我以下错误：

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

我很确定这不会那么简单，但我不知道如何解决它。有任何想法吗？蒂亚！

Answer 1

我的猜测是 CCCma 和 CRU 位于不同的网格上，因此当您尝试减去它们时会出现错误。您可能需要先将它们插入到同一个网格中（否则，iris 怎么知道您希望结果位于哪个网格上？）。

Answer 2

Iris 对于二元运算的立方体坐标匹配非常严格，有一个open issue 讨论是否以及如何使其更灵活地为版本 2 做好准备。与此同时，如果你的立方体形状相同，而你不介意加载数据，你可以这样做

Bias_CCCma = CCCma - CRU.data

如果您的立方体是不同的形状（即模型位于不同的网格上，正如 Jeremy 建议的那样）或者您不想加载数据，则需要注意以下几点：

1. 如果网格不同，那么您需要regrid 一个立方体来匹配另一个。
2. 对于减法运算，网格坐标名称需要匹配。如果您确信grid_latitude 和grid_longitude 与latitude 和longitude 的含义相同，那么您可以rename 您的一个立方体上的网格坐标。您还需要确保其他坐标元数据匹配（例如，var_name 通常是一个问题）。
3. 错误消息中出现的time 坐标几乎可以肯定是由于您在previous question 中确定的单位不匹配。我认为如果您重新排序代码以先进行时间平均然后取差值（二进制运算不太关心标量坐标），这个问题应该会消失。

Answer 3

谢谢大家的回答。最后，正如@RuthC 建议的那样，我需要先重新网格化数据。

所以代码变成了这样：

    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()