import matplotlib.pyplot as plt
#from mpl_toolkits.basemap import Basemap, cm
import matplotlib.mlab as mlab
import numpy as np
import pandas as pd
import numpy.ma as ma
import matplotlib.colors as colors
#from haversine import haversine
from matplotlib.colors import LinearSegmentedColormap
import cartopy.crs as ccrs
import matplotlib
import xarray as xr
import matplotlib.gridspec as gridspec # GRIDSPEC !
from matplotlib.colorbar import Colorbar
#import mplcursors
import cartopy.crs as ccrs
import cartopy
#import fiona
import shapely.geometry as sgeom
from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER
from matplotlib.image import BboxImage
from matplotlib.transforms import Bbox, TransformedBbox
from cartopy.feature import ShapelyFeature
from cartopy.io.shapereader import Reader
import matplotlib.ticker as mticker
from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER
from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter

def createMapProjections(lat0, lon0, region='Whole', fontsize=8, annotate=True, ax=None, linewidth=0.2, draw_grounding_line=True, regionLonLims=None):
    #m = Basemap(projection='hammer',lon_0=270)
    # plot just upper right quadrant (corners determined from global map).
    # keywords llcrnrx,llcrnry,urcrnrx,urcrnry used to define the lower
    # left and upper right corners in map projection coordinates.
    # llcrnrlat,llcrnrlon,urcrnrlon,urcrnrlat could be used to define
    # lat/lon values of corners - but this won't work in cases such as this
    # where one of the corners does not lie on the earth.
    #m1 = Basemap(projection='ortho',lon_0=-9,lat_0=60,resolution=None)
    
    #m = Basemap(projection='cyl', llcrnrlon=-120, llcrnrlat=-80, urcrnrlon=55, 
    #            urcrnrlat=-30, lat_0 = -74, lon_0 =-43);
    matplotlib.rcParams.update({'font.size': fontsize})    
    m1 = Basemap(projection='ortho', lat_0 =lat0, lon_0 =lon0, resolution='l');
    width = m1.urcrnrx - m1.llcrnrx
    height = m1.urcrnry - m1.llcrnry
    dist_x = 8199.491128244028
    w_factor = 0.3
    h_factor = 0.3
    if(region =='Whole'):
        m  = Basemap(projection='ortho', lat_0=lat0, lon_0=lon0, resolution='h', llcrnrx=-width*w_factor, llcrnry=-height*h_factor, urcrnrx=w_factor*width, urcrnry=h_factor*height)
        ## regionNames = ['Weddell Sea (WS)', 'Ross Sea (RS)', 'Prydz Bay (PB)', 'Cape Darnley (CD)', 'Knox Coast (KC)', 'Adelie Coast (AC)', 'Amundsen Sea (AS)', 'Belingshausen Sea (BS)']
        ## regionCoord = [[-45, -65, 40], [180 , -70, 0], [75, -50, -30], [65, -50, -25], [100, -55, 10], [140, -55, 40], [-100, -60, -10], [-90, -60, 90]]
        ## for i in range(len(regionNames)):
        ##     mxy = m(regionCoord[i][0], regionCoord[i][1] )
        ##     plt.annotate(regionNames[i], xy=(mxy[0], mxy[1]), rotation=regionCoord[i][2], fontsize=12)
    elif(region == 'Weddell'):
        m  = Basemap(projection='ortho', lat_0=lat0, lon_0=lon0, resolution='h', llcrnrx=-width*0.15, llcrnry=0.05*height, urcrnrx=0.*width, urcrnry=0.18*height)
    elif(region == 'Ross'):
        m  = Basemap(projection='ortho', lat_0=lat0, lon_0=180, resolution='h', llcrnrx=-width*0.075, llcrnry=height*0.08, urcrnrx=width*0.03, urcrnry=height*0.18)
    elif(region == 'CDP'):
        m = Basemap(projection='ortho', lat_0=lat0, lon_0=65, resolution='h', llcrnrx=-width*0.05, llcrnry=height*0.175, urcrnrx=width*0.05, urcrnry=height*0.23)
    elif(region == 'Prydz'):
        m = Basemap(projection='ortho', lat_0=lat0, lon_0=75, resolution='h', llcrnrx=-width*0.1, llcrnry=height*0.15, urcrnrx=width*0.07, urcrnry=height*0.25)
    elif(region == 'Knox'):
        m = Basemap(projection='ortho', lat_0=lat0, lon_0=106.5, resolution='h', llcrnrx=-width*0.07, llcrnry=height*0.17, urcrnrx=width*0.07, urcrnry=height*0.25)
    elif(region == 'Adelie'):
        m = Basemap(projection='ortho', lat_0=lat0, lon_0=140, resolution='h', llcrnrx=-width*0.07, llcrnry=height*0.17, urcrnrx=width*0.07, urcrnry=height*0.25)
    elif(region == 'Amundsen'):
        m = Basemap(projection='ortho', lat_0=lat0, lon_0=-110, resolution='h', llcrnrx=-width*0.1, llcrnry=height*0.12, urcrnrx=width*0.07, urcrnry=height*0.18)
    elif(region == 'Belingshausen'):
        m = Basemap(projection='ortho', lat_0=lat0, lon_0=-75, resolution='h', llcrnrx=-width*0.1, llcrnry=height*0.12, urcrnrx=width*0.07, urcrnry=height*0.25)
    elif(region == 'Global'):
        m = m1

    #m.drawmapboundary(linewidth=linewidth);
    if(draw_grounding_line == True):
        m.readshapefile("/media/data/Datasets/Shapefiles/AntarcticGroundingLine/GSHHS_f_L6", "GSHHS_f_L6", color='0.75', linewidth=0.1)
    #m.fillcontinents(color='#ddaa66');
    m.drawcoastlines(linewidth=linewidth)    

    if regionLonLims:
        ## Draw the regions longitudinal limits
        m.drawmeridians(regionLonLims, labels=[0,0,0,0], linewidth=0.5, color='b')    # labels: left,right,top,bottom

    if(annotate == True):
        parallels = np.arange(-80, -50+1, 5.)    
        m.drawparallels(parallels,labels=[0,0,0,0], linewidth=0.2) # labels: left,right,top,bottom
        meridians = np.arange(-180, 180, 20.)
        m.drawmeridians(meridians,labels=[1,1,1,1], linewidth=0.2)
        
        xy = [[-140, -55], [-140, -60] , [-140, -65], [-140, -70], [-140, -75] , [-140, -80]]
        xytext = np.arange(55, 81, 5)
        for i in range(len(xytext)):
            mxy = m(xy[i][0], xy[i][1])
            plt.annotate(str(xytext[i])+"$^o$S", xy=(mxy[0], mxy[1]), rotation=-45, fontsize=fontsize)
    return m


def createMapProjections_CRS(lat0, lon0, region='Whole', fontsize=8, annotate=True, ax=None, linewidth=0.2, draw_grounding_line=True, regionLonLims=None):
    matplotlib.rcParams.update({'font.size': fontsize})

    if(region =='Whole'):
        m = ccrs.Orthographic(central_longitude=lon0, central_latitude=lat0)
    elif(region == 'Weddell'):
        m  = Basemap(projection='ortho', lat_0=lat0, lon_0=lon0, resolution='h', llcrnrx=-width*0.15, llcrnry=0.05*height, urcrnrx=0.*width, urcrnry=0.18*height)
    elif(region == 'Ross'):
        m  = Basemap(projection='ortho', lat_0=lat0, lon_0=180, resolution='h', llcrnrx=-width*0.075, llcrnry=height*0.08, urcrnrx=width*0.03, urcrnry=height*0.18)
    elif(region == 'CDP'):
        m = ccrs.PlateCarree()
    elif(region == 'Prydz'):
        m = Basemap(projection='ortho', lat_0=lat0, lon_0=75, resolution='h', llcrnrx=-width*0.1, llcrnry=height*0.15, urcrnrx=width*0.07, urcrnry=height*0.25)
    elif(region == 'Knox'):
        m = Basemap(projection='ortho', lat_0=lat0, lon_0=106.5, resolution='h', llcrnrx=-width*0.07, llcrnry=height*0.17, urcrnrx=width*0.07, urcrnry=height*0.25)
    elif(region == 'Adelie'):
        m = Basemap(projection='ortho', lat_0=lat0, lon_0=140, resolution='h', llcrnrx=-width*0.07, llcrnry=height*0.17, urcrnrx=width*0.07, urcrnry=height*0.25)
    elif(region == 'Amundsen'):
        m = Basemap(projection='ortho', lat_0=lat0, lon_0=-110, resolution='h', llcrnrx=-width*0.1, llcrnry=height*0.12, urcrnrx=width*0.07, urcrnry=height*0.18)
    elif(region == 'Belingshausen'):
        m = Basemap(projection='ortho', lat_0=lat0, lon_0=-75, resolution='h', llcrnrx=-width*0.1, llcrnry=height*0.12, urcrnrx=width*0.07, urcrnry=height*0.25)
    elif(region == 'Global'):
        m = m1

    #m.drawmapboundary(linewidth=linewidth);
    if(draw_grounding_line == True):
        m.readshapefile("/media/data/Datasets/Shapefiles/AntarcticGroundingLine/GSHHS_f_L6", "GSHHS_f_L6", color='0.75', linewidth=0.1)
    #m.fillcontinents(color='#ddaa66');
    m.drawcoastlines(linewidth=linewidth)    

    if regionLonLims:
        ## Draw the regions longitudinal limits
        m.drawmeridians(regionLonLims, labels=[0,0,0,0], linewidth=0.5, color='b')    # labels: left,right,top,bottom

    if(annotate == True):
        parallels = np.arange(-80, -50+1, 5.)    
        m.drawparallels(parallels,labels=[0,0,0,0], linewidth=0.2) # labels: left,right,top,bottom
        meridians = np.arange(-180, 180, 20.)
        m.drawmeridians(meridians,labels=[1,1,1,1], linewidth=0.2)
        
        xy = [[-140, -55], [-140, -60] , [-140, -65], [-140, -70], [-140, -75] , [-140, -80]]
        xytext = np.arange(55, 81, 5)
        for i in range(len(xytext)):
            mxy = m(xy[i][0], xy[i][1])
            plt.annotate(str(xytext[i])+"$^o$S", xy=(mxy[0], mxy[1]), rotation=-45, fontsize=fontsize)
    return m



def getCellSize(m, nx=300, ny=300):
    width = m.urcrnrx - m.llcrnrx
    height = m.urcrnry - m.llcrnry
    left = [m.llcrnrx, m.llcrnry+height*0.5]
    right = [m.urcrnrx, m.urcrnry-height*0.5]
    bottom = [m.llcrnrx+width*0.5, m.llcrnry]
    top = [m.llcrnrx+width*0.5, m.urcrnry]
    
    left_lonlat = m(left[0], left[1], inverse=True)
    right_lonlat = m(right[0], right[1], inverse=True)
    top_lonlat = m(top[0], top[1], inverse=True)
    bot_lonlat = m(bottom[0], bottom[1], inverse=True)
    dist_x = haversine(list(left_lonlat[::-1]), list(right_lonlat[::-1]))
    dist_y = haversine(top_lonlat[::-1], bot_lonlat[::-1])

    return dist_x/float(nx), dist_y/float(ny)
    
def plotBotVarContourf(df,var="PSAL_ADJUSTED", units='PSU', cmin=33, cmax=35.5,
                         save=False, savename="savedFig.png", wd=7, ht=7, cmap='viridis', region='Whole', nmin=0, show=True, cx=10, cy=10, levs=[]):

    matplotlib.rcParams.update({'font.size': 8})        # setting fontsize for plot elements            
    plt.figure(1, figsize=(wd,ht));
    gs = gridspec.GridSpec(1, 2, width_ratios=[1, 0.05])
    mapax = plt.subplot(gs[0, 0])

    lat0 = -90
    lon0 = 0
    m = createMapProjections(lat0, lon0, region=region)
    width = m.urcrnrx - m.llcrnrx
    height = m.urcrnry - m.llcrnry
    left = [m.llcrnrx, m.llcrnry+height*0.5]
    right = [m.urcrnrx, m.urcrnry-height*0.5]
    bottom = [m.llcrnrx+width*0.5, m.llcrnry]
    top = [m.llcrnrx+width*0.5, m.urcrnry]
    
    left_lonlat = m(left[0], left[1], inverse=True)
    right_lonlat = m(right[0], right[1], inverse=True)
    top_lonlat = m(top[0], top[1], inverse=True)
    bot_lonlat = m(bottom[0], bottom[1], inverse=True)
    dist_x = haversine(list(left_lonlat[::-1]), list(right_lonlat[::-1]))
    dist_y = haversine(top_lonlat[::-1], bot_lonlat[::-1])
    nx = dist_x // cx
    ny = dist_x // cy
    print(nx, ny)

    botIndex = df.groupby('PROFILE_NUMBER').tail(1).index
    
    z = df.loc[botIndex, var].values
    xlons, ylats = m(df.loc[botIndex , 'LONGITUDE'].values, df.loc[botIndex,'LATITUDE'].values)

    Xgrid = np.linspace(m.llcrnrx, m.urcrnrx, nx)
    Ygrid = np.linspace(m.llcrnry, m.urcrnry, ny)
    XX, YY = np.meshgrid(Xgrid, Ygrid)
    zi = np.zeros((len(Ygrid), len(Xgrid)))
    ni = np.zeros((len(Ygrid), len(Xgrid)))
    
    for i in range(len(z)):
        if(np.isnan(z[i]) == True):
            pass
        else:
            zi[np.argmin(np.abs(Ygrid - ylats[i])), np.argmin(np.abs(Xgrid - xlons[i]))] += z[i]
            ni[np.argmin(np.abs(Ygrid - ylats[i])), np.argmin(np.abs(Xgrid - xlons[i]))] += 1

    #wz = np.where(zi == 0) # where zero
    #ni[wz[0], wz[1]] = 1
    zi[:, :] = zi[:, :]/ni[:, :]
    #zi = ma.array(zi)
    #zi_masked = ma.masked_where(zi == 0, zi)
    #ni_masked = ma.masked_less(ni, nmin)
    #zi_masked = ma.masked_array(data=zi, mask=(zi_masked.mask|ni_masked.mask))
    
    m.drawmapboundary();
    m.readshapefile("/media/data/Datasets/Shapefiles/AntarcticGroundingLine/GSHHS_f_L6", "GSHHS_f_L6", color='0.75', linewidth=0.1)
    m.drawcoastlines(linewidth=0.2)

    if not levs:
        levs = np.arange(cmin, cmax, 0.2)
    CF = m.contourf(XX, YY, zi, vmin=cmin, vmax=cmax, cmap=cmap, levels=levs, extend='max', ax=mapax)
    colorbarax = plt.subplot(gs[0,1])
    cbar = Colorbar(ax=colorbarax, mappable=CF, ticks=levs) 
    cbar.set_label(units)

    if(save==True):
        plt.savefig(savename, dpi=300, bbox_inches='tight')

    if(show == True):
        plt.show();
    else:
        plt.close();


def plotSurfVarContourf(df,var="PSAL_ADJUSTED", units='Cond.', cmin=33, cmax=35.5,
                         save=False, savename="savedFig.png",
                         wd=7, ht=7, cx=10, cy=10, cmap='viridis', nmin=0, region='Whole', levels=[0], show=True):
    plt.figure(1,figsize=(wd,ht));
            
    lat0 = -89
    lon0 = 0
    m  = createMapProjections(lat0, lon0, region=region)
    width = m.urcrnrx - m.llcrnrx
    height = m.urcrnry - m.llcrnry
    left = [m.llcrnrx, m.llcrnry+height*0.5]
    right = [m.urcrnrx, m.urcrnry-height*0.5]
    bottom = [m.llcrnrx+width*0.5, m.llcrnry]
    top = [m.llcrnrx+width*0.5, m.urcrnry]
    
    left_lonlat = m(left[0], left[1], inverse=True)
    right_lonlat = m(right[0], right[1], inverse=True)
    top_lonlat = m(top[0], top[1], inverse=True)
    bot_lonlat = m(bottom[0], bottom[1], inverse=True)
    dist_x = haversine(list(left_lonlat[::-1]), list(right_lonlat[::-1]))
    dist_y = haversine(top_lonlat[::-1], bot_lonlat[::-1])
    nx = dist_x // cx
    ny = dist_x // cy
    print(nx, ny)
    surfIndex = df.groupby('PROFILE_NUMBER').head(1).index
    
    z = df.loc[surfIndex, var].values
    xlons, ylats = m(df.loc[surfIndex , 'LONGITUDE'].values, df.loc[surfIndex,'LATITUDE'].values)

    Xgrid = np.linspace(m.llcrnrx, m.urcrnrx, nx)
    Ygrid = np.linspace(m.llcrnry, m.urcrnry, ny)
    XX, YY = np.meshgrid(Xgrid, Ygrid)
    zi = np.zeros((len(Ygrid), len(Xgrid)))
    ni = np.zeros((len(Ygrid), len(Xgrid)))
    
    for i in range(len(z)):
        zi[np.argmin(np.abs(Ygrid - ylats[i])), np.argmin(np.abs(Xgrid - xlons[i]))] += z[i]
        ni[np.argmin(np.abs(Ygrid - ylats[i])), np.argmin(np.abs(Xgrid - xlons[i]))] += 1

    wz = np.where(zi == 0) # where zero
    ni[wz[0], wz[1]] = 1 #setting 0 counts to 1 to avoid div by zero
    zi[:, :] = zi[:, :]/ni[:, :] #averaging over the no of samples from each grid cell
    zi = ma.array(zi)
    zi_masked = ma.masked_where(zi == 0, zi)
    ni_masked = ma.masked_less(ni, nmin)
    zi_masked = ma.masked_array(data=zi, mask=(zi_masked.mask|ni_masked.mask))

    if(len(levels) == 1):
        ticks = np.arange(cmin, cmax, 0.2)
    else:
        ticks = levels

    CF = m.contourf(XX, YY, zi_masked, vmin=cmin, vmax=cmax, cmap=cmap, levels=ticks, extend='both'); #, cmap='viridis'
    cbar = m.colorbar(CF, pad='12%', boundaries=ticks, ticks=ticks, spacing='uniform') 
    cbar.set_label(units)
    
    if(save==True):
        plt.savefig(savename, dpi=300)
    if(show == True):
        plt.show();
    else:
        plt.close();


