!conda install -y netcdf4

from netCDF4 import Dataset, num2date, date2num
from numpy import * 
import matplotlib.pyplot as plt 
%matplotlib inline

from ipywidgets import interact, interactive, fixed
import ipywidgets as widgets 

x = linspace(0, 1, 100) # generates a hundred values between 0 and 1
f = 2
a = 3

plt.plot(x, sin(2*pi*x*f))

def pltsin(f):
    plt.plot(x, sin(2*pi*x*f))

pltsin(0.5)

"""
Explanation: Widgets and Interactions
End of explanation
"""


interact(pltsin, f=(1, 10, 0.2), x = (1, 10, 0.2))

def pltsina(f, a):
    plt.plot(x, a*sin(2*pi*x*f))
    plt.ylim(-10.5, 10.5)

interact(pltsina, f=(1, 10, 0.2), a = (1, 10, 0.2))


"""
Explanation: Add to the function to allow amplitude to be varied and aadd in an additional slider to vary both f and a

may want to limit ylim
End of explanation
"""


f=Dataset ('ncep-data/air.sig995.2013.nc') # get individual data set out of the right folder 

air = f.variables['air'] # get variable

plt.imshow(air[0,:,:]) # display first timestep

# Create function to browse through the days 

def sh(time):
    plt.imshow(air[time,:,:])

# Now make it interactive

interact(sh, time=(0, 355, 1))

# Browse variable 

def sh(time =0, var='air',  year = '2013'):
    f=Dataset('ncep-data/'+var+'.sig995.'+year+'.nc')
    vv=f.variables[var]
    plt.imshow(vv[time,:,:])

#Give a list of variables 

variabs =['air', 'uwnd', 'vwnd', 'rhum']
year = ['2013', '2014', '2015']

# Now interact with it 

interact(sh, time=(0, 355, 1), year = year, var=variabs)

help(sh)

from mpl_toolkits.basemap import Basemap

# create north polar sterographic projection 

m=Basemap(projection='npstere', boundinglat=60, lon_0=0, resolution ='l')
m.fillcontinents(color='gray', lake_color='gray')
m.drawparallels(arange(-80.,81.,20.))
m.drawmeridians(arange(-180.,181.,20.))
m.drawmapboundary(fill_color='white')

# Set up some variables 
lon = f.variables['lon'][:]
lat = f.variables['lat'][:]
lon, lat = meshgrid(lon, lat)
x, y = m(lon, lat)

def sh(time =0, var='air',  year = '2013'):
    f=Dataset('ncep-data/'+var+'.sig995.'+year+'.nc')
    vv=f.variables[var]
    tt=f.variables['time']
    dd=num2date(tt[time], tt.units)
    m.fillcontinents(color='gray', lake_color='gray')
    m.drawparallels(arange(-80.,81.,20.))
    m.drawmeridians(arange(-180.,181.,20.))
    m.drawmapboundary(fill_color='white')
    cs = m.contourf(x, y, vv[time,:,:]-273.15)

interact(sh, year=year, time=(0,355,1), var=variabs)

my_map = Basemap (projection='merc', lat_0=0, lon_0=30,
                 resolution='h', area_thresh=1000.0,
                 llcrnrlon=29, llcrnrlat=-1,
                 urcrnrlon=31, urcrnrlat=1)
# area threshold states how rivers etc look - scale, resolution sets resolution, llcrnlon etc sets box, 
# lat and lon decide where you look 
my_map.drawcoastlines()

my_map.drawcountries()
my_map.fillcontinents(color='coral')

my_map.drawmapboundary()

my_map.drawmeridians(arange(0,360,30))
my_map.drawparallels(arange(-90, 90, 30))

lon=30
lat=0

x,y=my_map(lon, lat)
my_map.plot(x, y, 'bo', markersize=7.2)
plt.show() # here the function that decides actually plots 

# This just lets the output of the following code samples 
# display inline on this page, at an appropirate size
from pylab import rcParams

# Create a simple basemap 

my_map = Basemap (projection='ortho', lat_0=50, lon_0=0,
                 resolution='l', area_thresh=1000.0)

my_map.drawcoastlines()

my_map.drawcountries()
my_map.fillcontinents(color='red', lake_color='gray')


plt.show()

"""
Explanation: Climate data
End of explanation
"""


#Check the first few lats and longs

import csv

# Open the earthquake data file.
filename = '1.0_week.csv'

# Create empty lists for the latitudes and longitudes.
lats, lons, mags = [], [], []

# Read through the entire file, skip the first line,
#  and pull out just the lats and lons.
with open(filename) as f:
    # Create a csv reader object.
    reader = csv.reader(f)
    
    # Ignore the header row.
    next(reader)
    
    # Store the latitudes and longitudes in the appropriate lists.
    for row in reader:
        lats.append(float(row[1]))
        lons.append(float(row[2]))
        mags.append(float(row[4]))
        
# Display the first 5 lats and lons.
print('lats', lats[0:5])
print('lons', lons[0:5])
print('mags', mags[0:5])

### And now create a plot of these on a map projection

import csv

# Open the earthquake data file.
filename = '1.0_week.csv'

# Create empty lists for the latitudes and longitudes.
lats, lons, mags = [], [], []

# Read through the entire file, skip the first line,
#  and pull out just the lats and lons.
with open(filename) as f:
    # Create a csv reader object.
    reader = csv.reader(f)
    
    # Ignore the header row.
    next(reader)
    
    # Store the latitudes and longitudes in the appropriate lists.
    for row in reader:
        lats.append(float(row[1]))
        lons.append(float(row[2]))
        mags.append(float(row[4]))
        
# --- Build Map ---
from mpl_toolkits.basemap import Basemap
import matplotlib.pyplot as plt
import numpy as np
 
eq_map = Basemap(projection='robin', resolution = 'l', area_thresh = 1000.0,
              lat_0=52, lon_0=0)
eq_map.drawcoastlines()
eq_map.drawcountries()
eq_map.fillcontinents(color = 'coral')
eq_map.drawmapboundary()
eq_map.drawmeridians(np.arange(0, 360, 30))
eq_map.drawparallels(np.arange(-90, 90, 30))
 
min_marker_size = 1
for lon, lat, mags in zip(lons, lats, mags):
    x,y = eq_map(lon, lat) 
    msize = mags * min_marker_size  
    eq_map.plot(x, y, , markersize=msize)
    if mags >= 5.0 
        eqcolor = 'r'
    elif: mags >= 1.0 and <= 3.0  
        eqcolor = 'g'
    elif: <= 1.0
        eqcolor = 'y
    eq_map.plot(x, y, eqcolor, markersize=msize)
    
plt.show()



"""
Explanation: Plotting some live (ish) earthquake data...
Download the data first: http://earthquake.usgs.gov/earthquakes/feed/v1.0/summary/1.0_week.csv
This will download a file locally- move it into your working directory. Alternatively, use the historic dataset provided in this repo.
End of explanation
"""


x,y

"""
Explanation: This is great but one cool enhancement would be to make the size of the point represent the magnitude of the earthquake.
Here's one way to do it:
Read the magnitudes into a list along with their respective lat and long
Loop through the list, plotting one point at a time
As the magnitudes start at 1.0, you can just use the magnitude directly as the scale factor
To get the marker size, multiply the magnitude by the smallest dot you want on the map.
Add an extra enhancement of colour:
make small earthquakes
See if you can get similar data, perhaps for Whale sightings, and plot those on a map.
You might even have some of your own data to plot..
End of explanation
"""


