#!/usr/bin/env python
# coding: utf-8

# * This notebook was to modify the banks area in the SoG into land

# In[1]:


from salishsea_tools import (bathy_tools, nc_tools,viz_tools,tidetools)
from numpy import loadtxt
from __future__ import division, print_function
from salishsea_tools.nowcast import figures
import scipy.io as sio
import matplotlib.pyplot as plt
import netCDF4 as nc
import numpy as np
get_ipython().run_line_magic('matplotlib', 'inline')


# In[2]:


get_ipython().system('cp /ocean/jieliu/research/meopar/river-treatment/bathy_meter_SalishSea6.nc  bathy_meter_SalishSea13.nc')


# In[3]:


bathy = nc.Dataset('/ocean/jieliu/research/meopar/river-treatment/bathy_meter_SalishSea6.nc', 'r')
bathyy = bathy.variables['Bathymetry'][:]


# In[4]:


bathy_9 = nc.Dataset('bathy_meter_SalishSea9.nc','r+' )
bathyy_9 = bathy_9.variables['Bathymetry']                                            
lon_9 = bathy_9.variables['nav_lon']
lat_9 = bathy_9.variables['nav_lat']


# In[5]:


fig = bathy_tools.plot_colourmesh(
    bathy, 'NEMO - Salish Sea Bathymetry', 
    axis_limits=(-123.4, -122.7, 49., 49.7), colour_map='spectral', bins=30)


# In[6]:


bath, X, Y = tidetools.get_bathy_data(bathy)
fig, ax = plt.subplots(1, 1, figsize=(10, 8))
viz_tools.set_aspect(ax)
cmap = plt.get_cmap('winter_r')
cmap.set_bad('burlywood')
mesh = ax.pcolormesh(bath, cmap=cmap)
fig.colorbar(mesh)
plt.axis((280, 340, 370, 460))


# In[8]:


## to put the range of lat, lon indices and bathymetry into list
lat_ind = []
lon_ind = []
remain_edit_bathy = []
for j in np.arange(370, 460):
    for i in np.arange(280, 340):
        if  (bathyy[j,i] <=4):
            #lat_ind.append(j)
            #lon_ind.append(i)
            #remain_edit_bathy.append(bathyy_9[j,i])
            bathyy_9[j,i] = np.ma.masked_values(bathyy_9[j,i], 4)
            #print (j,i)


# In[9]:


bathy_9.close()


# In[47]:


a= nc.Dataset('bathy_meter_SalishSea9.nc')


# In[48]:


b= a.variables['Bathymetry']
b[427,288]


# #Create bathymetry3 with no banks of short river for Nancy 

# In[10]:


get_ipython().system('cp /ocean/jieliu/research/meopar/nemo-forcing/grid/bathy_meter_SalishSea2.nc  bathy_meter_SalishSea3.nc')


# In[11]:


bathy = nc.Dataset('/ocean/jieliu/research/meopar/nemo-forcing/grid/bathy_meter_SalishSea2.nc', 'r')
bathyy = bathy.variables['Bathymetry'][:]


# In[12]:


bathy_3 = nc.Dataset('bathy_meter_SalishSea3.nc','r+' )
bathyy_3 = bathy_3.variables['Bathymetry']                                            
lon_3 = bathy_3.variables['nav_lon']
lat_3 = bathy_3.variables['nav_lat']


# In[13]:


bath, X, Y = tidetools.get_bathy_data(bathy)
fig, ax = plt.subplots(1, 1, figsize=(10, 8))
viz_tools.set_aspect(ax)
cmap = plt.get_cmap('winter_r')
cmap.set_bad('burlywood')
mesh = ax.pcolormesh(bath, cmap=cmap)
fig.colorbar(mesh)
plt.axis((280, 310, 350, 480))


# In[14]:


## to put the range of lat, lon indices and bathymetry into list
lat_ind = []
lon_ind = []
remain_edit_bathy = []
for j in np.arange(350, 480):
    for i in np.arange(280, 310):
        if  (bathyy_3[j,i] <=4):
            #lat_ind.append(j)
            #lon_ind.append(i)
            #remain_edit_bathy.append(bathyy_3[j,i])
            bathyy_3[j,i] = np.ma.masked_values(bathyy_3[j,i], 4)
            #print (j,i)


# In[15]:


fig = bathy_tools.plot_colourmesh(
    bathy_3, 'NEMO - Salish Sea Bathymetry', 
    axis_limits=(-123.4, -122.7, 49., 49.7), colour_map='spectral', bins=30)


# In[16]:


## Long river channel without banks
bathy_9 = nc.Dataset('bathy_meter_SalishSea9.nc' )
fig = bathy_tools.plot_colourmesh(
    bathy_9, 'NEMO - Salish Sea Bathymetry', 
    axis_limits=(-123.4, -122.7, 49., 49.7), colour_map='spectral', bins=30)


# In[17]:


bathy_3.close()


# In[ ]: