#!/usr/bin/env python
# coding: utf-8

# In[1]:


get_ipython().run_line_magic('load_ext', 'load_style')
get_ipython().run_line_magic('load_style', 'talk.css')


# # Iris and Cartopy

# [iris](http://scitools.org.uk/iris/index.html) and [cartopy](http://scitools.org.uk/cartopy/index.html) are developed by the [UK Met. Office](http://scitools.org.uk/).
# 
# + [iris](http://scitools.org.uk/iris/index.html) is a Python package for analysing and visualising meteorological and oceanographic data sets
# + [cartopy](http://scitools.org.uk/cartopy/index.html) is a Python package for advanced map generation with a simple matplotlib interface.
# <br>
# <br>

# In[2]:


import numpy as np
import matplotlib
import matplotlib.pyplot as plt
import datetime as dt


# In[3]:


get_ipython().run_line_magic('matplotlib', 'inline')


# In[4]:


matplotlib.rcParams['figure.figsize'] = (10.0, 8.0)


# In[5]:


import iris


# In[6]:


import cartopy.crs as ccrs


# In[7]:


import iris.quickplot as qplt


# ### loading one day of Tropical Rainfall Measurement Mission (TRMM) rainfall (mm)

# In[8]:


get_ipython().system('ls ../data/')


# In[9]:


date = dt.datetime(2014,07,23)


# In[10]:


fname = '../data/3B42RT_daily.{}.nc'.format(date.strftime("%Y.%m.%d"))


# In[11]:


trmm = iris.load_cube(fname)


# In[14]:


print(trmm)


# In[16]:


lats = trmm.coord('latitude').points
lons = trmm.coord('longitude').points


# In[17]:


lats


# In[18]:


plt.imshow(trmm[0].data)


# ## quick and dirty mapping using the cartopy `qplt` function

# In[19]:


proj = ccrs.PlateCarree(central_longitude=-180.0)


# In[20]:


import cartopy.crs as ccrs
import matplotlib.pyplot as plt

ax = plt.axes(projection=proj)
ax.coastlines()

qplt.contourf(trmm[0],  np.arange(5,50,5), cmap=plt.get_cmap('Blues'), extend='max')

plt.show()


# ## if you want more control, use the `matplotlib` interface

# The line below makes accessible some `features` that you may want to add to your map  
# 
# see [The cartopy Feature interface doc](http://scitools.org.uk/cartopy/docs/latest/matplotlib/feature_interface.html) for more information

# In[21]:


import cartopy.feature as cfeature 


# In[22]:


f = plt.figure()

ax = plt.axes(projection=proj)

im = ax.contourf(lons, lats, trmm[0].data, np.arange(5,50,5),
             transform=ccrs.PlateCarree(), cmap=plt.get_cmap('Blues'), extend='max')

ax.coastlines(linewidth=1.5)
ax.gridlines(crs=proj, draw_labels=True)
cb = plt.colorbar(im, orientation='horizontal', pad=0.05)
cb.set_label('TRMM rainfall for {}: mm/day'.format(date.strftime('%Y/%m/%d')), fontsize=14)

ax.add_feature(cfeature.LAND, alpha=0.3); 


# In[23]:


f = plt.figure()

ax = plt.axes(projection=proj)

im = ax.contourf(lons, lats, trmm[0].data, np.arange(5,50,5),
             transform=ccrs.PlateCarree(), cmap=plt.get_cmap('Blues'), extend='max')

ax.coastlines(lw=0.5)
ax.gridlines(crs=proj, draw_labels=True)
cb = plt.colorbar(im, orientation='horizontal', pad=0.05)
cb.set_label('TRMM rainfall for {}: mm/day'.format(date.strftime('%Y/%m/%d')), fontsize=14)

ax.stock_img()


# In[ ]: