import matplotlib.pyplot as plt
import cartopy.crs as ccrs
from mpl_toolkits.basemap import Basemap
%matplotlib inline
import matplotlib.colors as colors
from numpy import array
from numpy import max
import numpy as np

plt.figure()
filename = "data/Justin/icepmag_20190118.csv"
frame = pd.read_csv(filename)

frame.head()
# replace -9999 by NaN (missing values)
frame.replace(-9999.0,np.NaN, inplace=True)
# remove rows where lon or lat are missing
frame.dropna(subset=['lon', 'lat'], inplace=True)
# drop duplicates
latlons = frame[['lat', 'lon']].drop_duplicates()

map = Basemap(llcrnrlon=335,llcrnrlat=63.0,urcrnrlon=348.,urcrnrlat=67.,
             resolution='i', projection='tmerc', lat_0 = 64.5, lon_0 = 340)


x,y = map(latlons['lon'].values,latlons['lat'].values)
map.hexbin(x,y, gridsize=20, mincnt=1, cmap='summer', norm=colors.LogNorm())


map.drawcoastlines()
plt.show()

import matplotlib
import matplotlib.pyplot as plt
from matplotlib import colors as c
%matplotlib inline
from mpl_toolkits.basemap import Basemap, shiftgrid
import numpy as np
import netCDF4
from ipywidgets import interact

mpl.rcParams['figure.figsize'] = [20., 10.]
def drawmap(data, title):
    ax = plt.axes(projection=ccrs.PlateCarree(central_longitude=180))
    data.z.plot.pcolormesh(
                   transform=ccrs.PlateCarree())
    ax.set_title(data.title)
    ax.coastlines()
    ax.gridlines()
    
def myanimate(i, filenames):
    dset = xr.open_dataset(filenames[i])
    new_contour = drawmap(dset, 'data %03d'%(i) ) 
    dset.close()
   

@interact(depth=(0,20))
def finteract(depth):
     ca = myanimate(depth, filenames)

import matplotlib.animation as animation

def drawmap(ax, data, title):
    data.z.plot.pcolormesh(add_colorbar=False,
                   transform=ccrs.PlateCarree())
    ax.set_title(data.title)
    ax.coastlines()
    ax.gridlines()
    
def myanimate(i, ax, filenames):
    ax.clear()
    dset = xr.open_dataset(filenames[i])
    new_contour = drawmap(ax, dset, 'data %03d'%(i) ) 
    dset.close()

FFMpegWriter = animation.writers['ffmpeg']
metadata = dict(title='Grace data', artist='CEED workshop',
                comment='Movie for sequence of images')
writer = FFMpegWriter(fps=20, metadata=metadata)

# Create the figure

fig = plt.figure(figsize=[20,10])  # a new figure window
ax = fig.add_subplot(1, 1, 1)  # specify (nrows, ncols, axnum)
ax = plt.axes(projection=ccrs.PlateCarree(central_longitude=180))
    
ani = animation.FuncAnimation(fig, myanimate, frames=np.arange(21), 
    fargs=(ax, filenames), interval=100)
ani.save("grace_movie.mp4"))

from IPython.display import HTML
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import numpy as np
mpl.rcParams["animation.embed_limit"]=100

import matplotlib.animation as animation

import glob
def getint(name):
    _, num = name.split('P_')
    num, _ = num.split('.')
    return int(num)

filenames = glob.glob("data/Grace/grd_files_to_interp/GypsumP_*.grd")

filenames = sorted(filenames, key=getint)

def drawmap(ax, data, title):
    data.z.plot.pcolormesh(add_colorbar=False,
                   transform=ccrs.PlateCarree())
    ax.set_title(data.title)
    ax.coastlines()
    ax.gridlines()
    
def myanimate(i, ax, filenames):
    ax.clear()
    dset = xr.open_dataset(filenames[i])
    new_contour = drawmap(ax, dset, 'data %03d'%(i) ) 
    dset.close()

FFMpegWriter = animation.writers['ffmpeg']
metadata = dict(title='Grace data', artist='CEED workshop',
                comment='Movie for sequence of images')
writer = FFMpegWriter(fps=20, metadata=metadata)

# Create the figure
fig = plt.figure(figsize=[20,10])  # a new figure window
ax = fig.add_subplot(1, 1, 1)  # specify (nrows, ncols, axnum)
ax = plt.axes(projection=ccrs.PlateCarree(central_longitude=180))

ani = animation.FuncAnimation(fig, myanimate, frames=np.arange(21), 
    fargs=(ax, filenames), interval=100)

HTML(ani.to_html5_video())

HTML(ani.to_jshtml())

import folium

map = folium.Map(location=[65., 335.0], zoom_start=5 , tiles='Stamen Terrain')
map

map = folium.Map(location=[65., 335.0], zoom_start=5 , 
                 tiles='https://server.arcgisonline.com/ArcGIS/rest/services/World_Imagery/MapServer/tile/{z}/{y}/{x}.png',
                 attr = "IGN")
map

import pandas as pd
import numpy as np

filename = "/opt/uio/deep_python/data/Justin/icepmag_20190118.csv"
frame = pd.read_csv(filename)

frame.head()

# replace -9999 by NaN (missing values)
frame.replace(-9999.0,np.NaN, inplace=True)
# remove rows where lon or lat are missing
frame.dropna(subset=['lon', 'lat'], inplace=True)
# drop duplicates
latlons = frame[['lat', 'lon']].drop_duplicates()

latlons.head()

map = folium.Map(location=[65., 335.0], zoom_start=5 , tiles='Stamen Terrain')
# loop over each point and add a marker to our map
for lat, lon in zip(latlons['lat'],latlons['lon']):
    folium.Marker([lat, lon]).add_to(map)
map

# Same as before but add a popup showing lat, lon values when clicking on a marker
# For popup, you can use HTML syntax
map = folium.Map(location=[65., 335.0], zoom_start=5 , tiles='Stamen Terrain')
for lat, lon in zip(latlons['lat'],latlons['lon']):
    folium.Marker([lat, lon], popup='<b>'+str(lat)+','+str(lon)+'</b>').add_to(map)
map

# Same as before but add a popup showing lat, lon values when clicking on a marker
# For popup, you can use HTML syntax
# CircleMarker can be customized too with colors and radius
map = folium.Map(location=[65., 335.0], zoom_start=5 , tiles='Stamen Terrain')

for lat, lon, age in zip(latlons['lat'],latlons['lon'], latlons['age']):
    if np.isnan(age):
        age=0.5
    folium.CircleMarker([lat, lon], popup='<b>'+str(lat)+','+str(lon)+'</b>', radius=age , color='black',
                    fill_color='red', fill_opacity=0.5, opacity=0.4).add_to(map)
map

import branca.colormap as cm
linear = cm.LinearColormap(
    ['green', 'yellow', 'red'],
    vmin=0, vmax=16
)

map = folium.Map(location=[65., 335.0], zoom_start=5 , tiles='Stamen Terrain')

for lat, lon, age in zip(latlons['lat'],latlons['lon'], latlons['age']):
    if np.isnan(age):
        age=0.5
    folium.CircleMarker([lat, lon], popup='<b>'+str(lat)+','+str(lon)+'</b>', radius=age , color='black',
                    fill_color=linear(age), fill_opacity=0.5, opacity=0.4).add_to(map)
map

map.save('map_iceland.html')

import ipyleaflet

map = ipyleaflet.Map(center=(65., 335.0), zoom=5, basemap=ipyleaflet.basemaps.Hydda.Full)
map

# Add Markers. There is no pop-up!
map = ipyleaflet.Map(center=(65., 335.0), zoom=5, basemap=ipyleaflet.basemaps.Hydda.Full)
for lat, lon in zip(latlons['lat'],latlons['lon']):
    marker=ipyleaflet.Marker(location=(lat, lon))
    map.add_layer(marker)

map

INSTANCE_ID=""
map = folium.Map(location=[60, 10.7], zoom_start=7, tiles=None)

folium.raster_layers.WmsTileLayer(
    url='https://services.sentinel-hub.com/ogc/wms/' + INSTANCE_ID + "?showlogo=0&time=2018-07-18/2018-07-20",
    layers='TRUE_COLOR',
    fmt='image/png',
    overlay=False,
    control=True,
).add_to(map)

map

opm_attr='<a href="https://github.com/openplanetary/opm/wiki/OPM-Basemaps" target="blank">OpenPlanetaryMap</a>'

map = folium.Map(attr=f'NASA/MOLA |{opm_attr}', max_zoom=10, no_wrap=True,
           tiles='https://s3-eu-west-1.amazonaws.com/whereonmars.cartodb.net/mola-color/{z}/{x}/{-y}.png')

map

map = folium.Map([0,0],attr='Lunar_LOLA_Basemap', max_zoom=8, no_wrap=True,
           tiles='https://webgis3.wr.usgs.gov/arcgis/rest/services/Lunar_LOLA_Basemap/MapServer/tile/{z}/{y}/{x}.jp2')
map

map = folium.Map([0,0],attr='LOLA/USGS', max_zoom=5, 
           tiles='https://s3.amazonaws.com/opmbuilder/301_moon/tiles/w/hillshaded-albedo/{z}/{x}/{y}.png'
              )
map

import folium
map = folium.Map(location=[41, -70], zoom_start=5, tiles=None)


folium.raster_layers.TileLayer(
    tiles='http://{s}.google.com/vt/lyrs=s&x={x}&y={y}&z={z}',
    attr='google',
    name='google maps',
    max_zoom=20,
    subdomains=['mt0', 'mt1', 'mt2', 'mt3'],
    overlay=False,
    control=True,
).add_to(map)

folium.raster_layers.TileLayer(
    tiles='http://{s}.google.com/vt/lyrs=m&x={x}&y={y}&z={z}',
    attr='google',
    name='google street view',
    max_zoom=20,
    subdomains=['mt0', 'mt1', 'mt2', 'mt3'],
    overlay=False,
    control=True,
).add_to(map)


folium.raster_layers.WmsTileLayer(
    url='https://demo.boundlessgeo.com/geoserver/ows?',
    layers='nasa:bluemarble',
    name='bluemarble',
    fmt='image/png',
    overlay=False,
    control=True,
).add_to(map)


folium.raster_layers.WmsTileLayer(
    url='http://mesonet.agron.iastate.edu/cgi-bin/wms/nexrad/n0r.cgi',
    name='test',
    fmt='image/png',
    layers='nexrad-n0r-900913',
    attr=u'Weather data © 2012 IEM Nexrad',
    transparent=True,
    overlay=True,
    control=True,
).add_to(map)

folium.LayerControl().add_to(map)


map

import glob
def getint(name):
    _, num = name.split('P_')
    num, _ = num.split('.')
    return int(num)

filenames = glob.glob("/opt/uio/deep_python/data/Grace/grd_files_to_interp/GypsumP_*.grd")

filenames = sorted(filenames, key=getint)

import numpy as np
import xarray as xr
for ix, filename in enumerate(filenames):
    print(ix, filename)
    dset = xr.open_dataset(filename)
    dset=dset.expand_dims('depth')
    dset['depth'] = ('depth', [float(getint(filename))])
    if ix > 0:
        tmp = xr.concat((newarray, dset), dim='depth')
        newarray = tmp
    else:
        newarray = dset
    dset.close()

newarray.to_netcdf("/opt/uio/deep_python/data/Grace/Gypsump_all.grd")

depths=np.linspace(100.,2650.,500)
print(depths)
lats=np.linspace(-90.,90,newarray.lat.shape[0]*2)
lons=np.linspace(-180.,180,newarray.lon.shape[0]*2)

newarray=newarray.interp(depth=depths, lat=lats, lon=lons)
newarray.to_netcdf("/opt/uio/deep_python/data/Grace/Gypsump_all_interp.grd")