def plotDataDensity(df, units='Data Density',
                    save=False, savename="savedFig.png", wd=7, ht=7,
                    cx=10, cy=10, show=False, lat0 = -90, lon0=0, regionLonLims = [ -60. ,  -20. ,  0., 29., 37., 60. ,   70. , 75,  82. ,  87., 101. ,  112. ,
        135. ,  145. ,  160. ,  180. , -120. , -100.],
                    levels=[0, 10, 20, 30, 40, 50, 60, 100, 200, 500], region='Whole', plotBathy=False, fontsize=8,
                    bathyLevels=[0, -100, -200, -300, -400, -500, -600, -700, -800, -900, -1000, -1250, -1500, -1750, -2000, -2500, -3000, -3500, -4000, -4500, -5000, -6000, -7000]):
    matplotlib.rcParams.update({'font.size': fontsize})        # setting fontsize for plot elements            
    plt.figure(1, figsize=(wd,ht));
    gs = gridspec.GridSpec(2, 2, height_ratios=[1, 0.15], width_ratios=[1, 0.03])
    mapax = plt.subplot(gs[0:2, 0])
    m  = createMapProjections(lat0, lon0, region=region, fontsize=fontsize, regionLonLims= regionLonLims)

        
    datacolorbar = plt.subplot(gs[0, 1])
    #bathycolorbar = plt.subplot(gs[1, 1])
    
    matplotlib.rcParams.update({'font.size': fontsize})    
    

    width = m.urcrnrx - m.llcrnrx
    height = m.urcrnry - m.llcrnry
    left = [m.llcrnrx, m.llcrnry+height*0.5]
    right = [m.urcrnrx, m.urcrnry-height*0.5]
    bottom = [m.llcrnrx+width*0.5, m.llcrnry]
    top = [m.llcrnrx+width*0.5, m.urcrnry]
    
    left_lonlat = m(left[0], left[1], inverse=True)
    right_lonlat = m(right[0], right[1], inverse=True)
    top_lonlat = m(top[0], top[1], inverse=True)
    bot_lonlat = m(bottom[0], bottom[1], inverse=True)
    dist_x = haversine(list(left_lonlat[::-1]), list(right_lonlat[::-1]))
    dist_y = haversine(top_lonlat[::-1], bot_lonlat[::-1])
    nx = dist_x // cx
    ny = dist_x // cy
    #lat_0 = -60, lon_0 = -20,

    surfIndex = df.groupby('PROFILE_NUMBER').head(1).index
    
    xlons, ylats = m(df.loc[surfIndex , 'LONGITUDE'].values, df.loc[surfIndex,'LATITUDE'].values)

    Xgrid = np.linspace(m.llcrnrx, m.urcrnrx, nx)
    Ygrid = np.linspace(m.llcrnry, m.urcrnry, ny)
    XX, YY = np.meshgrid(Xgrid, Ygrid)
    ni = np.zeros((len(Ygrid), len(Xgrid)))
    
    if(plotBathy == True):
        bathy = xr.open_dataset('/media/data/Datasets/Bathymetry/GEBCO_2014_2D.nc')
        lonlen = len(bathy.lon)
        lonindices = np.arange(0, lonlen+1, 30)
        lonindices[-1] = lonindices[-1] - 1
        bathyS = bathy.isel(lon=lonindices, lat=np.arange(0, 3600, 5))
        clevs = np.array([-1000, -2000, -3000])[::-1]
        
        longrid, latgrid = np.meshgrid(bathyS.lon.values, bathyS.lat.values)
        cs = m.contour(longrid, latgrid, bathyS.elevation.where(bathyS.elevation <= 0).values,  latlon=True, levels=[-1000], colors="magenta", linestyle=":", linewidths=0.35, extend='neither', ax=mapax) #, , cmap='rainbow'   , levels=clevs,
        ## plt.figure(2)
        ## cf = plt.contourf(longrid, latgrid,bathyS.elevation.where(bathyS.elevation <= 0).values, levels=clevs, extend='min') #, , cmap='rainbow'   , levels=clevs,
        ## plt.figure(1)
        
        #cbar1 = Colorbar(ax = bathycolorbar, mappable = cs, orientation = 'vertical')
        #cbar1.ax.get_children()[0].set_linewidths(5)
        #cbar1.set_label('Depth (m)')
        
    for i in range(len(xlons)):
        row = np.argmin(np.abs(Ygrid - ylats[i]))
        col = np.argmin(np.abs(Xgrid - xlons[i]))
        ni[row , col] = ni[row, col] + 1
        
    ni = ma.array(ni)
    ni_masked = ma.masked_equal(ni ,0)

    #zi_masked = ma.masked_array(zi, 0)
    ticks = levels
    from matplotlib.colors import BoundaryNorm
    bnorm = BoundaryNorm(levels, ncolors=len(levels)-1, clip=False)

    cmap = LinearSegmentedColormap.from_list(name='linearCmap', colors=['mistyrose', 'salmon', 'darkred'], N=len(levels)-1) 
        
    CF = m.contourf(XX, YY, ni_masked, vmin=min(levels), vmax=max(levels), levels=levels, norm=bnorm,
                    cmap=cmap, extend='max', ax=mapax); #, cmap='viridis'
    cbar = Colorbar(ax = datacolorbar, mappable = CF, ticks=ticks, orientation='vertical') #boundaries=levels, 
    cbar.set_label("Data density")


    if(save==True):
        plt.savefig(savename, dpi=300, bbox_inches='tight')
    if(show == True):
        plt.show();
    else:
        plt.close();


def plotDepthOfDive(df, units='Data Density',
                    save=False, savename="savedFig.png", wd=7, ht=7,
                    nx=820, ny=820, deltac=50, show=False,
                    levels=[0, -50, -100, -200, -300, -400, -500, -750, -1000, -1500, -2000], nmin=30):
    
    plt.figure(figsize=(wd,ht));
    lat0 = -89
    lon0 = 0
            
    m  = createMapProjections(lat0, lon0)
    
    #lat_0 = -60, lon_0 = -20,

    botIndex = df.groupby('PROFILE_NUMBER').tail(1).index
    
    xlons, ylats = m(df.loc[botIndex , 'LONGITUDE'].values, df.loc[botIndex,'LATITUDE'].values)
    z = df.loc[botIndex, 'DEPTH'].values
    Xgrid = np.linspace(m.llcrnrx, m.urcrnrx, nx)
    Ygrid = np.linspace(m.llcrnry, m.urcrnry, ny)
    XX, YY = np.meshgrid(Xgrid, Ygrid)
    zi = np.zeros((len(Ygrid), len(Xgrid)))
    ni = np.zeros((len(Ygrid), len(Xgrid)))
        
    for i in range(len(xlons)):
        row = np.argmin(np.abs(Ygrid - ylats[i]))
        col = np.argmin(np.abs(Xgrid - xlons[i]))
        zi[row , col] = z[i]
        ni[row, col] += 1

    zi = ma.array(zi)
    zi_masked = ma.masked_equal(zi , 0)
    ni_masked = ma.masked_less(ni, nmin)
    zi_masked = ma.masked_array(data=zi, mask=(zi_masked.mask|ni_masked.mask))

    levels= levels[::-1]
    ticks = levels
    from matplotlib.colors import BoundaryNorm
    bnorm = BoundaryNorm(levels, ncolors=len(levels)-1, clip=False)

    cmap = LinearSegmentedColormap.from_list(name='linearCmap', colors=['darkblue', 'CornFlowerBlue', 'w'], N=len(levels)-1)
    CF = m.contourf(XX, YY, zi_masked, vmin=min(levels), vmax=0, levels=levels, norm=bnorm,
                    cmap=cmap, extend='max'); #, cmap='viridis'
    cbar = m.colorbar(CF, pad=0.6, ticks=ticks, spacing='uniform') #boundaries=levels, 
    
    if(save==True):
        plt.savefig(savename, dpi=300)
    if(show == True):
        plt.show();
    else:
        plt.close();


