Skip to content

Instantly share code, notes, and snippets.

@pelson

pelson/cube_creation.ipynb

Created February 14, 2017 01:32
Show Gist options
  • Save pelson/1b717fc77e6186357c9c03c4cc91828f to your computer and use it in GitHub Desktop.
Save pelson/1b717fc77e6186357c9c03c4cc91828f to your computer and use it in GitHub Desktop.
Download ZIP
How to create an Iris cube by hand, save to NetCDF, and manipulate time dimensions
Display the source blob
Display the rendered blob
Raw
cube_creation.ipynb
{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Q: How does one create an Iris cube from scratch and save it as a NetCDF file?"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"I was asked specifically if it was possible to make a cube of orthogonal time components (year, month, day), so I'll start there:"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"import iris.cube\n",
"import iris.coords"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"import numpy as np"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"First, we need an array of data..."
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"data = np.empty(shape=[30, 12, 4, 500, 1000], dtype=np.float64)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now we create the cube to contain that data:"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"unknown / (unknown) (-- : 30; -- : 12; -- : 4; -- : 500; -- : 1000)\n"
]
}
],
"source": [
"cube = iris.cube.Cube(data)\n",
"print(cube)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Clearly this dataset needs some metadata..."
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"air_temperature / (kelvin) (-- : 30; -- : 12; -- : 4; -- : 500; -- : 1000)\n",
" Attributes:\n",
" source: Data from xyz\n"
]
}
],
"source": [
"cube.rename('air_temperature')\n",
"cube.units = 'kelvin'\n",
"cube.attributes['source'] = 'Data from xyz'\n",
"print(cube)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"And some coordinate metadata:"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"air_temperature / (kelvin) (year: 30; -- : 12; -- : 4; -- : 500; -- : 1000)\n",
" Dimension coordinates:\n",
" year x - - - -\n",
" Attributes:\n",
" source: Data from xyz\n"
]
}
],
"source": [
"year = iris.coords.DimCoord(points=np.arange(1980, 2010),\n",
" long_name='year', units='1')\n",
"cube.add_dim_coord(year, data_dim=0)\n",
"print(cube)"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"month = iris.coords.DimCoord(points=np.arange(1, 13),\n",
" long_name='month', units='1')\n",
"cube.add_dim_coord(month, data_dim=1)\n",
"\n",
"day = iris.coords.DimCoord(points=[1, 9, 17, 25],\n",
" long_name='day', units='1')\n",
"cube.add_dim_coord(day, data_dim=2)"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"lat = iris.coords.DimCoord(points=np.linspace(-90, 90, 500),\n",
" standard_name='latitude', units='degrees')\n",
"cube.add_dim_coord(lat, data_dim=3)\n",
"\n",
"lon = iris.coords.DimCoord(points=np.linspace(-180, 180, 1000),\n",
" standard_name='longitude', units='degrees')\n",
"cube.add_dim_coord(lon, data_dim=4)"
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"air_temperature / (kelvin) (year: 30; month: 12; day: 4; latitude: 500; longitude: 1000)\n",
" Dimension coordinates:\n",
" year x - - - -\n",
" month - x - - -\n",
" day - - x - -\n",
" latitude - - - x -\n",
" longitude - - - - x\n",
" Attributes:\n",
" source: Data from xyz\n"
]
}
],
"source": [
"print(cube)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We may also want to combine the date information into a single auxilliary time coordinate. This will be a 3d coordinate of shape 30x12x4."
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"import cf_units\n",
"import netcdftime\n",
"\n",
"time_unit = cf_units.Unit('hours since 2000-01-01')\n",
"times = np.empty([30, 12, 4])\n",
"\n",
"for i_y, i_m, i_d in np.ndindex(*times.shape):\n",
" dt = netcdftime.datetime(year.points[i_y], month.points[i_m], day.points[i_d])\n",
" times[i_y, i_m, i_d] = time_unit.date2num(dt)\n",
"\n",
"time_coord = iris.coords.AuxCoord(times, standard_name='time', units=time_unit)\n",
"cube.add_aux_coord(time_coord, data_dims=[0, 1, 2])"
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"air_temperature / (kelvin) (year: 30; month: 12; day: 4; latitude: 500; longitude: 1000)\n",
" Dimension coordinates:\n",
" year x - - - -\n",
" month - x - - -\n",
" day - - x - -\n",
" latitude - - - x -\n",
" longitude - - - - x\n",
" Auxiliary coordinates:\n",
" time x x x - -\n",
" Attributes:\n",
" source: Data from xyz\n"
]
}
],
"source": [
"print(cube)"
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"/home/pelson/dev/iris/lib/iris/fileformats/netcdf.py:2272: IrisDeprecation: NetCDF default saving behaviour currently assigns the outermost dimensions to unlimited. This behaviour is to be deprecated, in favour of no automatic assignment. To switch to the new behaviour, set iris.FUTURE.netcdf_no_unlimited to True.\n",
" warn_deprecated(msg)\n"
]
}
],
"source": [
"iris.save(cube, '/data/pelson/delme.nc')"
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"netcdf delme {\r\n",
"dimensions:\r\n",
"\tyear = UNLIMITED ; // (30 currently)\r\n",
"\tmonth = 12 ;\r\n",
"\tday = 4 ;\r\n",
"\tlatitude = 500 ;\r\n",
"\tlongitude = 1000 ;\r\n",
"variables:\r\n",
"\tdouble air_temperature(year, month, day, latitude, longitude) ;\r\n",
"\t\tair_temperature:standard_name = \"air_temperature\" ;\r\n",
"\t\tair_temperature:units = \"kelvin\" ;\r\n",
"\t\tair_temperature:coordinates = \"time\" ;\r\n",
"\tint64 year(year) ;\r\n",
"\t\tyear:units = \"1\" ;\r\n",
"\t\tyear:long_name = \"year\" ;\r\n",
"\tint64 month(month) ;\r\n",
"\t\tmonth:units = \"1\" ;\r\n",
"\t\tmonth:long_name = \"month\" ;\r\n",
"\tint64 day(day) ;\r\n",
"\t\tday:units = \"1\" ;\r\n",
"\t\tday:long_name = \"day\" ;\r\n",
"\tdouble latitude(latitude) ;\r\n",
"\t\tlatitude:axis = \"Y\" ;\r\n",
"\t\tlatitude:units = \"degrees_north\" ;\r\n",
"\t\tlatitude:standard_name = \"latitude\" ;\r\n",
"\tdouble longitude(longitude) ;\r\n",
"\t\tlongitude:axis = \"X\" ;\r\n",
"\t\tlongitude:units = \"degrees_east\" ;\r\n",
"\t\tlongitude:standard_name = \"longitude\" ;\r\n",
"\tdouble time(year, month, day) ;\r\n",
"\t\ttime:units = \"hours since 2000-01-01\" ;\r\n",
"\t\ttime:standard_name = \"time\" ;\r\n",
"\t\ttime:calendar = \"gregorian\" ;\r\n",
"\r\n",
"// global attributes:\r\n",
"\t\t:_NCProperties = \"version=1|netcdflibversion=4.4.1|hdf5libversion=1.8.17\" ;\r\n",
"\t\t:source = \"Data from xyz\" ;\r\n",
"\t\t:Conventions = \"CF-1.5\" ;\r\n",
"}\r\n"
]
}
],
"source": [
"!ncdump -h /data/pelson/delme.nc"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Single time dimension approach\n",
"\n",
"Instead of having (orthogonal) dimensions for each of the time components, we can combine them into a single dimension, and use the API to pick out the values of interest:"
]
},
{
"cell_type": "code",
"execution_count": 16,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"unknown / (unknown) (-- : 1440; -- : 500; -- : 1000)\n"
]
}
],
"source": [
"data = np.empty(shape=[30 * 12 * 4, 500, 1000], dtype=np.float64)\n",
"\n",
"cube2 = iris.cube.Cube(data)\n",
"print(cube2)"
]
},
{
"cell_type": "code",
"execution_count": 17,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"air_temperature / (kelvin) (-- : 1440; latitude: 500; longitude: 1000)\n",
" Dimension coordinates:\n",
" latitude - x -\n",
" longitude - - x\n",
" Attributes:\n",
" source: Data from xyz\n"
]
}
],
"source": [
"cube2.rename('air_temperature')\n",
"cube2.units = 'kelvin'\n",
"cube2.attributes['source'] = 'Data from xyz'\n",
"\n",
"lat = iris.coords.DimCoord(points=np.linspace(-90, 90, 500),\n",
" standard_name='latitude', units='degrees')\n",
"cube2.add_dim_coord(lat, data_dim=1)\n",
"\n",
"lon = iris.coords.DimCoord(points=np.linspace(-180, 180, 1000),\n",
" standard_name='longitude', units='degrees')\n",
"cube2.add_dim_coord(lon, data_dim=2)\n",
"\n",
"print(cube2)"
]
},
{
"cell_type": "code",
"execution_count": 18,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"air_temperature / (kelvin) (time: 1440; latitude: 500; longitude: 1000)\n",
" Dimension coordinates:\n",
" time x - -\n",
" latitude - x -\n",
" longitude - - x\n",
" Attributes:\n",
" source: Data from xyz\n"
]
}
],
"source": [
"time_unit = cf_units.Unit('days since 2000-01-01')\n",
"times = []\n",
"for year in np.arange(1980, 2010):\n",
" for month in np.arange(1, 13):\n",
" for day in [1, 9, 17, 25]:\n",
" times.append(netcdftime.datetime(year, month, day))\n",
"\n",
"time = iris.coords.DimCoord(points=time_unit.date2num(times),\n",
" standard_name='time', units=time_unit)\n",
"cube2.add_dim_coord(time, data_dim=0)\n",
"print(cube2)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"But how do we access useful things, like the first day of each month? First, let's add some extra coordinate information:"
]
},
{
"cell_type": "code",
"execution_count": 19,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"air_temperature / (kelvin) (time: 1440; latitude: 500; longitude: 1000)\n",
" Dimension coordinates:\n",
" time x - -\n",
" latitude - x -\n",
" longitude - - x\n",
" Auxiliary coordinates:\n",
" day_of_month x - -\n",
" month x - -\n",
" year x - -\n",
" Attributes:\n",
" source: Data from xyz\n"
]
}
],
"source": [
"import iris.coord_categorisation as iccat\n",
"\n",
"iccat.add_day_of_month(cube2, 'time')\n",
"iccat.add_month(cube2, 'time')\n",
"iccat.add_year(cube2, 'time')\n",
"print(cube2)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"What do these coordinates look like, and how can we use them?"
]
},
{
"cell_type": "code",
"execution_count": 20,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"[ 1 9 17 25 1 9 17 25 1 9]\n"
]
}
],
"source": [
"print(cube2.coord('day_of_month').points[0:10])"
]
},
{
"cell_type": "code",
"execution_count": 21,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"air_temperature / (kelvin) (time: 360; latitude: 500; longitude: 1000)\n",
" Dimension coordinates:\n",
" time x - -\n",
" latitude - x -\n",
" longitude - - x\n",
" Auxiliary coordinates:\n",
" day_of_month x - -\n",
" month x - -\n",
" year x - -\n",
" Attributes:\n",
" source: Data from xyz\n"
]
}
],
"source": [
"print(cube2.extract(iris.Constraint(day_of_month=1)))"
]
},
{
"cell_type": "code",
"execution_count": 22,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"['Jan' 'Jan' 'Jan' 'Jan' 'Feb' 'Feb' 'Feb' 'Feb' 'Mar' 'Mar']\n"
]
}
],
"source": [
"print(cube2.coord('month').points[0:10])"
]
},
{
"cell_type": "code",
"execution_count": 23,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"air_temperature / (kelvin) (time: 120; latitude: 500; longitude: 1000)\n",
" Dimension coordinates:\n",
" time x - -\n",
" latitude - x -\n",
" longitude - - x\n",
" Auxiliary coordinates:\n",
" day_of_month x - -\n",
" month x - -\n",
" year x - -\n",
" Attributes:\n",
" source: Data from xyz\n"
]
}
],
"source": [
"print(cube2.extract(iris.Constraint(month='Jan')))"
]
},
{
"cell_type": "code",
"execution_count": 24,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"[1980 1980 1980 1980 1980 1980 1980 1980 1980 1980]\n"
]
}
],
"source": [
"print(cube2.coord('year').points[0:10])"
]
},
{
"cell_type": "code",
"execution_count": 25,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"air_temperature / (kelvin) (time: 48; latitude: 500; longitude: 1000)\n",
" Dimension coordinates:\n",
" time x - -\n",
" latitude - x -\n",
" longitude - - x\n",
" Auxiliary coordinates:\n",
" day_of_month x - -\n",
" month x - -\n",
" year x - -\n",
" Attributes:\n",
" source: Data from xyz\n"
]
}
],
"source": [
"print(cube2.extract(iris.Constraint(year=1980)))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We can combine this together to get the first of the month for all March, April & May in the 1980s..."
]
},
{
"cell_type": "code",
"execution_count": 26,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"air_temperature / (kelvin) (time: 30; latitude: 500; longitude: 1000)\n",
" Dimension coordinates:\n",
" time x - -\n",
" latitude - x -\n",
" longitude - - x\n",
" Auxiliary coordinates:\n",
" day_of_month x - -\n",
" month x - -\n",
" year x - -\n",
" Attributes:\n",
" source: Data from xyz\n",
"DimCoord([1980-03-01 00:00:00, 1980-04-01 00:00:00, 1980-05-01 00:00:00,\n",
" 1981-03-01 00:00:00, 1981-04-01 00:00:00, 1981-05-01 00:00:00,\n",
" 1982-03-01 00:00:00, 1982-04-01 00:00:00, 1982-05-01 00:00:00,\n",
" 1983-03-01 00:00:00, 1983-04-01 00:00:00, 1983-05-01 00:00:00,\n",
" 1984-03-01 00:00:00, 1984-04-01 00:00:00, 1984-05-01 00:00:00,\n",
" 1985-03-01 00:00:00, 1985-04-01 00:00:00, 1985-05-01 00:00:00,\n",
" 1986-03-01 00:00:00, 1986-04-01 00:00:00, 1986-05-01 00:00:00,\n",
" 1987-03-01 00:00:00, 1987-04-01 00:00:00, 1987-05-01 00:00:00,\n",
" 1988-03-01 00:00:00, 1988-04-01 00:00:00, 1988-05-01 00:00:00,\n",
" 1989-03-01 00:00:00, 1989-04-01 00:00:00, 1989-05-01 00:00:00], standard_name='time', calendar='gregorian')\n"
]
}
],
"source": [
"first_mam_80s = iris.Constraint(day_of_month=1, month=['Mar', 'Apr', 'May'],\n",
" year=lambda c: 1980 <= c.point < 1990)\n",
"print(cube2.extract(first_mam_80s))\n",
"print(cube2.extract(first_mam_80s).coord('time'))"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python [conda env:pre-proc]",
"language": "python",
"name": "conda-env-pre-proc-py"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.5.2"
}
},
"nbformat": 4,
"nbformat_minor": 1
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment