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']
var = f.variables['PH2']

print var.shape

b=var[630][:]
print b.shape

var.missing_value=0
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]

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

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.show()