# Set up url
url='https://coastwatch.pfeg.noaa.gov/erddap/griddap/erdAGsstamday_LonPM180.nc?sst[(2015-08-16T12:00:00Z):1:(2015-08-16T12:00:00Z)][(0.0):1:(0.0)][(38):1:(48)][(-68):1:(-50)]'

"""
Explanation: TUTORIAL - Dowload Satellite-Temperature, and Oxygen Climatology... and use it to force MUSEEL model
This demo requires that you download, from Brightspace, the following files and place them in your working directory:

model_Mussel_IbarraEtal2014.py


IMPORTANT!
You will need to import 2 new modules: basemap to make maps... and netCDF4 to hangle .nc files. To do this, open a terminal by going to "anaconda navigator" > click on "Environments" > in the "root" area, click on the "play" symbol > click on "Open Terminal"
In the terminal, write:
conda install basemap netcdf4

NOw open a browser and open google...
First Google "ERDDAP NOAA" and go to the 1st link. Then, on right top corner, click on "View a List of All #### Datasets"
...or simply go here: https://coastwatch.pfeg.noaa.gov/erddap/info/index.html?page=1&itemsPerPage=1000

Lets start with Satellite SST (Sea Surface Temperature)
Find (Control+F): "SST, POES AVHRR, GAC, Global, Day and Night (Monthly Composite), Lon+/-180"
Click on "data"
...or just go here: https://coastwatch.pfeg.noaa.gov/erddap/griddap/erdAGsstamday_LonPM180.html

Now, start by selecting:

Time from 2015-08-16 
Latitude from 38 to 48
Longiture from -67 to -50
File type: .nc (NetCDF)

Click "generate url"... 
Below, create a valiable called url and paste the url from ERDDAP
End of explanation
"""


import urllib

# store data in NetCDF file
fileName='sst.nc'
urllib.urlretrieve (url, fileName)

"""
Explanation: I'll do the next step...
End of explanation
"""


import netCDF4

# open NetCDF data in 
nc = netCDF4.Dataset(file)
ncv = nc.variables

"""
Explanation: Above I:

Imported a library to download data from urls
Created a FileName
Downloaded contents of "url" and save it in "FileName"

Once you run the code above... check to see if the file "sst.nc" was created.

Next, lets open the file we created (sst.nc) and load the "variables"
End of explanation
"""


ncv

"""
Explanation: Great!
Now type below ncv to see its contents
End of explanation
"""


type(ncv)

"""
Explanation: Lots of information here. But for now, lets gloss over and do some tests. Answer the questions below:
What type of object is ncv?
Hint: use the type function to answer this question
End of explanation
"""


ncv.keys()


"""
Explanation: What "keys" are contained in the ncv dictionary?
Hint: Use method .keys()
End of explanation
"""


lon = ncv['longitude'][:]
lat = ncv['latitude'][:]

"""
Explanation: Lets extract the latitues and longitudes into two arrays. I'll do it for you.
End of explanation
"""


lat

"""
Explanation: Check out the contents of lat
End of explanation
"""


lon

"""
Explanation: Now check out the contents of lon
End of explanation
"""


type(lat)

type(lon)

"""
Explanation: What type of object are lat and lon?
End of explanation
"""


lat.shape

lon.shape

"""
Explanation: What is their shape?
Hint: use .shape method
End of explanation
"""


# Make a grid
lons, lats = np.meshgrid(lon,lat)

"""
Explanation: Step below is tricky to explain. Meshgrid created a two-dimensional matrix (a mesh) of lats and lons
End of explanation
"""


lons

"""
Explanation: Take a look inside lons
End of explanation
"""


type(lons)

"""
Explanation: What type of object is lons?
End of explanation
"""


lons.shape

"""
Explanation: What is the shape of lons?
hint: use .shape method
End of explanation
"""


lats.shape

"""
Explanation: What is the shape pf lats?
End of explanation
"""


%matplotlib inline

import matplotlib.pyplot as plt
fig, (ax1, ax2) = plt.subplots(1,2)
ax1.pcolor(lats)
ax1.set_title('lats')
ax2.pcolor(lons)
ax2.set_title('lons')

