# Importing libraries we will need.
import netCDF4
%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np

# Open the netCDF file
f = netCDF4.Dataset('data/rtofs_glo_3dz_f006_6hrly_reg3.nc', 'r')

# Getting the n-dimensional data
tempv = f.variables['temperature']
depth = f.variables['Depth']

print("The temperature variable...\n")
# Note the temperature variable has time, depth, y and x dimensions
print(tempv)
print("The dimensions...\n")
print(tempv.dimensions)

# Continued from previous cell..

# Our goal is temperature as a function of depth so slicing along the depth axis
# at a specific time and place
temp = tempv[0, :, 123, 486]

# Masked arrays are arrays that have bad values idenitifed by the mask array.
print("The masked array containing the temperature data...")
print(repr(temp))

# trick for filtering out good values
x = temp[~temp.mask] 
y = depth[~temp.mask]

print("Plotting...")
# plot and show data
plt.plot(y, x)

# close netCDF
f.close()

# Adding text, adjusting borders, and figure size
# Get figure hook to manipulate our plot
fig = plt.figure()
desc ='Figure 1. Temperature as a function of ocean depth as\n'\
    'predicted by the RTOFS model'
# Adding our description
plt.figtext(.5,.15,desc,fontsize=12,ha='center')
#adjust margin
fig.subplots_adjust(bottom=0.35)
#adjust figure size
fig.set_size_inches(7,7)

# Improve axes
# Get axis hook to manipulate our plot
ax = fig.add_subplot(111)
# Add axis labels
ax.set_xlabel('Depth (m)', fontweight='bold')
# \u00b0 : degree symbol
ax.set_ylabel(u'Temperature (\u00b0C)', fontweight='bold')
# Don't show top and right axis
ax.spines["right"].set_visible(False)
ax.spines["top"].set_visible(False)
# Define ticks
ax.tick_params(axis='both', direction='out')
ax.get_xaxis().tick_bottom() # remove unneeded ticks 
ax.get_yaxis().tick_left()

# Getting the data as we did before
f = netCDF4.Dataset('data/rtofs_glo_3dz_f006_6hrly_reg3.nc', 'r')
tempv = f.variables['temperature']
depth = f.variables['Depth']
temp = tempv[0,:,123,486]
x = temp[~temp.mask] #trick for getting data
y = depth[~temp.mask]

# Plotting line with triangle markers, and red line color.
plt.plot(y, x, marker=r'^', color='r', markersize=10, clip_on=False, linewidth=2.0)
plt.show()
f.close()

# Importing CartoPy
import cartopy

# Works with matplotlib's built-in transform support.
fig = plt.figure(figsize=(10, 8))
ax = fig.add_subplot(1, 1, 1, projection=cartopy.crs.Mercator())

# Sets the extent to cover the whole globe
ax.set_global()

# Adds standard background map
ax.stock_img()

fig = plt.figure(figsize=(10, 8))
ax = fig.add_subplot(1, 1, 1, projection=cartopy.crs.PlateCarree())

# Add variety of features
ax.add_feature(cartopy.feature.LAND)
ax.add_feature(cartopy.feature.OCEAN)
ax.add_feature(cartopy.feature.COASTLINE)

# Can also supply matplotlib kwargs
ax.add_feature(cartopy.feature.BORDERS, linestyle=':')
ax.add_feature(cartopy.feature.LAKES, alpha=0.5)
ax.add_feature(cartopy.feature.RIVERS)

# Set limits in lat/lon space
ax.set_extent([-140, -60, 25, 50])

fig = plt.figure(figsize=(10, 8))
ax = fig.add_subplot(1, 1, 1, projection=cartopy.crs.PlateCarree())

ax.add_feature(cartopy.feature.LAND)
ax.add_feature(cartopy.feature.OCEAN)
ax.add_feature(cartopy.feature.COASTLINE)
ax.add_feature(cartopy.feature.BORDERS, linestyle=':')

# Grab state borders
state_borders = cartopy.feature.NaturalEarthFeature(category='cultural', name='admin_1_states_provinces_lines', scale='50m',
                                                    facecolor='none') 

ax.add_feature(state_borders, linestyle="--", edgecolor='blue')

# Set limits in lat/lon space
ax.set_extent([-140, -60, 25, 50])

fig = plt.figure(figsize=(12,8))
ax = fig.add_subplot(1, 1, 1, projection=cartopy.crs.Robinson())
ax.set_global()
ax.stock_img()
ax.coastlines()

# Plot circle with and without specifying the transform
ax.plot(-105.295175, 40.013176, 'o')
ax.plot(-105.295175, 40.013176, 'ro', transform=cartopy.crs.Geodetic())

# Plot line between two lat/lon points. One using equi-rectangular and one geodetic
plt.plot([-0.08, 132.], [51.53, 43.17], 'r', transform=cartopy.crs.PlateCarree())
plt.plot([-0.08, 132.], [51.53, 43.17], 'g', transform=cartopy.crs.Geodetic())

# Open and read netCDF variables
nc = netCDF4.Dataset('data/rtofs_glo_3dz_f006_6hrly_reg3.nc', 'r')
tempv = nc.variables['temperature']
lat = nc.variables['Latitude'][:]
lon = nc.variables['Longitude'][:]
data = tempv[0, 0, :, :]

# Set up a stereographic projection
proj = cartopy.crs.Stereographic(central_latitude=60., central_longitude=-120.)

# Construct figure
fig = plt.figure(figsize=(10,8))
ax = fig.add_subplot(1, 1, 1, projection=proj)

# Setting the plot size and text
plt.figtext(.5, .15, 'Figure 1. Sea surface temperature as predicted by the RTOFS model', fontsize=12, ha='center')

# define color map
cmap = plt.cm.hsv

# Nice high-level, human-readable abstractions for dealing with maps.
ax.add_feature(cartopy.feature.LAND)
ax.add_feature(cartopy.feature.OCEAN)
ax.coastlines(zorder=3)
ax.gridlines()

#Color-filled contour plot
cs = ax.contourf(lon, lat, data, 50, cmap=cmap, transform=cartopy.crs.PlateCarree(), zorder=2)

# Color bar
cbar = plt.colorbar(cs)
cbar.set_label('$^{o}C$')

nc.close()

from datetime import datetime
date = datetime(2014, 9, 15, 0) # date to plot.

# open dataset.
dataset = netCDF4.Dataset('http://www.ncdc.noaa.gov/thredds/dodsC/OISST-V2-AVHRR_agg')
timevar = dataset.variables['time']
timeindex = netCDF4.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 lats and lons (representing centers of grid boxes).
lats = dataset.variables['lat'][:]
lons = dataset.variables['lon'][:]

# create figure, axes instances.
fig = plt.figure(figsize=(20, 10))
proj = cartopy.crs.Robinson(central_longitude=-105.)
ax = fig.add_axes([0.05, 0.05, 0.9, 0.9], projection=proj)

# plot sst, then ice with imshow--this relies on evenly spaced points
bounds = (lons[0], lons[-1], lats[0], lats[-1])
im1 = ax.imshow(sst, extent=bounds, interpolation='nearest',
                cmap=plt.cm.jet, transform=cartopy.crs.PlateCarree())

# Add gridlines and coastlines
ax.gridlines()
ax.coastlines('50m')

# add colorbar
plt.colorbar(im1)

# add a title.
ax.set_title('Analysis for %s' % date)

# specify date to plot.
yyyy=1993; mm=3; dd=14; hh=0

# set OpenDAP server URL.
URLbase = "http://nomads.ncdc.noaa.gov/thredds/dodsC/modeldata/cmd_pgbh/"
URL = URLbase + "{year}/{year}{month:02}/{year}{month:02}{day:02}/pgbh00.gdas.{year}{month:02}{day:02}{hour:02}.grb2".format(
    year=yyyy, month=mm, day=dd, hour=hh)
data = netCDF4.Dataset(URL)

# read lats,lons
latitudes = data.variables['lat'][:]
longitudes = data.variables['lon'][:]

# get sea level pressure and 10-m wind data.
# mult slp by 0.01 to put in units of hPa.
slpin = 0.01 * data.variables['Pressure_msl'][:].squeeze()
uin = data.variables['U-component_of_wind_height_above_ground'][:].squeeze()
vin = data.variables['V-component_of_wind_height_above_ground'][:].squeeze()

# make 2-d grid of lons, lats
lons, lats = np.meshgrid(longitudes, latitudes)

# make orthographic basemap.
proj = cartopy.crs.Orthographic(central_longitude=-60., central_latitude=60.)

# create figure, add axes
fig = plt.figure(figsize=(20, 10))
ax = fig.add_axes([0.1, 0.1, 0.8, 0.8], projection=proj)
ax.set_global()

# set desired contour levels.
clevs = np.arange(960., 1061., 8.)

# plot SLP contours.
ax.contour(lons, lats, slpin, clevs, linewidths=0.5, colors='k', transform=cartopy.crs.PlateCarree())
ax.contourf(lons, lats, slpin, clevs, cmap=plt.cm.RdBu_r, transform=cartopy.crs.PlateCarree())

stride = 15
ax.quiver(lons[::stride, ::stride], lats[::stride, ::stride],
          uin[::stride, ::stride], vin[::stride, ::stride],
          transform=cartopy.crs.PlateCarree())

# draw coastlines and gridlines
ax.coastlines(linewidth=1.5)
ax.gridlines()

date = datetime(yyyy, mm, dd, hh)
ax.set_title('SLP and Wind Vectors ' + str(date))