1. 首頁
  2. 問答
  3. Python Matplotlib兩個NetCDF數據集之間的區別

Python Matplotlib兩個NetCDF數據集之間的區別

2017-09-18 34 views
0

我試圖在設定的地理區域上映射氣候模擬數據和觀測數據之間的差異。Python Matplotlib兩個NetCDF數據集之間的區別

只創建了氣候模擬的地圖，我用這個代碼

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.

我很肯定這是不會這麼簡單，但我不知道如何解決它。有任何想法嗎？ TIA！

來源

2017-09-18 ErikaAWT

+0

哪一行代碼爲您提供了錯誤？你能寫一個更簡單的腳本形式嗎？它給出了同樣的錯誤，我們可以在本地運行？ –

回答

0

我的猜測是CCCma和CRU在不同的網格上，所以當你試圖減去它們時，你會得到一個錯誤。你可能需要先插入它們到同一個網格（否則，iris怎麼知道你想要結果放在哪個網格上？）。

來源

2017-09-18 17:03:37

0

虹膜是非常嚴格的配對二進制操作的立方體座標，並且有一個open issue討論是否以及如何使它更靈活地準備版本2.同時，如果你的立方體是相同的形狀，你不介意加載數據，你可以只是做

Bias_CCCma = CCCma - CRU.data

如果你的立方體形狀不同（即模型在不同的網格，傑里米建議），或者你不希望加載數據，有有幾件事要看：

  1. 如果網格是不同的，那麼你將需要regrid其中一個立方體匹配其他。
  2. 對於減法運算，網格座標名稱需要匹配。如果您確信grid_latitude和grid_longitude與latitude和longitude的含義相同，您可以在其中一個立方體上使用網格座標rename。您還需要確保其他座標元數據匹配（例如var_name通常是個問題）。
  3. time座標出現在您的錯誤信息幾乎肯定是由於您在previous question中確定的單位不匹配。我認爲這個問題應該消失，如果你重新排序你的代碼首先進行時間平均，然後消除差異（二進制操作不關心標量座標）。

來源

2017-09-18 18:11:56 RuthC

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-03 11:53:14 ErikaAWT

相關問題

 相關問題