def plotMaxVarContourf(df,var="PSAL_ADJUSTED", units='Cond.', cmin=33, cmax=35.5,
                         save=False, savename="savedFig.png", show=True,
                         wd=7, ht=7, cx=10, cy=10, cmap='viridis', nmin=30, region='Whole', levels=[0], isvar=True):
    plt.figure(1,figsize=(wd,ht));
            
    lat0 = -89
    lon0 = 0
    m  = createMapProjections(lat0, lon0, region=region)
    width = m.urcrnrx - m.llcrnrx
    height = m.urcrnry - m.llcrnry
    left = [m.llcrnrx, m.llcrnry+height*0.5]
    right = [m.urcrnrx, m.urcrnry-height*0.5]
    bottom = [m.llcrnrx+width*0.5, m.llcrnry]
    top = [m.llcrnrx+width*0.5, m.urcrnry]
    
    left_lonlat = m(left[0], left[1], inverse=True)
    right_lonlat = m(right[0], right[1], inverse=True)
    top_lonlat = m(top[0], top[1], inverse=True)
    bot_lonlat = m(bottom[0], bottom[1], inverse=True)
    dist_x = haversine(list(left_lonlat[::-1]), list(right_lonlat[::-1]))
    dist_y = haversine(top_lonlat[::-1], bot_lonlat[::-1])
    nx = dist_x // cx
    ny = dist_x // cy
    print(nx, ny)
    
    surfIndex = df.groupby('PROFILE_NUMBER')[var].transform(max) == df[var]
    
    if(isvar == False):
        var = 'DEPTH'
        levels= np.array([0, -100, -200, -300, -400, -500, -600, -700, -800])[::-1]
    
    z = df.loc[surfIndex, var].values
    xlons, ylats = m(df.loc[surfIndex , 'LONGITUDE'].values, df.loc[surfIndex,'LATITUDE'].values)

    Xgrid = np.linspace(m.llcrnrx, m.urcrnrx, nx)
    Ygrid = np.linspace(m.llcrnry, m.urcrnry, ny)
    XX, YY = np.meshgrid(Xgrid, Ygrid)
    zi = np.zeros((len(Ygrid), len(Xgrid)))
    ni = np.zeros((len(Ygrid), len(Xgrid)))
    
    for i in range(len(z)):
        zi[np.argmin(np.abs(Ygrid - ylats[i])), np.argmin(np.abs(Xgrid - xlons[i]))] += z[i]
        ni[np.argmin(np.abs(Ygrid - ylats[i])), np.argmin(np.abs(Xgrid - xlons[i]))] += 1

    wz = np.where(zi == 0) # where zero
    ni[wz[0], wz[1]] = 1 #setting 0 counts to 1 to avoid div by zero
    zi[:, :] = zi[:, :]/ni[:, :] #averaging over the no of samples from each grid cell
    zi = ma.array(zi)
    zi_masked = ma.masked_where(zi == 0, zi)
    ni_masked = ma.masked_less(ni, nmin)
    zi_masked = ma.masked_array(data=zi, mask=(zi_masked.mask|ni_masked.mask))

    
    if(len(levels) == 1):
        ticks = np.arange(cmin, cmax, 0.2)
    else:
        ticks = levels
        cmin = min(levels)
        cmax = max(levels)

    CF = m.contourf(XX, YY, zi_masked, vmin=cmin, vmax=cmax, cmap=cmap, levels=ticks, extend='both'); #, cmap='viridis'
    cbar = m.colorbar(CF, pad='12%', boundaries=ticks, ticks=ticks, spacing='uniform') 
    cbar.set_label(units)
    
    if(save==True):
        plt.savefig(savename, dpi=300)
    if(show == True):
        plt.show();
    else:
        plt.close();


def plotSeaIceConc(ICE, timeind, latmin=-90, latmax=-60, lonmin=-180, lonmax=180, save=False, savename="savedFig.png", show=False,
                         wd=7, ht=7, nx=820, ny=820, cmap='viridis', nmin=30, region='Whole', levels=np.arange(0,1.1,0.1), isvar=True):

    sic = ICE.sic
    
    fig = plt.figure(figsize=(wd,ht));
            
    ax = plt.axes(projection=ccrs.Orthographic(0, -90))
    #cbar_ax = fig.add_axes([0, 0, 0.1, 0.1])

    sic_contour = sic.isel(time=timeind).plot.contourf(ax=ax, transform=ccrs.PlateCarree(), levels=np.arange(0,1.1,0.1));
    ax.set_extent([lonmin, lonmax, latmin, latmax], crs=ccrs.PlateCarree())
    def resize_colobar(event):
        plt.draw()

        posn = ax.get_position()
        cbar_ax.set_position([posn.x0 + posn.width + 0.01, posn.y0,
                              0.04, posn.height])

    #fig.canvas.mpl_connect('resize_event', resize_colobar)
    #plt.colorbar(sic_contour, cax=cbar_ax)
    #ax.set_global();
    ax.coastlines();

        

    if(show==True):
        plt.show();
    if(save == True):
        plt.savefig(savename)
    

def find_area(df, m, units='Data Density',
                    nx=820, ny=820, lat0 = -90, lon0=0,region='Whole', years=[]):
    if not years:
        years = np.sort(df.JULD.dt.year.unique())
    area = np.zeros(len(years))
    
    for j in range(len(years)):
        year_mask = df['JULD'].dt.year == years[j]
        surfIndex = df[year_mask].groupby('PROFILE_NUMBER').head(1).index

        xlons, ylats = m(df[year_mask].loc[surfIndex , 'LONGITUDE'].values, df[year_mask].loc[surfIndex,'LATITUDE'].values)

        Xgrid = np.linspace(m.llcrnrx, m.urcrnrx, nx)
        Ygrid = np.linspace(m.llcrnry, m.urcrnry, ny)
        XX, YY = np.meshgrid(Xgrid, Ygrid)
        ni = np.zeros((len(Ygrid), len(Xgrid)))

        for i in range(len(xlons)):
            row = np.argmin(np.abs(Ygrid - ylats[i]))
            col = np.argmin(np.abs(Xgrid - xlons[i]))
            ni[row , col] = 1

        ni = ma.array(ni)
        ni_masked = ma.masked_equal(ni ,0)

        cell_size_x, cell_size_y = getCellSize(m, nx=nx, ny=ny)

        number_of_cells = np.sum(ni_masked)

        area[j] = cell_size_x * cell_size_y * number_of_cells
    
    return area

def plotDataDensity2layered(ax, df,mask1,mask2, units='Data Density',
                    save=False, savename="savedFig.png", wd=7, ht=7,
                    cx=10, cy=10, show=False, lat0 = -90, lon0=0, annotate=True,
                    levels=[1,5,10,20,30,50], region='Whole', m=None, fontsize=8):  #, 90, 150, 250, 500, 1000, 2500]):
                        
    matplotlib.rcParams.update({'font.size': fontsize})        
    #fig = plt.figure(figsize=(wd,ht));
    if(m == None):
        m  = createMapProjections(lat0, lon0, region=region, annotate=annotate)

    width = m.urcrnrx - m.llcrnrx
    height = m.urcrnry - m.llcrnry
    left = [m.llcrnrx, m.llcrnry+height*0.5]
    right = [m.urcrnrx, m.urcrnry-height*0.5]
    bottom = [m.llcrnrx+width*0.5, m.llcrnry]
    top = [m.llcrnrx+width*0.5, m.urcrnry]
    
    left_lonlat = m(left[0], left[1], inverse=True)
    right_lonlat = m(right[0], right[1], inverse=True)
    top_lonlat = m(top[0], top[1], inverse=True)
    bot_lonlat = m(bottom[0], bottom[1], inverse=True)
    dist_x = haversine(list(left_lonlat[::-1]), list(right_lonlat[::-1]))
    dist_y = haversine(top_lonlat[::-1], bot_lonlat[::-1])
    nx = dist_x // cx
    ny = dist_x // cy
    #lat_0 = -60, lon_0 = -20,

    surfIndex1 = df[mask1].groupby('PROFILE_NUMBER').head(1).index
    xlons1, ylats1 = m(df[mask1].loc[surfIndex1 , 'LONGITUDE'].values, df[mask1].loc[surfIndex1,'LATITUDE'].values)
    surfIndex2 = df[mask2].groupby('PROFILE_NUMBER').head(1).index
    xlons2, ylats2 = m(df[mask2].loc[surfIndex2 , 'LONGITUDE'].values, df[mask2].loc[surfIndex2,'LATITUDE'].values)

    Xgrid = np.linspace(m.llcrnrx, m.urcrnrx, nx)
    Ygrid = np.linspace(m.llcrnry, m.urcrnry, ny)
    XX, YY = np.meshgrid(Xgrid, Ygrid)
    ni1 = np.zeros((len(Ygrid), len(Xgrid)))
    ni2 = np.zeros((len(Ygrid), len(Xgrid)))
    
    for i in range(len(xlons1)):
        row = np.argmin(np.abs(Ygrid - ylats1[i]))
        col = np.argmin(np.abs(Xgrid - xlons1[i]))
        ni1[row , col] = ni1[row, col] + 1
        
    for i in range(len(xlons2)):
        row = np.argmin(np.abs(Ygrid - ylats2[i]))
        col = np.argmin(np.abs(Xgrid - xlons2[i]))
        ni2[row , col] = ni2[row, col] + 1

    ni1 = ma.array(ni1)
    ni1_masked = ma.masked_equal(ni1 ,0)
    ni2 = ma.array(ni2)
    ni2_masked = ma.masked_equal(ni2 ,0)

    #zi_masked = ma.masked_array(zi, 0)
    ticks = levels
    from matplotlib.colors import BoundaryNorm
    bnorm = BoundaryNorm(levels, ncolors=len(levels)-1, clip=False);
    
    cmap1 = LinearSegmentedColormap.from_list(name='linearCmap', colors=['mistyrose', 'salmon', 'darkred'], N=len(levels)-1); 
    cmap2 = LinearSegmentedColormap.from_list(name='linearCmap', colors=['paleturquoise', 'deepskyblue', 'navy'], N=len(levels)-1);
        
    CF1 = m.contourf(XX, YY, ni1_masked, vmin=min(levels), vmax=max(levels), levels=levels, norm=bnorm,
                    cmap=cmap1, extend='max', ax=ax); #, cmap='viridis'
    CF2 = m.contourf(XX, YY, ni2_masked, vmin=min(levels), vmax=max(levels), levels=levels, norm=bnorm,
                    cmap=cmap2, extend='max', ax=ax); #, cmap='viridis'

    ## cbar2 = m.colorbar(CF2, ticks=ticks, spacing='uniform', pad=0.6, location='bottom') #cax=cbaxes2) #boundaries=levels,
    ## cbar2.set_label('Area II, beyond 75km from GL')
    ## cbar1 = m.colorbar(CF1,  ticks=ticks, spacing='uniform', pad=0.6) #cax=cbaxes1) #boundaries=levels,    
    ## cbar1.set_label('Area I, within 75km of GL')

    ## plt.tight_layout();
    ## if(save==True):
    ##     plt.savefig(savename, dpi=300, bbox_inches='tight')
    ## if(show == True):
    ##     plt.show();
    ## else:
    ##     plt.close();
    return CF1, CF2

def createShadedRegionPlot(regionLonLims=[], regionNames=[], regionNamesxy=[], wd=7, ht=7, lat0=-90, lon0=0, region="Whole",
                           save=False, savename="untitled.png", fontsize=8, annotate=True, linewidth=0.2, dpi=300):
    matplotlib.rcParams.update({'font.size': fontsize, 'figure.autolayout': True})
    if not regionLonLims:
        print("Error! No regions mask provided")
        return 0
    fig = plt.figure(figsize=(wd,ht));
    #fig.set_facecolor('0.75')
    m  = createMapProjections(lat0, lon0, region=region, fontsize=fontsize, annotate=annotate, linewidth=0.2, draw_grounding_line=False)
    m.drawmapboundary(fill_color="grey", zorder=1, linewidth=0.1)
    m.drawcoastlines(linewidth=0.1, zorder=3)

    ## Fill the regions with white color, but only where depth is below 0
    for i in range(len(regionLonLims)):
        lons = np.linspace(regionLonLims[i][0], regionLonLims[i][1], 10)
        lats = np.linspace(-80, -60, 10)
        x,y = m(lons, lats)
        longrid, latgrid = np.meshgrid(lons, lats)
        white = np.zeros((len(lats), len(lons))) #np.random.randn(len(lats), len(lons) )
        m.contourf(longrid, latgrid, white, vmin=0, vmax=1, cmap="Greys", latlon=True, zorder=2) # cmap="Greys"
        
    m.fillcontinents(color='grey', zorder=2);
    
    for i in range(len(regionNames)):
        mxy = m(regionNamesxy[i][0], regionNamesxy[i][1])
        plt.annotate(regionNames[i], xy=(mxy[0], mxy[1]), fontsize=fontsize)
    plt.tight_layout()
    if(save == True):
        plt.savefig(savename, dpi=dpi, bbox_inches="tight")
    
    plt.show()


def plot_station_locations(positions, title=' ', save=False, savename="untitled.png", wd=12, ht=12, region='Whole', plotBathy=True, m=None):
    x = positions[:,1]
    y = positions[:,0]
    
    lat0 = -90
    lon0 = 0
    
    plt.figure(1, figsize=(wd,ht));
    gs = gridspec.GridSpec(2, 1, height_ratios=[1, 0.025])
    mapax = plt.subplot(gs[0, 0])
    if not m:
        m  = createMapProjections(lat0, lon0, region=region)

    m.scatter(x, y, latlon=True, color='r', s=1, marker='.')
    m.drawcoastlines()
    if(plotBathy == True):
        bathy = xr.open_dataset('/media/data/Datasets/Bathymetry/GEBCO_2014_2D.nc')
        lonlen = len(bathy.lon)
        lonindices = np.arange(0, lonlen+1, 30)
        lonindices[-1] = lonindices[-1] - 1
        bathyS = bathy.isel(lon=lonindices, lat=np.arange(0, 3600, 5))
        clevs = np.array([-100, -500, -1000, -1500, -2000, -3000])[::-1]
        
        longrid, latgrid = np.meshgrid(bathyS.lon.values, bathyS.lat.values)
        cs = m.contour(longrid, latgrid, bathyS.elevation.where(bathyS.elevation <= 0).values,  latlon=True, levels=clevs, linewidths=0.2, extend='min', ax=mapax) #, , cmap='rainbow'   , levels=clevs,
        ## plt.figure(2)
        ## cf = plt.contourf(longrid, latgrid,bathyS.elevation.where(bathyS.elevation <= 0).values, levels=clevs, extend='min') #, , cmap='rainbow'   , levels=clevs,
        ## plt.figure(1)
        bathycolorbar = plt.subplot(gs[1, 0])
        cbar1 = Colorbar(ax = bathycolorbar, mappable = cs, orientation = 'horizontal')
        cbar1.ax.get_children()[0].set_linewidths(5)
        cbar1.set_label('Depth (m)')

    if(save == True):
        plt.savefig(savename, dpi=300)#, bbox_inches='tight')
    plt.show()


def plotProfileNumberContours(df,var="PROFILE_NUMBER", plotBathy=True,
                         save=False, savename="savedFig.png",
                         wd=7, ht=7, cx=10, cy=10, cmap='viridis', nmin=0, region='Whole', levels=[0], show=True):
    plt.close(1)
    plt.figure(1,figsize=(wd,ht));
    gs = gridspec.GridSpec(2, 1, height_ratios=[1, 0.05])
    mapax = plt.subplot(gs[0,0])
        
    lat0 = -90
    lon0 = 0

    m  = createMapProjections(lat0, lon0, region=region)
    bathycolorbar = plt.subplot(gs[1, 0])
        
    width = m.urcrnrx - m.llcrnrx
    height = m.urcrnry - m.llcrnry
    left = [m.llcrnrx, m.llcrnry+height*0.5]
    right = [m.urcrnrx, m.urcrnry-height*0.5]
    bottom = [m.llcrnrx+width*0.5, m.llcrnry]
    top = [m.llcrnrx+width*0.5, m.urcrnry]
    
    left_lonlat = m(left[0], left[1], inverse=True)
    right_lonlat = m(right[0], right[1], inverse=True)
    top_lonlat = m(top[0], top[1], inverse=True)
    bot_lonlat = m(bottom[0], bottom[1], inverse=True)
    dist_x = haversine(list(left_lonlat[::-1]), list(right_lonlat[::-1]))
    dist_y = haversine(top_lonlat[::-1], bot_lonlat[::-1])
    nx = dist_x // cx
    ny = dist_x // cy
    print(nx, ny)
    surfIndex = df.groupby('PROFILE_NUMBER').head(1).index
    
    z = df.loc[surfIndex, var].values
    xlons, ylats = m(df.loc[surfIndex , 'LONGITUDE'].values, df.loc[surfIndex,'LATITUDE'].values)

    if(plotBathy == True):
        bathy = xr.open_dataset('/media/data/Datasets/Bathymetry/GEBCO_2014_2D.nc')
        lonlen = len(bathy.lon)
        lonindices = np.arange(0, lonlen+1, 30)
        lonindices[-1] = lonindices[-1] - 1
        bathyS = bathy.isel(lon=lonindices, lat=np.arange(0, 3600, 5))
        clevs = np.array([-100, -500, -1000, -1500, -2000, -3000])[::-1]
        
        longrid, latgrid = np.meshgrid(bathyS.lon.values, bathyS.lat.values)
        cs = m.contour(longrid, latgrid, bathyS.elevation.where(bathyS.elevation <= 0).values,  latlon=True, levels=clevs, linewidths=0.2, extend='min', ax=mapax) #, , cmap='rainbow'   , levels=clevs,
        ## plt.figure(2)
        ## cf = plt.contourf(longrid, latgrid,bathyS.elevation.where(bathyS.elevation <= 0).values, levels=clevs, extend='min') #, , cmap='rainbow'   , levels=clevs,
        ## plt.figure(1)
        cbar1 = Colorbar(ax = bathycolorbar, mappable = cs, orientation = 'horizontal')
        cbar1.ax.get_children()[0].set_linewidths(5)
        cbar1.set_label('Depth (m)')

    
    if(save==True):
        plt.savefig(savename, dpi=300)
    if(show == True):
        plt.show();
    else:
        plt.close();


def plot_scalar_field(sclr, lons, lats, m=None, bathy=None, plotBathy=True, save=False, savename="untitled.png", levs=None, drawMeridians=True, meridians=None, pickPoints=True):
    plt.close(1)
    fig = plt.figure(1, figsize=(10,10))
    gs = gridspec.GridSpec(2, 2, height_ratios=[1, 0.05], width_ratios=[1, 0.05])
    mapax = plt.subplot(gs[0,0])
    colorbar_ax = plt.subplot(gs[1,0])
    bathycolorbar = plt.subplot(gs[0, 1])
    
    X, Y = np.meshgrid(lons, lats)
    
    if not m:
        m = createMapProjections(-90, 0, region="Whole")
        
    CF = m.contourf(X, Y, sclr, ax=mapax, latlon=True, levels=levs, extend="both")
    cbar1 = Colorbar(ax = colorbar_ax, mappable = CF, orientation = 'horizontal')
    
    if(plotBathy == True):
        if not bathy:
            bathy = xr.open_dataset('/media/data/Datasets/Bathymetry/GEBCO_2014_2D.nc')
        lonlen = len(bathy.lon)
        lonindices = np.arange(0, lonlen+1, 30)
        lonindices[-1] = lonindices[-1] - 1
        bathyS = bathy.isel(lon=lonindices, lat=np.arange(0, 3600, 5))
        clevs = np.array([-100, -500, -1000, -1500, -2000, -3000])[::-1]
        
        longrid, latgrid = np.meshgrid(bathyS.lon.values, bathyS.lat.values)
        cs = m.contour(longrid, latgrid, bathyS.elevation.where(bathyS.elevation <= 0).values,  latlon=True, 
                       levels=clevs, linewidths=0.2, extend='min', ax=mapax, colors='k' ) #, cmap="RdYlBu"
        cbar1 = Colorbar(ax = bathycolorbar, mappable = cs, orientation = 'vertical')
        cbar1.ax.get_children()[0].set_linewidths(5)
        cbar1.set_label('Depth (m)')
        
    if drawMeridians:
        if not meridians:
            meridians = np.arange(-180, 180, 20)
        m.drawmeridians(meridians,labels=[1,1,1,1], linewidth=0.2, ax=mapax)
    if save:
        plt.savefig(savename, dpi=600)
    m.drawcoastlines(ax=mapax)
    plt.show()    

    if pickPoints:
        cid = fig.canvas.mpl_connect('key_press_event', lambda event: onpick3(event, m))

        
def onpick3(event, m):
    print( m(event.xdata, event.ydata, inverse=True))

    
def plot_region_bathy(region=None, wd=7, ht=5, bathy=None, save=False, savename="untitled.png", levs=[]):
    plt.close(1)
    fig = plt.figure(1, figsize=(wd, ht))
    gs = gridspec.GridSpec(4, 2, width_ratios=[1, 0.02])

    if(region == "CDP"):
        llcrnrlon, llcrnrlat = 60, -68
        urcrnrlon, urcrnrlat = 71, -65.5
        map_proj = ccrs.PlateCarree() #ccrs.Orthographic(65, -90)
    if(region == "WPB"):
        llcrnrlon, llcrnrlat = 69, -70
        urcrnrlon, urcrnrlat = 80, -66
        map_proj = ccrs.PlateCarree() #ccrs.Orthographic(65, -90)
    if(region == "EPB"):
        llcrnrlon, llcrnrlat = 69, -70
        urcrnrlon, urcrnrlat = 83, -64.5
        map_proj = ccrs.PlateCarree() #ccrs.Orthographic(65, -90)
    if(region == "LAC"):
        llcrnrlon, llcrnrlat = 80, -68
        urcrnrlon, urcrnrlat = 89, -64.5
        map_proj = ccrs.PlateCarree() #ccrs.Orthographic(65, -90)
    if(region == "KC"):
        llcrnrlon, llcrnrlat = 99, -68
        urcrnrlon, urcrnrlat = 112, -63
        map_proj = ccrs.PlateCarree() #ccrs.Orthographic(65, -90)
    if(region == "AC"):
        llcrnrlon, llcrnrlat = 133, -68
        urcrnrlon, urcrnrlat = 147, -64
        map_proj = ccrs.PlateCarree() #ccrs.Orthographic(65, -90)
    if(region == "RS"):
        llcrnrlon, llcrnrlat = 160, -79
        urcrnrlon, urcrnrlat = 180, -71
        map_proj = ccrs.PlateCarree() #ccrs.Orthographic(65, -90)
    if(region == "AS"):
        llcrnrlon, llcrnrlat = -122, -75.5
        urcrnrlon, urcrnrlat = -98, -70
        map_proj = ccrs.PlateCarree() #ccrs.Orthographic(65, -90)
    if(region == "BS"):
        llcrnrlon, llcrnrlat = -102, -71.5
        urcrnrlon, urcrnrlat = -58, -60
        map_proj = ccrs.PlateCarree() #ccrs.Orthographic(65, -90)        
    if(region == "WS"):
        llcrnrlon, llcrnrlat = -50, -82
        urcrnrlon, urcrnrlat = -20, -72
        map_proj = ccrs.PlateCarree() #ccrs.Orthographic(65, -90)
    if(region == "PMC"):
        llcrnrlon, llcrnrlat = -19, -74
        urcrnrlon, urcrnrlat = 0, -65
        map_proj = ccrs.PlateCarree() #ccrs.Orthographic(65, -90)
    if(region == "PHC"):
        llcrnrlon, llcrnrlat = 28, -71
        urcrnrlon, urcrnrlat = 38, -65
        map_proj = ccrs.PlateCarree() #ccrs.Orthographic(65, -90)                
    else:
        map_proj = ccrs.PlateCarree()
    
    
    mapax = plt.subplot(gs[:,0], projection=map_proj)
    mapax.set_extent([llcrnrlon, urcrnrlon, llcrnrlat, urcrnrlat], crs=ccrs.PlateCarree())
    colorbar_ax = plt.subplot(gs[1:3,1])

    try:
        if bathy.variables:
            pass
    except:
        bathy = xr.open_dataset(MEOPDIR+'/Datasets/Bathymetry/GEBCO_2014_2D.nc')
    bathyS = bathy.sel(lon=slice(llcrnrlon-0.5, urcrnrlon+0.5), lat=slice(llcrnrlat-0.5, urcrnrlat+0.5))

    cs2 = mapax.contourf(bathyS.lon, bathyS.lat, bathyS.elevation.where(bathyS.elevation <= 0).values,  levels=levs, extend="min", transform=ccrs.PlateCarree()) #, , cmap='rainbow'   , levels=clevs,
    ## plt.figure(2)
    ## cf = plt.contourf(longrid, latgrid,bathyS.elevation.where(bathyS.elevation <= 0).values, levels=clevs, extend='min') #, , cmap='rainbow'   , levels=clevs,
    ## plt.figure(1)
    cbar = Colorbar(ax = colorbar_ax, mappable = cs2, orientation = 'vertical')
    cbar.ax.get_children()[0].set_linewidths(5)
    cbar.set_label('Depth (m)')
    MEOPDIR = "/media/data"
    shpfile = MEOPDIR+"/Datasets/Shapefiles/AntarcticGroundingLine/GSHHS_f_L6.shp"
    with fiona.open(shpfile) as records:
        geometries = [sgeom.shape(shp['geometry']) for shp in records]
    mapax.add_geometries(geometries, ccrs.PlateCarree(), edgecolor='gray', facecolor='none')
    
    ISedgefname = MEOPDIR+"/Datasets/Shapefiles/AntIceShelfEdge/ne_10m_antarctic_ice_shelves_lines.shp"
    ISe_feature = ShapelyFeature(Reader(ISedgefname).geometries(), ccrs.PlateCarree(), facecolor='none', edgecolor="k")
    mapax.add_feature(ISe_feature, zorder=3)

    mapax.set_extent([llcrnrlon, urcrnrlon, llcrnrlat, urcrnrlat])        
    gl = mapax.gridlines(crs=ccrs.PlateCarree(), draw_labels=True, linewidth=1, color='gray', alpha=0.5, linestyle='--')
    gl.ylabels_right = False
    gl.xlabels_top = False
    gl.xformatter = LONGITUDE_FORMATTER
    gl.yformatter = LATITUDE_FORMATTER
    gl.xlabel_style = {'size': 8, 'color': 'k'}
    gl.ylabel_style = {'size': 8, 'color': 'k'}

    if save:
        plt.savefig(savename, dpi=300)

    plt.show()



def plotDataDensity_NIS_DIS(df1, df2, units='Data Density', save=False, savename="savedFig.png", wd=7, ht=7, show=True, mapax = None, subplotlabel=None, levels=[0, 10, 20, 30, 40, 50, 60, 100, 200, 500], region='Whole', plotBathy=True, fontsize=8, DATADIR = "/media/sda7", dx=0.5, dy=0.5, region_lons=None, xticks = np.arange(-180, 181, 20), yticks = np.arange(-80, -59, 5), xlabels_top=False, xlabels_bottom=False, ylabels_left=False, ylabels_right=False, draw_lonlat_label=True):

    df = [df1, df2]    
    matplotlib.rcParams.update({'font.size': fontsize})        # setting fontsize for plot elements            
    plt.figure(1, figsize=(wd,ht));
    external_mapaxis = True
    if not mapax:
        print("creating local mapaxis")
        external_mapaxis = False
        gs = gridspec.GridSpec(1, 3, height_ratios=[1], width_ratios=[1, 0.03, 0.03])
        mapax = plt.subplot(gs[0, 0], projection = ccrs.Orthographic(central_longitude=0, central_latitude=-90) )
        datacolorbar1 = plt.subplot(gs[0, 1])
        #datacolorbar2 = plt.subplot(gs[0, 2])

    shpfile = DATADIR+"/Datasets/Shapefiles/AntarcticGroundingLine/GSHHS_f_L6.shp"
    with fiona.open(shpfile) as records:
        geometries = [sgeom.shape(shp['geometry']) for shp in records]
    ISedgefname = DATADIR+"/Datasets/Shapefiles/AntIceShelf/ne_10m_antarctic_ice_shelves_polys.shp"
    ISe_feature = ShapelyFeature(Reader(ISedgefname).geometries(), 
                                 ccrs.PlateCarree(), linewidth=0.2,
                                 facecolor='none', 
                                 edgecolor="k")

    mapax.add_geometries(geometries, ccrs.PlateCarree(), edgecolor='0.25', facecolor='0.7',alpha=0.25, linewidth=0.2)
    mapax.add_feature(ISe_feature, zorder=3)


    gl = mapax.gridlines(crs=ccrs.PlateCarree(), draw_labels=False, zorder=2,
                  linewidth=0.5, color='gray', alpha=1, linestyle='--')
    
    gl.xlocator = mticker.FixedLocator(xticks) 
    gl.ylocator = mticker.FixedLocator(yticks) 
    gl.xformatter = LONGITUDE_FORMATTER
    gl.yformatter = LATITUDE_FORMATTER
    gl.xlabels_top = xlabels_top
    gl.xlabels_bottom = xlabels_bottom
    gl.ylabels_left = ylabels_left
    gl.ylabels_right = ylabels_right

    if region_lons:
        gl_regions = mapax.gridlines(crs=ccrs.PlateCarree(), draw_labels=False, linewidth=0.5, color='blue', alpha=1, linestyle=':', zorder=3)
        #gl_regions_lowZ = mapax.gridlines(crs=ccrs.PlateCarree(), draw_labels=False, linewidth=0.5, color='blue', alpha=1, linestyle='--', zorder=1)
        gl_regions.xlocator = mticker.FixedLocator(region_lons)
        gl_regions.ylocator = mticker.FixedLocator(np.arange(-80, -59, 5) )
        gl_regions.xformatter = LONGITUDE_FORMATTER
        gl_regions.yformatter = LATITUDE_FORMATTER
        gl_regions.ylines = False
    
    matplotlib.rcParams.update({'font.size': fontsize})    
    
    color_scheme = [['paleturquoise', 'deepskyblue', 'navy'], ['mistyrose', 'salmon', 'darkred']]
    #color_scheme = ['paleturquoise', 'deepskyblue', 'navy']
    CF = []
    for r in range(1):
        surfIndex = df[r].groupby('PROFILE_NUMBER').head(1).index

        xlons, ylats = df[r].loc[surfIndex , 'LONGITUDE'].values, df[r].loc[surfIndex,'LATITUDE'].values

        Xgrid = np.arange(-180, 180.1, dx)
        Ygrid = np.arange(-80, -60.1, dy)
        XX, YY = np.meshgrid(Xgrid, Ygrid)
        ni = np.zeros((len(Ygrid), len(Xgrid)))


        for i in range(len(xlons)):
            row = np.argmin(np.abs(Ygrid - ylats[i]))
            col = np.argmin(np.abs(Xgrid - xlons[i]))
            ni[row , col] = ni[row, col] + 1

        ni = ma.array(ni)
        ni_masked = ma.masked_equal(ni ,0)

        #zi_masked = ma.masked_array(zi, 0)
        ticks = levels
        from matplotlib.colors import BoundaryNorm
        bnorm = BoundaryNorm(levels, ncolors=len(levels)-1, clip=False)

        cmap = LinearSegmentedColormap.from_list(name='linearCmap', colors=color_scheme[1], N=len(levels)-1) 

        CF.append(mapax.pcolormesh(XX, YY, ni_masked, vmin=min(levels), vmax=max(levels), norm=bnorm, 
                                   cmap=cmap, transform = ccrs.PlateCarree(), zorder=1))
    
    
    if(plotBathy == True):
        bathyS = xr.open_dataset('/media/hdd2/SOSE_1_12/bathyS.nc')
        cs = mapax.contour(bathyS.lon, bathyS.lat, bathyS.elevation.where(bathyS.elevation <= 0).values,  levels=[-3000], colors="b", linestyle=":", linewidths=0.25, transform = ccrs.PlateCarree())
        
    if not external_mapaxis:
        print("no external mapaxis")
        cbar1 = Colorbar(ax = datacolorbar1, mappable = CF[0], ticks=ticks, orientation='vertical') 
    #cbar2 = Colorbar(ax = datacolorbar2, mappable = CF[1], ticks=ticks, orientation='vertical')
    #cbar1.ax.set_yticklabels("")
    #cbar2.set_label("Data density")


    if draw_lonlat_label:
        for l in np.arange(0, 360, 20):
            if( (l == 80) or (l == 100) ):
                text_lat = -80
            else:
                text_lat  = -62.5
            mapax.text(l, text_lat, str(l)+"$^{\circ}$", transform=ccrs.PlateCarree() )
    if not subplotlabel:
        pass
    else:
        mapax.text(-134, -49, subplotlabel, transform = ccrs.PlateCarree() )
    if(save==True):
        plt.savefig(savename, dpi=300, bbox_inches='tight')
    if mapax:
        return CF[0]
    else:
        plt.show();

        

def plot_fields_orthographic(field, longitude_coord, latitude_coord, vmin, vmax, cmap = 'viridis', save=False, savename="untitled.png", norm=None,text_lon=65,
                             colorbar_label="", region_lons=None, mapax=None, plot_colorbar=True, xlim=None, ylim=None, xticks=np.arange(-180, 181, 30), yticks=np.arange(-80, -59, 5), xlabels_top=False, xlabels_bottom=False, ylabels_left=False, ylabels_right=False, plot_latlon_labels=True, cl=0):
    DATADIR = "/data/SO12/sio-anarayanan/shapefiles/"
    passing_mapaxis = True
    if not mapax:
        passing_mapaxis = False
        plt.close(1)
        plt.figure(1, figsize=(7,7) )
        gs = gridspec.GridSpec(5, 2, width_ratios=[1, 0.05], wspace=0.05)
    
        mapax = plt.subplot(gs[:,0], projection = ccrs.Orthographic(central_longitude=cl, central_latitude=-90) )
        


    shpfile = DATADIR+"AntarcticGroundingLine/GSHHS_f_L6.shp"
    # with fiona.open(shpfile) as records:
    #     geometries = [sgeom.shape(shp['geometry']) for shp in records]
    GLfeature = ShapelyFeature(Reader(shpfile).geometries(), ccrs.PlateCarree(), linewidth=0.2, facecolor='0.5', edgecolor='0.25', zorder=10)
    ISedgefname = "/data/SO12/sio-anarayanan/shapefiles/AntIceShelfEdge/ne_10m_antarctic_ice_shelves_lines.shp" #DATADIR+"AntIceShelfEdge/ne_10m_antarctic_ice_shelves_polys.shp"
    ISe_feature = ShapelyFeature(Reader(ISedgefname).geometries(), 
                                 ccrs.PlateCarree(), linewidth=0.2, zorder=6,
                                 facecolor='azure', 
                                 edgecolor="k")

    #mapax.add_geometries(geometries, ccrs.PlateCarree(), edgecolor='0.25', facecolor='0.7',alpha=0.25, linewidth=0.2)
    mapax.add_feature(GLfeature, zorder=10)
    mapax.add_feature(ISe_feature, zorder=10)
    
    CF = mapax.pcolormesh(longitude_coord, latitude_coord, 
                          field, transform=ccrs.PlateCarree(), cmap= cmap, vmin = vmin, vmax = vmax, zorder=1, norm=norm)

    if plot_latlon_labels:
        for l in yticks[:]:
            text_lon = text_lon
            mapax.text(text_lon, l, str(l)+"$^{\circ}$S", transform=ccrs.PlateCarree(), zorder=4, rotation=90 )

    if not passing_mapaxis:
        colorbar_ax = plt.subplot(gs[1:-1, 1])
        Colorbar(mappable = CF, ax = colorbar_ax)
        colorbar_ax.set_ylabel(colorbar_label)

    gl = mapax.gridlines(crs=ccrs.PlateCarree(central_longitude=cl), draw_labels=True, zorder=4,
                  linewidth=0.5, color='gray', alpha=1, linestyle='--')
    
    
    if region_lons:
        gl_regions = mapax.gridlines(crs=ccrs.PlateCarree(central_longitude=cl), draw_labels=False, linewidth=0.5, color='blue', alpha=1, linestyle=':', zorder=3)
        #gl_regions_lowZ = mapax.gridlines(crs=ccrs.PlateCarree(), draw_labels=False, linewidth=0.5, color='blue', alpha=1, linestyle='--', zorder=1)
        gl_regions.xlocator = mticker.FixedLocator(region_lons)
        gl_regions.ylocator = mticker.FixedLocator(np.arange(-80, -59, 5) )
        gl_regions.xformatter = LONGITUDE_FORMATTER
        gl_regions.yformatter = LATITUDE_FORMATTER
        gl_regions.ylines = False

    gl.xlocator = mticker.FixedLocator(xticks)
    gl.ylocator = mticker.FixedLocator(yticks)
    gl.xformatter = LONGITUDE_FORMATTER
    gl.yformatter = LATITUDE_FORMATTER
    
    
    gl.top_labels = xlabels_top
    gl.bottom_labels = xlabels_bottom
    gl.left_labels = ylabels_left
    gl.right_labels = ylabels_right
    
    if save:
        plt.savefig(savename, dpi=150, bbox_inches="tight")

    if passing_mapaxis:
        return CF
    plt.show()



    
def plot_fields_polarStereographic(field, longitude_coord, latitude_coord, vmin, vmax, cmap = 'viridis', save=False, savename="untitled.png",
                             colorbar_label="", region_lons=None, mapax=None, plot_colorbar=True, xlim=None, ylim=None):
    DATADIR = "/media/aditya/AdityaNarayanan"
    passing_mapaxis = True
    if not mapax:
        passing_mapaxis = False
        plt.close(1)
        plt.figure(1, figsize=(7,7) )
        gs = gridspec.GridSpec(5, 2, width_ratios=[1, 0.05], wspace=0.05)
    
        mapax = plt.subplot(gs[:,0], projection = ccrs.SouthPolarStereo(central_longitude=0, globe=None) )
    print("plotting pcolormesh")
    CF = mapax.pcolormesh(longitude_coord, latitude_coord, 
                          field,  cmap= cmap, vmin = vmin, vmax = vmax)
    print("completed pcolormesh")
    gl = mapax.gridlines(crs=ccrs.PlateCarree(), draw_labels=False, zorder=2,
                  linewidth=0.5, color='gray', alpha=1, linestyle='--')
    
    gl.xlocator = mticker.FixedLocator(np.arange(-180, 181, 20))
    gl.ylocator = mticker.FixedLocator(np.arange(-80, -59, 5))
    gl.xformatter = LONGITUDE_FORMATTER
    gl.yformatter = LATITUDE_FORMATTER
    print("completed gridlines")
    if region_lons:
        gl_regions = mapax.gridlines(crs=ccrs.PlateCarree(), draw_labels=False, linewidth=0.5, color='blue', alpha=1, linestyle=':', zorder=3)
        #gl_regions_lowZ = mapax.gridlines(crs=ccrs.PlateCarree(), draw_labels=False, linewidth=0.5, color='blue', alpha=1, linestyle='--', zorder=1)
        gl_regions.xlocator = mticker.FixedLocator(region_lons)
        gl_regions.ylocator = mticker.FixedLocator(np.arange(-80, -59, 5) )
        gl_regions.xformatter = LONGITUDE_FORMATTER
        gl_regions.yformatter = LATITUDE_FORMATTER
        gl_regions.ylines = False
    print("plotting shape files")
    shpfile = DATADIR+"/Datasets/Shapefiles/AntarcticGroundingLine/GSHHS_f_L6.shp"
    # with fiona.open(shpfile) as records:
    #     geometries = [sgeom.shape(shp['geometry']) for shp in records]
    GLfeature = ShapelyFeature(Reader(shpfile).geometries(), ccrs.PlateCarree(), linewidth=0.2, facecolor='0.7', edgecolor='0.25')
    ISedgefname = DATADIR+"/Datasets/Shapefiles/AntIceShelf/ne_10m_antarctic_ice_shelves_polys.shp"
    ISe_feature = ShapelyFeature(Reader(ISedgefname).geometries(), 
                                 ccrs.PlateCarree(), linewidth=0.2,
                                 facecolor='none', 
                                 edgecolor="k")
    print("completed shape files")
    #mapax.add_geometries(geometries, ccrs.PlateCarree(), edgecolor='0.25', facecolor='0.7',alpha=0.25, linewidth=0.2)
    mapax.add_feature(GLfeature, zorder=3)
    mapax.add_feature(ISe_feature, zorder=3)

    # for l in np.arange(0, 360, 20):
    #     if( (l == 80) or (l == 100) ):
    #         text_lat = -62.5
    #     else:
    #         text_lat  = -62.5
    #     mapax.text(l, text_lat, str(l)+"$^{\circ}$", transform=ccrs.PlateCarree() )    

    for l in np.arange(0, 360, 20):
        if( (l == 80) or (l == 100) ):
            text_lat = -80
        else:
            text_lat  = -62.5
        mapax.text(l, text_lat, str(l)+"$^{\circ}$", transform=ccrs.PlateCarree() )
    
    if not passing_mapaxis:
        colorbar_ax = plt.subplot(gs[1:-1, 1])
        Colorbar(mappable = CF, ax = colorbar_ax)
        colorbar_ax.set_ylabel(colorbar_label)
    
    if save:
        plt.savefig(savefig, dpi=600)

    if passing_mapaxis:
        return CF
    plt.show()




def plot_scatter_orthographic(field, longitude_coord, latitude_coord, vmin, vmax, cmap = 'viridis', save=False, savename="untitled.png",
                             colorbar_label="", region_lons=None, mapax=None, plot_colorbar=True, xlim=None, ylim=None, xticks=np.arange(-180, 181, 20), yticks=np.arange(-80, -59, 5), xlabels_top=False, xlabels_bottom=False, ylabels_left=False, ylabels_right=False):
    DATADIR = "/media/sda7/"
    passing_mapaxis = True
    if not mapax:
        passing_mapaxis = False
        plt.close(1)
        plt.figure(1, figsize=(7,7) )
        gs = gridspec.GridSpec(5, 2, width_ratios=[1, 0.05], wspace=0.05)
    
        mapax = plt.subplot(gs[:,0], projection = ccrs.Orthographic(central_longitude=0, central_latitude=-90) )
        

    gl = mapax.gridlines(crs=ccrs.PlateCarree(), draw_labels=False, zorder=2,
                  linewidth=0.5, color='gray', alpha=1, linestyle='--')
    
    gl.xlocator = mticker.FixedLocator(xticks)
    gl.ylocator = mticker.FixedLocator(yticks)
    gl.xformatter = LONGITUDE_FORMATTER
    gl.yformatter = LATITUDE_FORMATTER
    gl.xlabels_top = xlabels_top
    gl.xlabels_bottom = xlabels_bottom
    gl.ylabels_left = ylabels_left
    gl.ylabels_right = ylabels_right
    
    if region_lons:
        gl_regions = mapax.gridlines(crs=ccrs.PlateCarree(), draw_labels=False, linewidth=0.5, color='blue', alpha=1, linestyle=':', zorder=3)
        #gl_regions_lowZ = mapax.gridlines(crs=ccrs.PlateCarree(), draw_labels=False, linewidth=0.5, color='blue', alpha=1, linestyle='--', zorder=1)
        gl_regions.xlocator = mticker.FixedLocator(region_lons)
        gl_regions.ylocator = mticker.FixedLocator(np.arange(-80, -59, 5) )
        gl_regions.xformatter = LONGITUDE_FORMATTER
        gl_regions.yformatter = LATITUDE_FORMATTER
        gl_regions.ylines = False

    shpfile = DATADIR+"/Datasets/Shapefiles/AntarcticGroundingLine/GSHHS_f_L6.shp"
    # with fiona.open(shpfile) as records:
    #     geometries = [sgeom.shape(shp['geometry']) for shp in records]
    GLfeature = ShapelyFeature(Reader(shpfile).geometries(), ccrs.PlateCarree(), linewidth=0.2, facecolor='0.9', edgecolor='0.25', zorder=3)
    ISedgefname = DATADIR+"/Datasets/Shapefiles/AntIceShelf/ne_10m_antarctic_ice_shelves_polys.shp"
    ISe_feature = ShapelyFeature(Reader(ISedgefname).geometries(), 
                                 ccrs.PlateCarree(), linewidth=0.2, zorder=3,
                                 facecolor='azure', 
                                 edgecolor="k")

    #mapax.add_geometries(geometries, ccrs.PlateCarree(), edgecolor='0.25', facecolor='0.7',alpha=0.25, linewidth=0.2)
    mapax.add_feature(GLfeature, zorder=3)
    mapax.add_feature(ISe_feature, zorder=3)
    
    CF = mapax.scatter(longitude_coord, latitude_coord, 
                          c=field, transform=ccrs.PlateCarree(), cmap= cmap, vmin = vmin, vmax = vmax, zorder=1)

    # for l in np.arange(0, 360, 20):
    #     if( (l == 80) or (l == 100) ):
    #         text_lat = -62.5
    #     else:
    #         text_lat  = -62.5
    #     mapax.text(l, text_lat, str(l)+"$^{\circ}$", transform=ccrs.PlateCarree() )    

    if not passing_mapaxis:
        colorbar_ax = plt.subplot(gs[1:-1, 1])
        Colorbar(mappable = CF, ax = colorbar_ax)
        colorbar_ax.set_ylabel(colorbar_label)
    
    if save:
        plt.savefig(savefig, dpi=600)

    if passing_mapaxis:
        return CF
    plt.show()