"""
Explanation: Perhaps it is best if we plot lats and lons so you get a better idea of what the are:
End of explanation
"""


# Get SST
sst = ncv['sst'][0,0,:,:]

"""
Explanation: Next, lets extract Sea Surface Temperature (SST) from our variable ncv
End of explanation
"""


sst

"""
Explanation: Take a look inside sst
End of explanation
"""


type(sst)

"""
Explanation: What type of object is sst?
End of explanation
"""


sst.shape

"""
Explanation: What is the shape of sst?
End of explanation
"""


year = 2015
month = 8
minlat = 38
maxlat = 48
minlon = -67
maxlon = -50
isub = 0.5

minday = str(year)+'-'+str(month).zfill(2)+'-11T12:00:00Z'
maxday = str(year)+'-'+str(month+1).zfill(2)+'-11T12:00:00Z'

base_url='http://coastwatch.pfeg.noaa.gov/erddap/griddap/erdAGsstamday_LonPM180.nc?'
query='sst[('+minday+'):'+str(isub)+':('+maxday+')][(0.0):'+str(isub)+':(0.0)][('+str(minlat)+'):'+str(isub)+':('+str(maxlat)+')][('+str(minlon)+'):'+str(isub)+':('+str(maxlon)+')]'
url = base_url+query

# store data in NetCDF file
file='satellite_data_tempfile.nc'
urllib.urlretrieve (url, file)

# open NetCDF data in 
nc = netCDF4.Dataset(file)
ncv = nc.variables

lon = ncv['longitude'][:]
lat = ncv['latitude'][:]
lons, lats = np.meshgrid(lon,lat)
sst = ncv['sst'][0,0,:,:]

"""
Explanation: Let me put it all together, so that you can see how it looks when all is in one place... run the code below to import sst
End of explanation
"""


import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemap

# Create map (Miller Projection)
m = Basemap(projection='mill', llcrnrlat=minlat,urcrnrlat=maxlat,llcrnrlon=minlon, urcrnrlon=maxlon,resolution='i')
fig1 = plt.figure(figsize=(13,13))
cs = m.pcolormesh(lons,lats,sst,cmap=plt.cm.jet,latlon=True)
m.drawcoastlines()
m.drawmapboundary()

"""
Explanation: Mapping
Lets make a map!!!!
End of explanation
"""


# Create map (Othogonal Projection)
m = Basemap(projection='ortho', lat_0=(maxlat+minlat)/2,lon_0=(maxlon+minlon)/2,resolution='l')
fig1 = plt.figure(figsize=(13,13))
cs = m.pcolormesh(lons,lats,sst,cmap=plt.cm.jet,latlon=True)
m.drawcoastlines()
m.drawmapboundary()

"""
Explanation: See how easy is to change projection and make another map
End of explanation
"""


# Create map (kav7 Projection)
m = Basemap(projection='kav7', lon_0=(maxlon+minlon)/2,resolution='l')
fig1 = plt.figure(figsize=(13,13))
cs = m.pcolormesh(lons,lats,sst,cmap=plt.cm.jet,latlon=True)
m.drawcoastlines()
m.drawmapboundary()

"""
Explanation: Lets do a map using "kav7" projection
End of explanation
"""


# Create map (Lamber Comformal Conic Projection)
m = Basemap(llcrnrlon=-80.5,llcrnrlat=40.,urcrnrlon=-40.566,urcrnrlat=44.352,\
            rsphere=(6378137.00,6356752.3142),\
            resolution='i',area_thresh=1000.,projection='lcc',\
            lat_1=45.,lon_0=-107.)
fig1 = plt.figure(figsize=(13,13))
cs = m.pcolormesh(lons,lats,sst,cmap=plt.cm.jet,latlon=True)
m.drawcoastlines()
m.drawmapboundary()
plt.title('Satellite SST (Monthly composite for '+str(year)+'/'+str(month)+')')
cbar = plt.colorbar(orientation='horizontal', extend='both')
cbar.ax.set_xlabel('Sea Surface Temperature (oC)')

"""
Explanation: One more using "Lamber Comformal Conic Projection". Here I added some extra labels and colorbar
End of explanation
"""


fig1.savefig('myplot.png')

"""
Explanation: This is how you save your figure:
End of explanation
"""


import pandas as pd

minday='2010-01-01T12:00:00Z'
maxday='2012-12-31T12:00:00Z'
lat = 43
lon = -62
isub = 1

base_url='http://coastwatch.pfeg.noaa.gov/erddap/griddap/erdAGsstamday_LonPM180.csv?'
query='sst[('+minday+'):'+str(isub)+':('+maxday+')][(0.0):'+str(isub)+':(0.0)][('+str(lat)+'):'+str(isub)+':('+str(lat)+')][('+str(lon)+'):'+str(isub)+':('+str(lon)+')]'
url = base_url+query
print "Downloading Satellite Data from: POES, AVHRR and GAC (Time-series)..."

data = pd.read_csv(url)
sst = data['sst'][1:].astype(float)
ax = sst.plot()
ax.set_xlabel('months')
ax.set_ylabel('Sea Surface Temperature (oC)')

"""
Explanation: Time-series of SST
Rather than a map (snapshopt in time), lets import a time-series of a single spot
End of explanation
"""


# Resize array and Eliminate Gaps by Interpolating in between
import scipy.interpolate as intrp
import numpy as np

months = len(sst)

x = np.arange(0, months)

f = intrp.interp1d(x, sst, kind='linear', fill_value='extrapolate' )

days = 365 * 3
dt = 0.01
NoSTEPS = int(days / dt)
newx = np.linspace(0,days/30,NoSTEPS) # Makes and vector array of equally spaced numbers from zero to "days"

new_sst = f(newx)

# Get rid of nans
nans, x= np.isnan(new_sst), lambda z: z.nonzero()[0]
new_sst[nans]= np.interp(x(nans), x(~nans), new_sst[~nans])

"""
Explanation: The next part is tricky to explain... it re-shapes the original time-series to make it of the same size of the model output (you need to state "days" and "dt")... and it also "patches" any missing data. It is all done using interpolation and extrapolation.
End of explanation
"""


new_sst

"""
Explanation: Take a look inside new_sst
End of explanation
"""


new_sst.shape

"""
Explanation: What is the shape of new_sst?
End of explanation
"""


import matplotlib.pyplot as plt
fig, (ax) = plt.subplots(1,1)
ax.plot(new_sst,'b-')

"""
Explanation: Lets do a quick plot:
End of explanation
"""


import model_Mussel_IbarraEtal2014 as MusselModel

days, dt, par, InitCond = MusselModel.load_defaults()
output = MusselModel.run_model(days,dt,InitCond,par)
MusselModel.plot_model(output)

"""
Explanation: Focing MUSSEL model with SST
Make sure model_Mussel_IbarraEtal2014 is in your working directory.
End of explanation
"""


import model_Mussel_IbarraEtal2014_SST as MusselModel
reload(MusselModel)

days, dt, par, InitCond = MusselModel.load_defaults()
output = MusselModel.run_model(days,dt,InitCond,par,new_sst)
MusselModel.plot_model(output)

"""
Explanation: Now... open model_Mussel_IbarraEtal2014.py in spyder and save it with a different name (e.g. model_Mussel_IbarraEtal2014_SST.py
Then do the next 3 minor changes:
(1) Indicate sst as input
Line 65 BEFORE: def run_model(days,dt,InitCond,par):
Line 65 AFTER: def run_model(days,dt,InitCond,par,sst):
(2) Initialize Temp with sst
Line 107 BEFORE: Temp = InitCond['Temp']
Line 107 AFTER: Temp = sst
(3) Make Temp as function of time: Temp[t]
Line 122 and 122 BEFORE:
L_Temp[t] = min(max(0.,1.-np.exp(-par['KTempL']*(Temp-par['TempL']))), \
                     max(0.,1.+((1.-np.exp(par['KTempH']*Temp))/(np.exp(par['KTempH']*par['TempH'])-1.))))
Line 122 and 122 AFTER:
```
        L_Temp[t] = min(max(0.,1.-np.exp(-par['KTempL'](Temp[t]-par['TempL']))), \
                     max(0.,1.+((1.-np.exp(par['KTempH']Temp[t]))/(np.exp(par['KTempH']*par['TempH'])-1.))))
Lets run model_Mussel_IbarraEtal2014_SST using our downloaded satellite sst
End of explanation
"""


