<jupyter_start><jupyter_text>Handling netCDF files for simple climate modelsIn this notebook we give a brief introduction to iris, the library we use for our analysis, before giving a demonstration of some of the key functionality of `netcdf_scm`.<jupyter_code># NBVAL_IGNORE_OUTPUT
from os.path import join
    
import numpy as np
import iris
import iris.coord_categorisation
import iris.quickplot as qplt
import iris.analysis.cartography
import matplotlib
import matplotlib.pyplot as plt

from netcdf_scm.iris_cube_wrappers import SCMCube, MarbleCMIP5Cube
plt.style.use('bmh') 
%matplotlib inline
DATA_PATH_TEST = join("..", "tests", "test-data")
DATA_PATH_TEST_MARBLE_CMIP5_ROOT = join(DATA_PATH_TEST, "marble-cmip5")<jupyter_output><empty_output><jupyter_text>Loading a cube Loading with irisHere we show how to load a cube directly using iris.<jupyter_code>tas_file = join(
    DATA_PATH_TEST_MARBLE_CMIP5_ROOT,
    "cmip5",
    "1pctCO2",
    "Amon",
    "tas",
    "CanESM2",
    "r1i1p1",
    "tas_Amon_CanESM2_1pctCO2_r1i1p1_189201-190312.nc",
)
# NBVAL_IGNORE_OUTPUT
# Ignore output as the warnings are likely to change with
# new iris versions
tas_iris_load = SCMCube()
# you need this in case your cube has multiple variables
variable_constraint = iris.Constraint(
    cube_func=(lambda c: c.var_name == np.str("tas"))
)

