In [179]:
%matplotlib inline
This demo requires that you download, from Brightspace, the following files and place them in your working directory:
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 netcdf4NOw 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:
Click "generate url"...
Below, create a valiable called url and paste the url from ERDDAP
In [180]:
# 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)]'
I'll do the next step...
In [181]:
import urllib
# store data in NetCDF file
fileName='sst.nc'
urllib.urlretrieve (url, fileName)
Out[181]:
Above I:
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"
In [182]:
import netCDF4
# open NetCDF data in
nc = netCDF4.Dataset(file)
ncv = nc.variables
Great!
Now type below ncv to see its contents
In [183]:
ncv
Out[183]:
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
In [184]:
type(ncv)
Out[184]:
What "keys" are contained in the ncv dictionary?
Hint: Use method .keys()
In [185]:
ncv.keys()
Out[185]:
Lets extract the latitues and longitudes into two arrays. I'll do it for you.
In [186]:
lon = ncv['longitude'][:]
lat = ncv['latitude'][:]
Check out the contents of lat
In [187]:
lat
Out[187]:
Now check out the contents of lon
In [188]:
lon
Out[188]:
What type of object are lat and lon?
In [189]:
type(lat)
Out[189]:
In [190]:
type(lon)
Out[190]:
What is their shape?
Hint: use .shape method
In [191]:
lat.shape
Out[191]:
In [192]:
lon.shape
Out[192]:
Step below is tricky to explain. Meshgrid created a two-dimensional matrix (a mesh) of lats and lons
In [193]:
# Make a grid
lons, lats = np.meshgrid(lon,lat)
Take a look inside lons
In [194]:
lons
Out[194]:
What type of object is lons?
In [195]:
type(lons)
Out[195]:
What is the shape of lons?
hint: use .shape method
In [196]:
lons.shape
Out[196]:
What is the shape pf lats?
In [197]:
lats.shape
Out[197]:
Perhaps it is best if we plot lats and lons so you get a better idea of what the are:
In [198]:
%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')
Out[198]:
Next, lets extract Sea Surface Temperature (SST) from our variable ncv
In [199]:
# Get SST
sst = ncv['sst'][0,0,:,:]
Take a look inside sst
In [200]:
sst
Out[200]:
What type of object is sst?
In [201]:
type(sst)
Out[201]:
What is the shape of sst?
In [202]:
sst.shape
Out[202]:
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
In [203]:
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,:,:]
In [204]:
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()
Out[204]:
See how easy is to change projection and make another map
In [205]:
# 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()
Out[205]:
Lets do a map using "kav7" projection
In [206]:
# 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()
Out[206]:
One more using "Lamber Comformal Conic Projection". Here I added some extra labels and colorbar
In [207]:
# 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)')
Out[207]:
This is how you save your figure:
In [208]:
fig1.savefig('myplot.png')
In [209]:
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)')
Out[209]:
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.
In [210]:
# 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])
Take a look inside new_sst
In [211]:
new_sst
Out[211]:
What is the shape of new_sst?
In [212]:
new_sst.shape
Out[212]:
Lets do a quick plot:
In [213]:
import matplotlib.pyplot as plt
fig, (ax) = plt.subplots(1,1)
ax.plot(new_sst,'b-')
Out[213]:
In [214]:
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)
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
In [215]:
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)
What are the differences between the model forced with SST and the one without forcing?
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
Below is how you use the url to get data (using pandas)
In [216]:
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)
What type of object is data?
In [217]:
type(data)
Out[217]:
Take a look into data
In [218]:
data
Out[218]:
What are the "keys" of data?
In [219]:
data.keys()
Out[219]:
Take a look into data['disOxygen_an']
In [220]:
data['disOxygen_an']
Out[220]:
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.
In [221]:
Oxy = data['disOxygen_an'][1:]
Oxy = Oxy.append(data['disOxygen_an'][1:])
Oxy = Oxy.append(data['disOxygen_an'][1:])
Take a look inside Oxy
In [222]:
Oxy
Out[222]:
In [223]:
#Change pandas object to nparray
Oxy = Oxy.values
In [224]:
Oxy
Out[224]:
In [225]:
Oxy = Oxy.astype(float)
In [226]:
Oxy
Out[226]:
In [227]:
# 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])
In [228]:
import matplotlib.pyplot as plt
fig, (ax) = plt.subplots(1,1)
ax.plot(new_Oxy,'b-')
Out[228]:
In [229]:
new_Oxy.shape
new_Oxy
Out[229]:
In [230]:
# Oxy = Oxy.astype(float)
In [231]:
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)
In [232]:
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)
In [ ]: