from mpl_toolkits.basemap import Basemap
import sys
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline

plt.ioff() #Turn interactive plotting off

pos=[(0,90),
     (-36*4,25),(-36*3,-25),(-36*2,25),(-36,-25),(0,25),
     (36,-25),(36*2,25),(36*3,-25),(36*4,25),(36*5,-25),
     (36,-90)]

for n in range(0,len(pos)):
    sys.stdout.write("%d "%(n+1))
    plt.figure(figsize=(4,4))
    map = Basemap(width=9.8E6, height=9.8E6, \
                projection='gnom', lon_0=pos[n][0], lat_0=pos[n][1])
    map.drawmapboundary(fill_color='aqua')
    map.drawcoastlines()
    map.fillcontinents(color='orange',lake_color='aqua')
    map.drawparallels(np.arange(-90,90,10))
    map.drawmeridians(np.arange(-180,180,10))
    plt.axis('off')
    plt.savefig('map_' + "%02d"%(n+1) + '.png', dpi=100,  bbox_inches='tight', pad_inches=0)
    plt.clf()      # clear memmory
    plt.close()
    

from IPython.display import Image, HTML, display
from glob import glob
imagesList=''.join( ["<img style='width: 120px; margin: 0px; float: left; border: 1px solid black;' src='%s' />" % str(s) 
                     for s in sorted(glob('map_*.png')) ])
display(HTML(imagesList))

# Image size is 310x310 when figsize is set to 4,4 and axis not drawn
s0=155
c0=155

c1 = np.cos(2*np.pi/5) * s0
c2 = np.cos(np.pi/5) * s0 + 4   # +4 to get small ovelays
s1 = np.sin(2*np.pi/5) * s0
s2 = np.sin(4*np.pi/5) * s0

import PIL
import PIL.ImageDraw

def mask_pentagon(filein, fileout, up=True):

    # read image as RGB and add alpha (transparency)
    im = PIL.Image.open(filein).convert("RGBA")

    # convert to numpy (for convenience)
    imArray = np.asarray(im)

    # create mask from polygon ABCDE A(0,1), B(s1,c1), C(s2,-c2), D(-s2,-c2), E(-s1,c1)
    # image origin (0,0) is located at upper left corner with positive coordinates downward
    if up:
        A=(s0, 0)
        B=(s0+s1, c0-c1)
        C=(s0+s2, c0+c2)
        D=(s0-s2, c0+c2)
        E=(s0-s1, c0-c1)
    else:
        A=(s0-s2, c0-c2)
        B=(s0+s2, c0-c2)
        C=(s0+s1, c0+c1)
        D=(s0, c0+c0)
        E=(s0-s1, c0+c1) 
    polygon = [A,B,C,D,E]  
    maskIm = PIL.Image.new('L', (imArray.shape[1], imArray.shape[0]), 0)
    PIL.ImageDraw.Draw(maskIm).polygon(polygon, outline=1, fill=1)
    mask = np.array(maskIm)

    # assemble new image (uint8: 0-255)
    newImArray = np.empty(imArray.shape,dtype='uint8')

    # colors (three first columns, RGB)
    newImArray[:,:,:3] = imArray[:,:,:3]

    # transparency (4th column)
    newImArray[:,:,3] = mask*255

    # back to Image from numpy
    newIm = PIL.Image.fromarray(newImArray, "RGBA")
    newIm.save(fileout)

for i in range(1,13) :
    mask_pentagon("map_%02d"%i + ".png", "pentagon_map_%02d"%i + ".png", up=i%2)

imagesList=''.join( ["<img style='width: 120px; margin: 0px; float: left; border: 1px solid black; background: #CCC' src='%s' />" % str(s) 
                     for s in sorted(glob('pentagon_map_*.png')) ])
display(HTML(imagesList))

import netCDF4

f=netCDF4.Dataset("https://prodn.idris.fr/thredds/dodsC/STORE/rfry938/ORCA025-PIS2DIC/MBG/Analyse/TS_MO/surf/ORCA025_1958_2010_1M_PH2_surf.nc")
print f.variables
lats = f.variables['nav_lat']
lons = f.variables['nav_lon']
times = f.variables['time_counter']
var = f.variables['PH2']

print var.shape
var.missing_value=0   # correct missing value

print var.long_name

times = f.variables['time_counter']
print times.shape, times.units, times.calendar

from netCDF4 import num2date, date2num, date2index
dates = num2date(times[:], times.units, calendar=times.calendar)

datesStr = [date.strftime('%Y/%m/%d') for date in dates]
print times[-1], datesStr[-1]

plt.figure(figsize=(10,10))

map = Basemap(projection='ortho', lat_0=-20, lon_0=-120, resolution='c')
map.drawcoastlines(linewidth=0.25)
map.drawcountries(linewidth=0.25)
map.drawmeridians(np.arange(0, 360, 30))
map.drawparallels(np.arange(-90, 90, 30))

x, y = map(lons[:], lats[:])

# colormap: http://matplotlib.org/users/colormaps.html
# extend:  "neither", "both", "min", "max"
timeIndex = -1
map.contourf(x, y, 1E9*var[timeIndex][0], levels=np.arange(7.5,10.5,.2), extend='both', cmap=plt.cm.CMRmap_r) 
plt.colorbar()
plt.title("Surface hydrogen ion concentration [H+] (nmol/kg) - " + datesStr[timeIndex])

plt.savefig('ORCA0.25_pH.png', dpi=100,  bbox_inches='tight', pad_inches=0)

plt.ioff() #Turn interactive plotting off

pos=[(0,90),(-36*4,25),(-36*3,-25),(-36*2,25),(-36,-25),(0,25),
     (36,-25),(36*2,25),(36*3,-25),(36*4,25),(36*5,-25),(36,-90)]

for n in range(0,len(pos)):
    sys.stdout.write("%d "%(n+1))
    plt.figure(figsize=(4,4))
    map = Basemap(width=9.8E6, height=9.8E6, \
                projection='gnom', lon_0=pos[n][0], lat_0=pos[n][1])
    map.drawmapboundary(fill_color='aqua')
    map.drawcoastlines()
    map.drawparallels(np.arange(-90,90,10))
    map.drawmeridians(np.arange(-180,180,10))
    x, y = map(lons[:], lats[:])
    map.contourf(x, y, 1E9*var[timeIndex][0], levels=np.arange(7.5,10.5,.2), extend='both', cmap=plt.cm.CMRmap_r)
    plt.axis('off')
    plt.savefig('mapvar_' + "%02d"%(n+1) + '.png', dpi=100,  bbox_inches='tight', pad_inches=0)
    plt.clf()    # clear memory
    plt.close()

for i in range(1,13) :
    mask_pentagon("mapvar_%02d"%i + ".png", "pentagon_mapvar_%02d"%i + ".png", up=i%2)
    
imagesList=''.join( ["<img style='width: 120px; margin: 0px; float: left; border: 1px solid black; background: #CCC' src='%s' />" % str(s) 
                     for s in sorted(glob('pentagon_mapvar_*.png')) ])
display(HTML(imagesList))