newarray = xr.open_dataset("/opt/uio/deep_python/data/Grace/Gypsump_all_interp.grd")

img = np.zeros((newarray.depth.shape[0],newarray.lat.shape[0], newarray.lon.shape[0]))
img[:,:,:] = newarray.where((newarray.z<=-0.5) | (newarray.z >= 0.5))['z'][:,:]

bounds = np.array(
        [-180, 180, -90, 90, 1, 110]
    ).astype("int")
	
import k3d
volume = k3d.volume(img[0:200,::,::] , 
                    compression_level=9,height=1200, bounds = bounds, 
					color_range=[-4,4], color_map=np.array(k3d.basic_color_maps.Jet , dtype=np.float64))

plot = k3d.plot(height=1200)
plot += volume
plot.display()

import pandas as pd
import numpy as np

filename = "data/Isabel/Edmundson_data_python_DEEP.xlsx"
frame = pd.read_excel(filename)

frame.head()

import plotly.offline as py
import plotly.graph_objs as go

py.init_notebook_mode()

data_to_plot = [go.Scatter(x=frame.index, y=frame["HC_column"])]
py.iplot(data_to_plot)

# To add a title, etc.
layout = go.Layout(
    title='Test ',
    xaxis=dict(
        title='index',
        titlefont=dict(
            size=18,
            color='#7f7f7f'
        )
    ),
    yaxis=dict(
        title='HC_column',
        titlefont=dict(
            size=18,
            color='#7f7f7f'
        )
    )
)
fig = go.Figure(data=data_to_plot, layout=layout)
py.iplot(fig)

py.plot(fig, filename = 'plotly_example.html', auto_open=False)

import os
os.environ['PROJ_LIB']= "/opt/tljh/user/share/proj"
import plotly.offline as py
import cartopy as ccrs
import numpy as np
import xarray as xr
import numpy as np
import matplotlib as mpl
mpl.use('Agg')
from mpl_toolkits.basemap import Basemap
# For plotly
py.init_notebook_mode()

# convert degrees to radians
def degree2radians(degree):
    #convert degrees to radians
    return degree*np.pi/180.

# get depth from filename
def getint(name):
    _, num = name.split('P_')
    num, _ = num.split('.')
    return int(num)

# maps the points of coords (lon, lat) to points onto the  sphere of radius radius
def mapping_map_to_sphere(lon, lat, radius=1):
    lon=np.array(lon, dtype=np.float64)
    lat=np.array(lat, dtype=np.float64)
    lon=degree2radians(lon)
    lat=degree2radians(lat)
    xs=radius*np.cos(lon)*np.cos(lat)
    ys=radius*np.sin(lon)*np.cos(lat)
    zs=radius*np.sin(lat)
    return xs, ys, zs

# Functions converting coastline/country polygons to lon/lat traces
def polygons_to_traces(poly_paths, N_poly, m):
    ''' 
    pos arg 1. (poly_paths): paths to polygons
    pos arg 2. (N_poly): number of polygon to convert
    '''
    # init. plotting list
    lons=[]
    lats=[]

    for i_poly in range(N_poly):
        poly_path = poly_paths[i_poly]
        
        # get the Basemap coordinates of each segment
        coords_cc = np.array(
            [(vertex[0],vertex[1]) 
             for (vertex,code) in poly_path.iter_segments(simplify=False)]
        )
        
        # convert coordinates to lon/lat by 'inverting' the Basemap projection
        lon_cc, lat_cc = m(coords_cc[:,0],coords_cc[:,1], inverse=True)
    
        
        lats.extend(lat_cc.tolist()+[None]) 
        lons.extend(lon_cc.tolist()+[None])
        
       
    return lons, lats

# Function generating coastline lon/lat 
def get_coastline_traces(m):
    poly_paths = m.drawcoastlines().get_paths() # coastline polygon paths
    N_poly = 91  # use only the 91st biggest coastlines (i.e. no rivers)
    cc_lons, cc_lats= polygons_to_traces(poly_paths, N_poly, m)
    return cc_lons, cc_lats