import pandas as pd

url = 'https://coastwatch.pfeg.noaa.gov/erddap/griddap/nodcWoa09mon1t_LonPM180.csv?temperature_an[(0000-01-16):1:(0000-12-16T00:00:00Z)][(0.0):1:(0)][(43.5):1:(43.5)][(-61.5):1:(-61.5)],salinity_an[(0000-01-16):1:(0000-12-16T00:00:00Z)][(0.0):1:(0)][(43.5):1:(43.5)][(-61.5):1:(-61.5)],disOxygen_an[(0000-01-16):1:(0000-12-16T00:00:00Z)][(0.0):1:(0)][(43.5):1:(43.5)][(-61.5):1:(-61.5)]'

data = pd.read_csv(url)

"""
Explanation: What are the differences between the model forced with SST and the one without forcing?
Forcing with Oxygen Climatology
Here we need to go to the ERDDAP to download climatologies from the World Ocean Atlas 2009: https://coastwatch.pfeg.noaa.gov/erddap/griddap/nodcWoa09mon1t_LonPM180.html
Here we'll download all data for a single spot (like sst above)
Using the sliders, select all dates for latitude = 43.5 and longitude = -61.5

Select only:  temperature_an, salinity_an and disOxygen_an
File type: .csv
Click on "Generate URL"

Below is how you use the url to get data (using pandas)
End of explanation
"""


type(data)

"""
Explanation: What type of object is data?
End of explanation
"""


data


"""
Explanation: Take a look into data
End of explanation
"""


data.keys()

"""
Explanation: What are the "keys" of data?
End of explanation
"""


data['disOxygen_an']

"""
Explanation: Take a look into data['disOxygen_an']
End of explanation
"""


Oxy = data['disOxygen_an'][1:]
Oxy = Oxy.append(data['disOxygen_an'][1:])
Oxy = Oxy.append(data['disOxygen_an'][1:])

"""
Explanation: We only have 1 year of data, but we need a 3-year long time-series. So we'll copy-paste the same data 3 times to create a what we need.
End of explanation
"""


Oxy

#Change pandas object to nparray
Oxy = Oxy.values

Oxy

Oxy = Oxy.astype(float)

Oxy

# Resize array and Eliminate Gaps by Interpolating in between
import scipy.interpolate as intrp
import numpy as np

months = len(Oxy)

x = np.arange(0, months)

print len(x)

f = intrp.interp1d(x, Oxy, kind='linear', fill_value='extrapolate' )

days = 365 * 3
dt = 0.01
NoSTEPS = int(days / dt)
newx = np.linspace(0,days/30,NoSTEPS) # Makes and vector array of equally spaced numbers from zero to "days"

new_Oxy = f(newx)

# Get rid of nans
nans, x= np.isnan(new_Oxy), lambda z: z.nonzero()[0]
new_Oxy[nans]= np.interp(x(nans), x(~nans), new_Oxy[~nans])

import matplotlib.pyplot as plt
fig, (ax) = plt.subplots(1,1)
ax.plot(new_Oxy,'b-')


new_Oxy.shape
new_Oxy

# Oxy = Oxy.astype(float)

import model_Mussel_IbarraEtal2014_SST_Oxy as MusselModel
reload(MusselModel)

days, dt, par, InitCond = MusselModel.load_defaults()
output = MusselModel.run_model(days,dt,InitCond,par,new_sst,new_Oxy)
MusselModel.plot_model(output)


days, dt, par, InitCond = MusselModel.load_defaults()
par['OxyL']  = 2
par['KOxyL'] = 0.04
output = MusselModel.run_model(days,dt,InitCond,par,new_sst,new_Oxy)
MusselModel.plot_model(output)

"""
Explanation: Take a look inside Oxy
End of explanation
"""