tas_iris_load.cube = iris.load_cube(tas_file, variable_constraint)<jupyter_output>Using or importing the ABCs from 'collections' instead of from 'collections.abc' is deprecated, and in 3.8 it will stop working
Missing CF-netCDF measure variable 'areacella', referenced by netCDF variable 'tas'
Using or importing the ABCs from 'collections' instead of from 'collections.abc' is deprecated, and in 3.8 it will stop working<jupyter_text>The warning tells us that we need to add the areacella as a measure variable to our cube. Doing this manually everytime involves finding the areacella file, loading it, turning it into a cell measure and then adding it to the cube. This is a pain and involves about 100 lines of code. To make life easier, we wrap all of that away using `netcdf_scm`, which we will introduce in the next section.  Loading with `netcdf_scm`There are a couple of particularly annoying things involved with processing netCDF data. Firstly, the data is often stored in a folder hierarchy which can be fiddly to navigate. Secondly, the metadata is often stored separate from the variable cubes. Hence in `netcdf_scm`, we try to abstract the code which solves these two things away to make life a bit easier. This involves defining a cube in `netcdf_scm.iris_cube_wrappers`. The details can be read there, for now we just give an example.Our example uses `MarbleCMIP5Cube`. This cube is designed to work with the CMIP5 data on our server at University of Melbourne, which has been organised into a number of folders which are similar, but not quite identical, to the CMOR directory structure described in section 3.1 of the [CMIP5 Data Reference Syntax](https://cmip.llnl.gov/cmip5/docs/cmip5_data_reference_syntax_v1-00_clean.pdf). To facilitate our example, the test data in `DATA_PATH_TEST_MARBLE_CMIP5_ROOT` is organised in the same way. Loading with identifiersWith our `MarbleCMIP5Cube`, we can simply pass in the information about the data we want (experiment, model, ensemble member etc.) and it will load our desired cube using the `load_data_from_identifiers` method.<jupyter_code>tas = MarbleCMIP5Cube()
tas.load_data_from_identifiers(
    root_dir=DATA_PATH_TEST_MARBLE_CMIP5_ROOT,
    activity="cmip5",
    experiment="1pctCO2",
    mip_table="Amon",
    variable_name="tas",
    model="CanESM2",
    ensemble_member="r1i1p1",
    time_period="189201-190312",
    file_ext=".nc",
)<jupyter_output><empty_output><jupyter_text>We can verify that the loaded cube is exactly the same as the cube we loaded in the previous section (where we provided the full path).<jupyter_code># NBVAL_IGNORE_OUTPUT
assert tas.cube == tas_iris_load.cube<jupyter_output><empty_output><jupyter_text>We can have a look at our cube too (note that the broken cell measures representation is intended to be fixed in https://github.com/SciTools/iris/pull/3173).<jupyter_code>tas.cube<jupyter_output><empty_output><jupyter_text>Loading with filepathWith our `MarbleCMIP5Cube`, we can also pass in the filepath and the cube will determine the relevant attributes for us, as well as loading the other required cubes.<jupyter_code>example_path = join(
    DATA_PATH_TEST_MARBLE_CMIP5_ROOT,
    "cmip5",
    "1pctCO2",
    "Amon",
    "tas",
    "CanESM2",
    "r1i1p1",
    "tas_Amon_CanESM2_1pctCO2_r1i1p1_189201-190312.nc",
)
example_path
tas_from_path = MarbleCMIP5Cube()
tas_from_path.load_data_from_path(example_path)
tas_from_path.model<jupyter_output><empty_output><jupyter_text>We can also confirm that this cube is the same again.<jupyter_code># NBVAL_IGNORE_OUTPUT
assert tas_from_path.cube == tas_iris_load.cube<jupyter_output><empty_output><jupyter_text>Acting on the cubeOnce we have loaded our `SCMCube`, we can act on its `cube` attribute like any other Iris Cube. For example, we can add a year categorisation, take an annual mean and then plot the timeseries.<jupyter_code># NBVAL_IGNORE_OUTPUT
year_cube = tas.cube.copy()
iris.coord_categorisation.add_year(year_cube, 'time')
annual_mean = year_cube.aggregated_by(
    ['year'], # Do the averaging 
    iris.analysis.MEAN
)
global_annual_mean = annual_mean.collapsed(
    ['latitude', 'longitude'], 
    iris.analysis.MEAN,
    weights=iris.analysis.cartography.area_weights(annual_mean)
)
qplt.plot(global_annual_mean);<jupyter_output>Using DEFAULT_SPHERICAL_EARTH_RADIUS.
Using or importing the ABCs from 'collections' instead of from 'collections.abc' is deprecated, and in 3.8 it will stop working<jupyter_text>We can also take a time average and make a spatial plot.<jupyter_code># NBVAL_IGNORE_OUTPUT
time_mean = tas.cube.collapsed('time', iris.analysis.MEAN)
qplt.contourf(
    time_mean, 
    levels=np.arange(200, 320, 10), 
    extend='max',
);<jupyter_output><empty_output><jupyter_text>Add a spatial plot with coastlines.<jupyter_code># NBVAL_IGNORE_OUTPUT
# we ignore output here as CI sometimes has to 
# download the map file
qplt.contourf(
    time_mean, 
    levels=np.arange(200, 320, 10), 
    extend='max',
)
plt.gca().coastlines();<jupyter_output><empty_output><jupyter_text>SCM specificsFinally, we present the key functions of this package. These are directly related to processing netCDF files for simple climate models. Getting SCM timeseriesThe major one is `get_scm_timeseries`. This function wraps a number of steps:1. load the land surface fraction data1. combine the land surface fraction and latitude data to determine the hemisphere and land/ocean boxes1. cut the data into the relevant boxes1. take a time average in each box1. return it all as an OpenSCMDataframeAs you can imagine, we find it very useful to be able to abstract all these normally nasty steps away.<jupyter_code>tas_scm_timeseries = tas.get_scm_timeseries()
type(tas_scm_timeseries)
tas_scm_timeseries.tail()<jupyter_output><empty_output><jupyter_text>Having the data as a Dataframe makes it trivial to plot and work with.<jupyter_code>restricted_time_df = tas_scm_timeseries.filter(
    year=range(1895, 1901),
    region="*Ocean*"  # Try e.g. "*", "World", "*Land", "*Southern Hemisphere*" here
)
restricted_time_df.line_plot(
    x="time", 
    color="region", 
    figsize=(16, 9),
);
tas_scm_timeseries_annual_mean = tas_scm_timeseries.filter(
    region="World"
).timeseries().T
tas_scm_timeseries_annual_mean = tas_scm_timeseries_annual_mean.groupby(
    tas_scm_timeseries_annual_mean.index.map(lambda x: x.year)
).mean()
tas_scm_timeseries_annual_mean.head()
tas_scm_timeseries_annual_mean.plot(figsize=(16, 9));<jupyter_output><empty_output><jupyter_text>Getting SCM timeseries cubesAs part of the process above, we calculate all the timeseries as `iris.cube.Cube`'s. Extracting these intermediate cubes can be done with `get_scm_timeseries_cubes`. These intermediate cubes are useful as they contain all the metadata from the source cube in a slightly more friendly format than `SCMDataFrame`'s `metadata` attribute [Note: `SCMDataFrame`'s metadata handling is a work in progress].<jupyter_code>tas_scm_ts_cubes = tas.get_scm_timeseries_cubes()
# NBVAL_IGNORE_OUTPUT
print(tas_scm_ts_cubes["World"].cube)<jupyter_output>air_temperature                            (time: 144)
     Dimension coordinates:
          time                                  x
     Scalar coordinates:
          height: 2.0 m
          land_fraction: 0.30309096774652516
          land_fraction_northern_hemisphere: 0.4076169376537766
          land_fraction_southern_hemisphere: 0.19856502699902237
          latitude: 0.0 degrees, bound=(-90.0, 90.0) degrees
          longitude: 178.59375 degrees, bound=(-1.40625, 358.59375) degrees
    Scalar cell measures:
          cell_area
     Attributes:
          CCCma_data_licence: 1) GRANT OF LICENCE - The Government of Canada (Environment Canada) is...
          CCCma_parent_runid: IGA
          CCCma_runid: IDK
          CDI: Climate Data Interface version 1.9.7.1 (http://mpimet.mpg.de/cdi)
          CDO: Climate Data Operators version 1.9.7.1 (http://mpimet.mpg.de/cdo)
          Conventions: CF-1.4
          activity_id: cmip5
          associated_files: baseURL: http://cmip-pcmdi.lln[...]<jupyter_text>In particular, the `land_fraction*` auxillary co-ordinates provide useful information about the fraction of area that was assumed to be land in the crunching.<jupyter_code>tas_scm_ts_cubes["World"].cube.coords("land_fraction")
tas_scm_ts_cubes["World"].cube.coords("land_fraction_northern_hemisphere")<jupyter_output><empty_output><jupyter_text>Investigating the weightsAnother utility function is `get_scm_timeseries_weights`. This function is very similar to `get_scm_timeseries` but returns just the weights rather than a Dataframe. These weights are also area weighted.<jupyter_code>tas_scm_weights = tas.get_scm_timeseries_weights()
# NBVAL_IGNORE_OUTPUT
plt.figure(figsize=(18, 18))
no_rows = 3
no_cols = 4

total_panels = no_cols * no_rows
rows_plt_comp = no_rows*100
cols_plt_comp = no_cols*10
for i, (label, weights) in enumerate(tas_scm_weights.items()):
    if label == "World":
        index = (int((no_rows + 1) / 2))
        plt.subplot(no_rows, 1, index)
    else:
        if label == "World|Northern Hemisphere":
            index = 1
        elif label == "World|Southern Hemisphere":
            index = 1 + (no_rows - 1) * no_cols
        elif label == "World|Land":
            index = 2
        elif label == "World|Ocean":
            index = 2 + (no_rows - 1) * no_cols
        else:
            index = 3
            if "Ocean" in label:
                index += 1
            if "Southern Hemisphere" in label:
                index += (no_rows - 1) * no_cols

        plt.subplot(no_rows, no_cols, index)
        
    weight_cube = tas.cube.collapsed("time", iris.analysis.MEAN)
    weight_cube.data = weights
    qplt.contourf(
        weight_cube, 
        extend='max',
    )
    plt.title(label)
    plt.gca().coastlines()
    

plt.tight_layout()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.08254842722587524, 'fraction_numerical': 0.02319175404300693, 'mean_word_length': 4.460941290635613, 'pattern_counts': {'":': 5, '<': 58, '<?xml version=': 0, '>': 58, 'https://': 2, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 2, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '15', 'dir': 'jupyter-structured-clean-dedup', 'id': '154404', 'license': 'BSD-2-Clause', 'n_tokens_mistral': 3974, 'n_tokens_neox': 3694, 'n_words': 1097, 'path': 'notebooks/demo.ipynb', 'repo_name': 'lewisjared/netcdf-scm', 'size': '11254'}