# Function generating country lon/lat 
def get_country_traces(m):
    poly_paths = m.drawcountries().get_paths() # country polygon paths
    N_poly = len(poly_paths)  # use all countries
    country_lons, country_lats= polygons_to_traces(poly_paths, N_poly, m)
    return country_lons, country_lats

filename = "/opt/uio/deep_python/data/Grace/grd_files_to_interp/GypsumP_100.grd"

depth = getint(filename)

df = xr.open_dataset(filename)
df

lon=df.lon.values
lat=df.lat.values
z=df.z.values
# min and max values (used for plotting)
vmin=z.min()
vmax=z.max()

# Make shortcut to Basemap object, 
# not specifying projection type for this example
map = Basemap() 

cc_lons, cc_lats=get_coastline_traces(map)
country_lons, country_lats=get_country_traces(map)

#concatenate the lon/lat for coastlines and country boundaries:
lons=cc_lons+[None]+country_lons
lats=cc_lats+[None]+country_lats

xs, ys, zs=mapping_map_to_sphere(lons, lats, radius=1.01)# here the radius is slightly greater than 1 
                                                         #to ensure lines visibility; otherwise (with radius=1)
                                                         # some lines are hidden by contours colors
        
boundaries=dict(type='scatter3d',
               x=xs,
               y=ys,
               z=zs,
               mode='lines',
               line=dict(color='black', width=1)
                )

# Take lat, lon from field
clons=np.array(lon.tolist()+[180], dtype=np.float64)
clats=np.array(lat, dtype=np.float64)
# Create a meshgrid from lats and lons
clons, clats=np.meshgrid(clons, clats)

# Project our (lat,lon,z) to our sphere
XS, YS, ZS=mapping_map_to_sphere(clons, clats)

# generate values for the field Z
# Wrap over longitudes
nrows, ncolumns=clons.shape
VALUES=np.zeros(clons.shape, dtype=np.float64)
VALUES[:, :ncolumns-1]=np.copy(np.array(z,  dtype=np.float64))
VALUES[:, ncolumns-1]=np.copy(z[:, 0])

# Create a sphere: colors are taken from our field Z that was mapped to our sphere(VALUES)
sphere=dict(type='surface',
            x=XS, 
            y=YS, 
            z=ZS,
            colorscale='Jet',
            surfacecolor=VALUES,
            cmin=vmin, 
            cmax=vmax,
            colorbar=dict(thickness=20, len=0.75, ticklen=4, title= 'Z'),
            hoverinfo='text')

# Turn off all axis
noaxis=dict(showbackground=False,
            showgrid=False,
            showline=False,
            showticklabels=False,
            ticks='',
            title='',
            zeroline=False)

# Define layout i.e. title, figure size, axis (x, y, z), camera, etc.
layout3d=dict(title='Z<br>Depth '+ str(getint(filename)) + 'm',
              font=dict(family='Balto', size=14),
              width=700, 
              height=700,
              scene=dict(xaxis=noaxis, 
                         yaxis=noaxis, 
                         zaxis=noaxis,
                         aspectratio=dict(x=1,
                                          y=1,
                                          z=1),
                         camera=dict(eye=dict(x=1.15, 
                                     y=1.15, 
                                     z=1.15)
                                    )
            ),
            paper_bgcolor='rgba(255,255,255, 1.0)'  
           )

fig=dict(data=[sphere, boundaries], layout=layout3d)
py.iplot(fig, filename='z-map2sphere')

import ipyvolume as ipv

Z = df.z.values.astype('float64')
from matplotlib import cm
colormap = cm.jet

znorm = VALUES - VALUES.min()
znorm /= znorm.ptp()
znorm.min(), znorm.max()
color = colormap(znorm)
ipv.figure()
mesh = ipv.plot_surface(XS, YS, ZS, color=color[...,:3])
countries = ipv.pylab.plot(xs,ys,zs, color='black')
ipv.pylab.style.box_off()
ipv.pylab.style.axes_off()
ipv.show()
ipv.pylab.save("ipyvolume_sphere3D.html")