0

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

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

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.

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

4

3 回答 3

0

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

Bias_CCCma = CCCma - CRU.data

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

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

我的猜测是 CCCma 和 CRU 在不同的网格上，所以当你试图减去它们时，你会得到一个错误。您可能需要先将它们插入同一个网格（否则，如何iris知道您希望结果位于哪个网格上？）。

于 2017-09-18T17:03:37.420 回答
0

谢谢大家的答案。最后，正如@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()
于 2017-10-03T11:53:14.627 回答