!/usr/bin/env python

from mpl_toolkits.basemap import Basemap
from netCDF4 import Dataset, date2index
import numpy as np
import matplotlib.pyplot as plt
from datetime import datetime

start_year = 2014
end_year = 2015
save_folder = "/library/Webserver/Documents/FE550/HW4-Map-Murray/data/"
#month = 1
for year in xrange(start_year, end_year):
    for month in xrange(1, 13):

        date = datetime(year,month,1,0) # date to plot.

        # open dataset.
        dataset = \
        Dataset('http://www.ncdc.noaa.gov/thredds/dodsC/OISST-V2-AVHRR_agg')
        timevar = dataset.variables['time']
        timeindex = date2index(date,timevar) # find time index for desired date.
        # read sst.  Will automatically create a masked array using
        # missing_value variable attribute. 'squeeze out' singleton dimensions.
        sst = dataset.variables['sst'][timeindex,:].squeeze()
        # read ice.
        ice = dataset.variables['ice'][timeindex,:].squeeze()
        # read lats and lons (representing centers of grid boxes).
        lats = dataset.variables['lat'][:]
        lons = dataset.variables['lon'][:]
        lons, lats = np.meshgrid(lons,lats)
        # create figure, axes instances.
        fig = plt.figure()
        ax = fig.add_axes([0.05,0.05,0.9,0.9])
        # create Basemap instance.
        # coastlines not used, so resolution set to None to skip
        # continent processing (this speeds things up a bit)
        m = Basemap(projection='kav7',lon_0=0,resolution=None)
        # draw line around map projection limb.
        # color background of map projection region.
        # missing values over land will show up this color.
        m.drawmapboundary(fill_color='0.3')
        # plot sst, then ice with pcolor
        im1 = m.pcolormesh(lons,lats,sst,shading='flat',cmap=plt.cm.jet,latlon=True)
        im2 = m.pcolormesh(lons,lats,ice,shading='flat',cmap=plt.cm.gist_gray,latlon=True)
        # draw parallels and meridians, but don't bother labelling them.
        m.drawparallels(np.arange(-90.,99.,30.))
        m.drawmeridians(np.arange(-180.,180.,60.))
        # add colorbar
        cb = m.colorbar(im1,"bottom", size="5%", pad="2%")
        # add a title.
        ax.set_title('SST and ICE analysis for %s'%date)
        #save this file with label based on year and month
        MAPOUT = (save_folder + str(year) + "-" + str(month).zfill(2) + "-ice.png")
        plt.savefig(MAPOUT)
        #plt.show()