Skip to content

Instantly share code, notes, and snippets.

@dfulu

dfulu/xarray_example.ipynb

Last active October 18, 2019 15:12
Show Gist options

  • Select an option

    Learn more about clone URLs
  • Save dfulu/e9f8eb753f5601532a6eeb9dcc5c5608 to your computer and use it in GitHub Desktop.

Select an option

Learn more about clone URLs
Save dfulu/e9f8eb753f5601532a6eeb9dcc5c5608 to your computer and use it in GitHub Desktop.
Download ZIP
xarray_example
Display the source blob
Display the rendered blob
Raw
xarray_example.ipynb
{
"cells": [
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"!ls /exports/csce/datastore/geos/users/s1205782/climate_data/"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"import numpy as np\n",
"import xarray as xr\n",
"import matplotlib.pyplot as plt\n",
"import cartopy.crs as ccrs"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Load in some data from netcdf"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# pick file\n",
"data_filename = \"/exports/csce/datastore/geos/users/s1205782/climate_data/HadCM3_initial/\" \\\n",
" +\"temp236_apm_Euroclim500_HadCM3_DRIFT_r1_080012_200011.nc\""
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# load data from netcdf\n",
"ds = xr.open_dataset(data_filename)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# quick view on the data\n",
"print(ds)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# rename some of the coords that end with 0\n",
"ds = ds.rename({'time0':'time',\n",
" 'latitude0':'latitude',\n",
" 'longitude0':'longitude'})\n",
"\n",
"print(ds)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# The dimension z10_height is redundant. Get rid of it\n",
"# No real need to do this, just tidying up really\n",
"ds = ds.squeeze('z10_height').drop('z10_height')"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# have a look at the time coordinates\n",
"print(ds.time)\n",
"# You could also select the time coordinates like this too \n",
"#print(ds['time'])"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Make some field plots of the data"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Select temperate data by its time index\n",
"# here we select the 0th and the 6th fields to get fields 6 months apart\n",
"tempfields = ds.temp236.isel(time=[0,6])\n",
"\n",
"\n",
"# temps are in Kelvin, convert to celcius instead\n",
"tempfields = tempfields - 273.15\n",
"\n",
"# Now we will make a plot of these fields\n",
"\n",
"# This is the map projection we want to plot *onto*\n",
"map_proj = ccrs.PlateCarree()\n",
"\n",
"p = tempfields.plot(\n",
" col='time', # make a new plot for each time coordinate\n",
" col_wrap=1, # stack the subplots above each other\n",
" subplot_kws={'projection': map_proj}) # the plot's projection\n",
"\n",
"# We have to set the map's options on all four axes\n",
"# loop over the axes\n",
"for ax in p.axes.flat:\n",
" # add coastlines\n",
" ax.coastlines()\n",
"\n",
"# show\n",
"plt.show()"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# this time we are just going to make a signle plot\n",
"\n",
"# select data from Decmber of year 1676\n",
"tempfield = ds.temp236.sel(time='1676-12') - 273.15\n",
"\n",
"# mapo projection\n",
"map_proj = ccrs.PlateCarree()\n",
"\n",
"# plot figure in a similar way as before\n",
"plt.figure(figsize = (14,6))\n",
"ax = plt.axes(projection=map_proj);\n",
"p = tempfield.plot(ax = ax, transform=map_proj)\n",
"ax.coastlines()\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## calculate area weighted global warming signal"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# find the cos latitudes. We weight each pixel by this\n",
"coslat = np.cos(np.deg2rad(ds.latitude)).where(~ds.temp236.isnull())\n",
"\n",
"# calculate weights to apply, making sure to normalise\n",
"weights = coslat / coslat.sum(dim=('latitude', 'longitude'))\n",
"\n",
"# multiply the temp fields by the weights at each time step and find mean temp\n",
"global_mean_st_weighted = (ds.temp236 * weights).sum(dim=('latitude', 'longitude'))\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# depending on the packages installed you might need to run this line\n",
"# global_mean_st_weighted['time'] = np.arange(len(global_mean_st_weighted.time))"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"plt.figure()\n",
"global_mean_st_weighted.plot(label='mean ST')\n",
"plt.legend()\n",
"plt.show()"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# next take a rolling mean of the temp\n",
"years_smoothing = 10\n",
"rolling_global_mean_st = global_mean_st_weighted.rolling(time=12*years_smoothing, center=True).mean()"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# make figure to plot average surface temp\n",
"plt.figure()\n",
"rolling_global_mean_st.plot(label='rolling {}yr mean ST'.format(years_smoothing))\n",
"plt.legend()\n",
"plt.show()"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# select just some of an area\n",
"tempfield = ds.temp236.sel(time='1676-12', latitude=slice(90,0), longitude=slice(0, 60)) - 273.15\n",
"\n",
"# mapo projection\n",
"map_proj = ccrs.PlateCarree()\n",
"\n",
"# plot figure in a similar way as before\n",
"plt.figure(figsize = (14,6))\n",
"ax = plt.axes(projection=map_proj);\n",
"p = tempfield.plot(ax = ax, transform=map_proj)\n",
"ax.coastlines()\n",
"plt.show()"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"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.7.3"
}
},
"nbformat": 4,
"nbformat_minor": 4
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment