# For ToolsUI screenshot
from IPython.display import Image
Image(filename='tu1.png')

Image(filename='tu2.png')


import netCDF4
f = netCDF4.Dataset('https://motherlode.ucar.edu/repository/opendap/41f2b38a-4e70-4135-8ff8-dbf3d1dcbfc1/entry.das')

grid_type_code = f.variables['grid_type_code'][0]
print 'grid type code is ' + str(grid_type_code)


grid_type = f.variables['grid_type']
print 'The grid type is ' + ''.join(grid_type[0,:])

nx, ny = f.variables['Nx'][0], f.variables['Ny'][0]
print nx, ny

la1, lo1 = f.variables['La1'][0], f.variables['Lo1'][0]
print la1, lo1


latin1, latin2 = f.variables['Latin1'][0], f.variables['Latin2'][0]
print latin1, latin2


lov = f.variables['LoV'][0]
print lov

dx, dy = f.variables['Dx'][0], f.variables['Dy'][0]
print dx, dy

from pyproj import Proj
prj = Proj(proj='lcc', R=6371200, lat_1=latin1, lat_2=latin2, lat_0=la1 , lon_0=lov)

prj(lo1, la1) # note, longitude is first!

fe, fn = prj(lo1, la1) # note, longitude is first!
prj = Proj(proj='lcc', R=6371200, lat_1=latin1, lat_2=latin2, lat_0=la1 , lon_0=lov, x_0=-fe, y_0=-fn)

prj(lo1, la1) # note, longitude is first!

urcrnrlon, urcrnrlat =  prj(nx*dx,ny*dy, inverse=True)

# And for parameter convenience 
llcrnrlon, llcrnrlat = lo1,la1

print llcrnrlat, llcrnrlon, urcrnrlat, urcrnrlon

from mpl_toolkits.basemap import Basemap

fig = plt.figure(figsize=(10,10))
bm = Basemap(llcrnrlon=llcrnrlon,llcrnrlat=llcrnrlat,urcrnrlon=urcrnrlon,
             urcrnrlat=urcrnrlat,
             projection='lcc',lat_1=latin1,lat_2=latin2,lon_0=lov,
             resolution ='l')
bm.drawcoastlines()
bm.drawcountries()
bm.drawmapboundary(fill_color='#99ffff')
bm.fillcontinents(color='#cc9966',lake_color='#99ffff')
parallels = np.arange(0.,81,10.)
bm.drawparallels(parallels,labels=[False,True,True,False])
meridians = np.arange(10.,351.,20.)
bm.drawmeridians(meridians,labels=[True,False,False,True])

secants = [60,30]
bm.drawparallels(secants,labels=[False,False,False,False],
    dashes=[1000, 1],linewidth=3,linestyle='solid',color = 'red')
central_meridian = [lov]
bm.drawmeridians(central_meridian,labels=[False,False,False,False],
    dashes=[1000, 1],linewidth=3,linestyle='solid',color = 'red')

plt.show()

# Get the temperature from the netCDF file
t = f.variables['th']

fig = plt.figure(figsize=(10,10))
bm = Basemap(llcrnrlon=llcrnrlon,llcrnrlat=llcrnrlat,urcrnrlon=urcrnrlon,
             urcrnrlat=urcrnrlat,
             projection='lcc',lat_1=latin1,lat_2=latin2,lon_0=lov,
             resolution ='l')

lons, lats = bm.makegrid(nx, ny) # get lat/lons of ny by nx evenly space grid.
x, y = bm(lons, lats) # compute map proj coordinates.
cs = bm.contourf(x,y,t[0][0], zorder=1)

bm.drawcoastlines()
bm.drawcountries()
bm.drawmapboundary(fill_color='#99ffff')
bm.fillcontinents(color='#cc9966',lake_color='#99ffff', zorder=0)
parallels = np.arange(0.,81,10.)
bm.drawparallels(parallels,labels=[False,True,True,False])
meridians = np.arange(10.,351.,20.)
bm.drawmeridians(meridians,labels=[True,False,False,True])
cbar = bm.colorbar(cs,location='bottom',pad="5%")
cbar.set_label('kelvin')
plt.show()