Updated use case

notebooks/demo/use-case_simple-vis_xarray-matplotlib_cmip6.ipynb

@@ -15,7 +15,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "filename = '/pool/data/CMIP6/data/CMIP/NCAR/CESM2/historical/r1i1p1f1/Amon/tas/gn/v20190308/tas_Amon_CESM2_historical_r1i1p1f1_gn_185001-201412.nc'"
+    "filename = '/work/ik1017/CMIP6/data/CMIP6/CMIP/NCAR/CESM2/historical/r1i1p1f1/Amon/tas/gn/v20190308/tas_Amon_CESM2_historical_r1i1p1f1_gn_185001-201412.nc'"
    ]
   },
   {
 %% Cell type:markdown id: tags:
 
 # Simple visualisation of a CMIP6 dataset with Python Xarray
 
   - this is an adoption of [NordicESMhub example](https://nordicesmhub.github.io/NEGI-Abisko-2019/training/CMIP6_example.html) to use CMIP6 data stored in the DKRZ data pool
 
 %% Cell type:code id: tags:
 
 ``` python
-filename = '/pool/data/CMIP6/data/CMIP/NCAR/CESM2/historical/r1i1p1f1/Amon/tas/gn/v20190308/tas_Amon_CESM2_historical_r1i1p1f1_gn_185001-201412.nc'
+filename = '/work/ik1017/CMIP6/data/CMIP6/CMIP/NCAR/CESM2/historical/r1i1p1f1/Amon/tas/gn/v20190308/tas_Amon_CESM2_historical_r1i1p1f1_gn_185001-201412.nc'
 ```
 
 %% Cell type:markdown id: tags:
 
 ### Import python packages
 
 %% Cell type:code id: tags:
 
 ``` python
 import xarray as xr
 import cartopy.crs as ccrs
 import matplotlib.pyplot as plt
 import numpy as np
 import cftime
 ```
 
 %% Cell type:markdown id: tags:
 
 ## Open dataset
 
 - Use `xarray` python package to analyze netCDF dataset
 - `open_dataset` allows to get all the metadata without loading data into memory.
 - with `xarray`, we only load into memory what is needed.
 
 %% Cell type:code id: tags:
 
 ``` python
 dset = xr.open_dataset(filename, decode_times=True, use_cftime=True)
 print(dset)
 ```
 
 %% Cell type:markdown id: tags:
 
 ### Get metadata corresponding to near-surface air temperature (tas)
 
 %% Cell type:code id: tags:
 
 ``` python
 print(dset['tas'])
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 dset.time.values
 ```
 
 %% Cell type:markdown id: tags:
 
 ### Select time
 
 - Select a specific time
 
 %% Cell type:code id: tags:
 
 ``` python
 dset['tas'].sel(time=cftime.DatetimeNoLeap(1850, 1, 15, 12, 0, 0, 0, 2, 15)).plot(cmap = 'coolwarm')
 ```
 
 %% Cell type:markdown id: tags:
 
 - select the nearest time. Here from 1st April 1950
 
 %% Cell type:code id: tags:
 
 ``` python
 dset['tas'].sel(time=cftime.DatetimeNoLeap(1850, 4, 1), method='nearest').plot(cmap='coolwarm')
 ```
 
 %% Cell type:markdown id: tags:
 
 ## Customize plot
 
 ### Set the size of the figure and add coastlines
 
 %% Cell type:code id: tags:
 
 ``` python
 fig = plt.figure(1, figsize=[30,13])
 
 # Set the projection to use for plotting
 ax = plt.subplot(1, 1, 1, projection=ccrs.PlateCarree())
 ax.coastlines()
 
 # Pass ax as an argument when plotting. Here we assume data is in the same coordinate reference system than the projection chosen for plotting
 # isel allows to select by indices instead of the time values
 dset['tas'].isel(time=0).plot.pcolormesh(ax=ax, cmap='coolwarm')
 ```
 
 %% Cell type:markdown id: tags:
 
 ### Change plotting projection
 
 %% Cell type:code id: tags:
 
 ``` python
 fig = plt.figure(1, figsize=[10,10])
 
 # We're using cartopy and are plotting in Orthographic projection
 # (see documentation on cartopy)
 ax = plt.subplot(1, 1, 1, projection=ccrs.Orthographic(0, 90))
 ax.coastlines()
 
 # We need to project our data to the new Orthographic projection and for this we use `transform`.
 # we set the original data projection in transform (here PlateCarree)
 dset['tas'].isel(time=0).plot(ax=ax, transform=ccrs.PlateCarree(), cmap='coolwarm')
 
 # One way to customize your title
 plt.title(dset.time.values[0].strftime("%B %Y"), fontsize=18)
 ```
 
 %% Cell type:markdown id: tags:
 
 ### Choose the extent of values
 - Fix your minimum and maximum values in your plot and
 - Use extend so values below the minimum and max
 
 %% Cell type:code id: tags:
 
 ``` python
 fig = plt.figure(1, figsize=[10,10])
 
 ax = plt.subplot(1, 1, 1, projection=ccrs.Orthographic(0, 90))
 ax.coastlines()
 
 # Fix extent
 minval = 240
 maxval = 300
 
 # pass extent with vmin and vmax parameters
 dset['tas'].isel(time=0).plot(ax=ax, vmin=minval, vmax=maxval, transform=ccrs.PlateCarree(), cmap='coolwarm')
 
 # One way to customize your title
 plt.title(dset.time.values[0].strftime("%B %Y"), fontsize=18)
 ```
 
 %% Cell type:markdown id: tags:
 
 ## Multiplots
 
 ### Faceting
 
 %% Cell type:code id: tags:
 
 ``` python
 proj_plot = ccrs.Orthographic(0, 90)
 
 p = dset['tas'].sel(time = dset.time.dt.year.isin([1850, 2014])).plot(x='lon', y='lat',
                                                                  transform=ccrs.PlateCarree(),
                                                                  aspect=dset.dims["lon"] / dset.dims["lat"],  # for a sensible figsize
                                                                  subplot_kws={"projection": proj_plot},
                                                                  col='time', col_wrap=6, robust=True, cmap='PiYG')
 # We have to set the map's options on all four axes
 for ax,i in zip(p.axes.flat,  dset.time.sel(time = dset.time.dt.year.isin([1850, 2014])).values):
     ax.coastlines()
     ax.set_title(i.strftime("%B %Y"), fontsize=18)
 ```
 
 %% Cell type:markdown id: tags:
 
 ### Combine plots with different projections
 
 %% Cell type:code id: tags:
 
 ``` python
 fig = plt.figure(1, figsize=[20,10])
 
 # Fix extent
 minval = 240
 maxval = 300
 
 # Plot 1 for Northern Hemisphere subplot argument (nrows, ncols, nplot)
 # here 1 row, 2 columns and 1st plot
 ax1 = plt.subplot(1, 2, 1, projection=ccrs.Orthographic(0, 90))
 
 # Plot 2 for Southern Hemisphere
 # 2nd plot
 ax2 = plt.subplot(1, 2, 2, projection=ccrs.Orthographic(180, -90))
 
 tsel = 0
 for ax,t in zip([ax1, ax2], ["Northern", "Southern"]):
     map = dset['tas'].isel(time=tsel).plot(ax=ax, vmin=minval, vmax=maxval,
                                            transform=ccrs.PlateCarree(),
                                            cmap='coolwarm',
                                            add_colorbar=False)
     ax.set_title(t + " Hemisphere \n" , fontsize=15)
     ax.coastlines()
     ax.gridlines()
 
 # Title for both plots
 fig.suptitle('Near Surface Temperature\n' + dset.time.values[tsel].strftime("%B %Y"), fontsize=20)
 
 
 cb_ax = fig.add_axes([0.325, 0.05, 0.4, 0.04])
 
 cbar = plt.colorbar(map, cax=cb_ax, extend='both', orientation='horizontal', fraction=0.046, pad=0.04)
 cbar.ax.tick_params(labelsize=25)
 cbar.ax.set_ylabel('K', fontsize=25)
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 ```