diff --git a/cool_maps/calc.py b/cool_maps/calc.py
index 3d69ed0..4b6ee36 100644
--- a/cool_maps/calc.py
+++ b/cool_maps/calc.py
@@ -3,13 +3,14 @@
 import numpy as np
 
 
-def calculate_ticks(extent, direction):
+def calculate_ticks(extent, direction, decimal_degrees=False):
     """
     Define major and minor tick locations and major tick labels
 
     Args:
         extent (tuple or list): extent (x0, x1, y0, y1) of the map in the given coordinate system.
         dirs (str): Tell function which bounds to calculate. Must be 'longitude' or 'latitude'.
+        decimal_degrees (bool, optional): label ticks in decimal degrees (True) or degree-minute-second (False, default)
 
     Returns:
         list: minor ticks
@@ -94,40 +95,91 @@ def calculate_ticks(extent, direction):
         major_int)
     major_ticks = major_ticks[major_ticks <= l1]
 
-    if major_int < 1:
-        d, m, _ = dd2dms(np.array(major_ticks))
-        if direction == 'longitude':
-            n = 'W' * sum(d < 0)
-            p = 'E' * sum(d >= 0)
-            dir = n + p
-            major_tick_labels = [str(np.abs(int(d[i]))) + u"\N{DEGREE SIGN}" + str(int(m[i])) + "'" + dir[i] for i in
-                                 range(len(d))]
-        elif direction == 'latitude':
-            n = 'S' * sum(d < 0)
-            p = 'N' * sum(d >= 0)
-            dir = n + p
-            major_tick_labels = [str(np.abs(int(d[i]))) + u"\N{DEGREE SIGN}" + str(int(m[i])) + "'" + dir[i] for i in
-                                 range(len(d))]
-        else:
-            major_tick_labels = [str(int(d[i])) + u"\N{DEGREE SIGN}" + str(int(m[i])) + "'" for i in range(len(d))]
+    if decimal_degrees:
+        major_tick_labels = major_ticks
     else:
-        d = major_ticks
-        if direction == 'longitude':
-            n = 'W' * sum(d < 0)
-            p = 'E' * sum(d >= 0)
-            dir = n + p
-            major_tick_labels = [str(np.abs(int(d[i]))) + u"\N{DEGREE SIGN}" + dir[i] for i in range(len(d))]
-        elif direction == 'latitude':
-            n = 'S' * sum(d < 0)
-            p = 'N' * sum(d >= 0)
-            dir = n + p
-            major_tick_labels = [str(np.abs(int(d[i]))) + u"\N{DEGREE SIGN}" + dir[i] for i in range(len(d))]
+        if major_int < 1:
+            d, m, _ = dd2dms(np.array(major_ticks))
+            if direction == 'longitude':
+                n = 'W' * sum(d < 0)
+                p = 'E' * sum(d >= 0)
+                dir = n + p
+                major_tick_labels = [str(np.abs(int(d[i]))) + u"\N{DEGREE SIGN}" + str(int(m[i])) + "'" + dir[i] for i in
+                                    range(len(d))]
+            elif direction == 'latitude':
+                n = 'S' * sum(d < 0)
+                p = 'N' * sum(d >= 0)
+                dir = n + p
+                major_tick_labels = [str(np.abs(int(d[i]))) + u"\N{DEGREE SIGN}" + str(int(m[i])) + "'" + dir[i] for i in
+                                    range(len(d))]
+            else:
+                major_tick_labels = [str(int(d[i])) + u"\N{DEGREE SIGN}" + str(int(m[i])) + "'" for i in range(len(d))]
         else:
-            major_tick_labels = [str(int(d[i])) + u"\N{DEGREE SIGN}" for i in range(len(d))]
+            d = major_ticks
+            if direction == 'longitude':
+                n = 'W' * sum(d < 0)
+                p = 'E' * sum(d >= 0)
+                dir = n + p
+                major_tick_labels = [str(np.abs(int(d[i]))) + u"\N{DEGREE SIGN}" + dir[i] for i in range(len(d))]
+            elif direction == 'latitude':
+                n = 'S' * sum(d < 0)
+                p = 'N' * sum(d >= 0)
+                dir = n + p
+                major_tick_labels = [str(np.abs(int(d[i]))) + u"\N{DEGREE SIGN}" + dir[i] for i in range(len(d))]
+            else:
+                major_tick_labels = [str(int(d[i])) + u"\N{DEGREE SIGN}" for i in range(len(d))]
 
     return minor_ticks, major_ticks, major_tick_labels
 
 
+def calculate_colorbar_ticks(vmin, vmax, c0=False):
+    """
+    Calculate tick locations for colorbar
+
+    Args:
+        vmin (float): minimum value of colorbar (should match vmin of object colorbar is mapped to)
+        vmax (float): maximum value of colorbar (should match vmax of object colorbar is mapped to)
+        c0 (bool): center values around 0 (for anomaly products)
+
+    Returns:
+        list: tick locations
+    """
+
+    if c0:
+        vmax=np.max((np.abs(vmin), np.abs(vmax)))
+        vmin=0
+        scale = 1
+        if vmax-vmin<0:
+            scale = 10**(-np.floor(np.log10(vmax-vmin))+1)
+
+        cbticks = np.arange(vmin, np.floor(vmax*scale+.99))
+        if len(cbticks)>3:
+            cbticks=cbticks[::int(np.floor(len(cbticks)/3))]
+        cbticks = cbticks/scale
+        i=np.diff(cbticks)[0]
+        if cbticks[-1]+i<=vmax:
+            cbticks=np.append(cbticks, cbticks[-1]+i)
+        i=np.diff(cbticks)[0]
+        cbticks=np.append(np.arange(-np.max(cbticks),0,i),cbticks)
+        return cbticks
+    
+    scale = 1
+    if vmax-vmin<4:
+        scale = 10**(-np.floor(np.log10(vmax-vmin))+1)
+
+    cbticks = np.arange(np.ceil(vmin*scale+.01), np.floor(vmax*scale+.99))
+    if len(cbticks)>5:
+        cbticks=cbticks[::int(np.floor(len(cbticks)/5))]
+    cbticks = cbticks/scale
+    i=np.diff(cbticks)[0]
+    if cbticks[0]-i>=vmin:
+        cbticks=np.append(cbticks[0]-i, cbticks)
+    if cbticks[-1]+i<=vmax:
+        cbticks=np.append(cbticks, cbticks[-1]+i)
+    
+    return cbticks
+
+
 def categorical_cmap(nc, nsc, cmap="tab10", continuous=False):
     """
     Expand your colormap by changing the alpha value (opacity) of each color.
diff --git a/cool_maps/colormaps.py b/cool_maps/colormaps.py
new file mode 100644
index 0000000..eb5a01f
--- /dev/null
+++ b/cool_maps/colormaps.py
@@ -0,0 +1,389 @@
+import cmocean
+import matplotlib.pyplot as plt
+import matplotlib.colors as mcolors
+import numpy as np
+
+"""
+Modified colormaps: colormaps with defined breakpoints; categorical colormaps
+
+Functions to modify colormaps with defined breakpoints.
+Mostly for oxygen, but who knows where this will take us.
+Individual colors can be removed from any version
+
+Based on our best research (Fed and NJ) so far, 
+and concentration/saturation equivalents based on summer 2023 ru28 and ru40 deployments
+(concentration to saturation conversions are variable so these should be used with a grain of salt):
+
+hypoxia = < 3 mg/L, = 93.75 umol/L, approx = 40% saturation
+low DO = < 5 mg/L, = 156.25 umol/L, approx = 65% saturation
+supersaturation = > 100% saturation, approx = 7.5 mg/L, = 234.375
+
+    cm_oxy_mod: original cmocean colormap red to gray to yellow
+    cm_rygg: red to yellow to gray to green
+    cm_rogg: red to orange to gray to green (gray shading reversed from original oxy)
+    cm_partialturbo_r: red to orange/yellow to gray to blue
+
+Function to define categorical colormaps based on defined number of categories
+
+    cm_categorical: categorical colormaps
+
+Function to modify cyclical colormaps
+
+    cm_annualcycle: modified hsv with blues at ends
+
+"""
+
+def cm_oxy_mod(vmin=2, vmax=9, breaks=[3,7.5], red=True, gray=True, yellow=True):
+    """
+    Modify cmocean oxy colormap with defined break points
+    red (dark to less dark) to gray (dark to light) to yellow (light to dark-ish)
+
+    Args:
+        vmin (float): colormap minimum
+        vmax (float): colormap maximum
+        breaks (list of floats): break-points within colormap to transition to new color
+        red, gray, yellow (bool): whether to include each color in new colormap (default True for all)
+
+    Raises:
+        ValueError: nans or not-numbers provided for vmin, vmax, or breaks
+        ValueError: values do not increase from vmin -> breaks -> vmax
+        ValueError: number of breakpoints does not match with number of colors in new map
+
+    Returns:
+        object: matplotlib colormap
+    """
+    all_nums = np.append(np.append(vmin,breaks),vmax)
+    # check that all values are real numbers
+    if not np.all(np.isreal(all_nums)) or np.any(np.isnan(all_nums)):
+        raise ValueError('All values must be real numbers.')
+    # check that all values increase in the right order
+    if np.any(np.diff(all_nums)<0):
+        raise ValueError('All values must be increasing from vmin -> breaks -> vmax')
+    # check that number of break values works with number of colors
+    if red+gray+yellow != len(breaks)+1:
+        raise ValueError(f'Number of breakpoints ({len(breaks)}) must be one less than number of colors ({red+gray+yellow})')
+    # define interval
+    ni = (vmax-vmin)/100
+    
+    b=1
+    if red:
+        nints = int(np.floor((all_nums[b]-all_nums[b-1])/ni))
+        cm0 = cmocean.cm.oxy(np.linspace(0,.19,nints))
+        if b==1:
+            cmfull = cm0
+        else:
+            cmfull = np.vstack((cmfull,cm0))
+        b+=1
+    if gray:
+        nints = int(np.floor((all_nums[b]-all_nums[b-1])/ni))
+        cm0 = cmocean.cm.oxy(np.linspace(.2,.79,nints))
+        if b==1:
+            cmfull = cm0
+        else:
+            cmfull = np.vstack((cmfull,cm0))
+        b+=1
+    if yellow:
+        nints = int(np.floor((all_nums[b]-all_nums[b-1])/ni))
+        cm0 = cmocean.cm.oxy(np.linspace(.8,1,nints))
+        if b==1:
+            cmfull = cm0
+        else:
+            cmfull = np.vstack((cmfull,cm0))
+        b+=1
+    newmap = mcolors.LinearSegmentedColormap.from_list('my_colormap',cmfull)
+    return newmap
+
+
+def cm_rygg(vmin=2, vmax=9, breaks=[3,5,7.5], red=True, yellow=True, gray=True, green=True):
+    """
+    Modify cmocean oxy colormap with defined break points
+    red (dark to less dark) to yellow (light to dark-ish)to gray (dark to light) to green (light to mid)
+    
+    Args:
+        vmin (float): colormap minimum
+        vmax (float): colormap maximum
+        breaks (list of floats): break-points within colormap to transition to new color
+        red, yellow, gray, green (bool): whether to include each color in new colormap (default True for all)
+        
+    Raises:
+        ValueError: nans or not-numbers provided for vmin, vmax, or breaks
+        ValueError: values do not increase from vmin -> breaks -> vmax
+        ValueError: number of breakpoints does not match with number of colors in new map
+
+    Returns:
+        object: matplotlib colormap
+    """
+    all_nums = np.append(np.append(vmin,breaks),vmax)
+    # check that all values are real numbers
+    if not np.all(np.isreal(all_nums)) or np.any(np.isnan(all_nums)):
+        raise ValueError('All values must be real numbers.')
+    # check that all values increase in the right order
+    if np.any(np.diff(all_nums)<0):
+        raise ValueError('All values must be increasing from vmin -> breaks -> vmax')
+    # check that number of break values works with number of colors
+    if red+yellow+gray+green != len(breaks)+1:
+        raise ValueError(f'Number of breakpoints ({len(breaks)}) must be one less than number of colors ({red+yellow+gray+green})')
+    # define interval
+    ni = (vmax-vmin)/100
+    
+    b=1
+    if red:
+        nints = int(np.floor((all_nums[b]-all_nums[b-1])/ni))
+        cm0 = cmocean.cm.oxy(np.linspace(0,.19,nints))
+        if b==1:
+            cmfull = cm0
+        else:
+            cmfull = np.vstack((cmfull,cm0))
+        b+=1
+    if yellow:
+        nints = int(np.floor((all_nums[b]-all_nums[b-1])/ni))
+        cm0 = cmocean.cm.oxy(np.linspace(.8,1,nints))
+        if b==1:
+            cmfull = cm0
+        else:
+            cmfull = np.vstack((cmfull,cm0))
+        b+=1
+    if gray:
+        nints = int(np.floor((all_nums[b]-all_nums[b-1])/ni))
+        cm0 = cmocean.cm.oxy(np.linspace(.2,.79,nints))
+        if b==1:
+            cmfull = cm0
+        else:
+            cmfull = np.vstack((cmfull,cm0))
+        b+=1
+    if green:
+        nints = int(np.floor((all_nums[b]-all_nums[b-1])/ni))
+        cm0 = cmocean.cm.algae(np.linspace(0,.5,nints))
+        if b==1:
+            cmfull = cm0
+        else:
+            cmfull = np.vstack((cmfull,cm0))
+        b+=1
+    
+    newmap = mcolors.LinearSegmentedColormap.from_list('my_colormap',cmfull)
+    return newmap
+
+def cm_rogg(vmin=2, vmax=9, breaks=[3,5,7.5], red=True, orange=True, gray=True, green=True):
+    """
+    Modify cmocean oxy colormap with defined break points
+    red (dark to less dark) to orange (dark-ish to light)to gray (light to dark) to green (mid to light)
+
+    Args:
+        vmin (float): colormap minimum
+        vmax (float): colormap maximum
+        breaks (list of floats): break-points within colormap to transition to new color
+        red, orange, gray, green (bool): whether to include each color in new colormap (default True for all)
+        
+    Raises:
+        ValueError: nans or not-numbers provided for vmin, vmax, or breaks
+        ValueError: values do not increase from vmin -> breaks -> vmax
+        ValueError: number of breakpoints does not match with number of colors in new map
+
+    Returns:
+        object: matplotlib colormap
+    """
+    all_nums = np.append(np.append(vmin,breaks),vmax)
+    # check that all values are real numbers
+    if not np.all(np.isreal(all_nums)) or np.any(np.isnan(all_nums)):
+        raise ValueError('All values must be real numbers.')
+    # check that all values increase in the right order
+    if np.any(np.diff(all_nums)<0):
+        raise ValueError('All values must be increasing from vmin -> breaks -> vmax')
+    # check that number of break values works with number of colors
+    if red+orange+gray+green != len(breaks)+1:
+        raise ValueError(f'Number of breakpoints ({len(breaks)}) must be one less than number of colors ({red+orange+gray+green})')
+    # define interval
+    ni = (vmax-vmin)/100
+    
+    b=1
+    if red:
+        nints = int(np.floor((all_nums[b]-all_nums[b-1])/ni))
+        cm0 = cmocean.cm.oxy(np.linspace(0,.19,nints))
+        if b==1:
+            cmfull = cm0
+        else:
+            cmfull = np.vstack((cmfull,cm0))
+        b+=1
+    if orange:
+        nints = int(np.floor((all_nums[b]-all_nums[b-1])/ni))
+        cm0 = plt.cm.Oranges(np.linspace(.8,.4,nints))
+        if b==1:
+            cmfull = cm0
+        else:
+            cmfull = np.vstack((cmfull,cm0))
+        b+=1
+    if gray:
+        nints = int(np.floor((all_nums[b]-all_nums[b-1])/ni))
+        cm0 = cmocean.cm.oxy(np.linspace(.7,.2,nints))
+        if b==1:
+            cmfull = cm0
+        else:
+            cmfull = np.vstack((cmfull,cm0))
+        b+=1
+    if green:
+        nints = int(np.floor((all_nums[b]-all_nums[b-1])/ni))
+        cm0 = cmocean.cm.algae(np.linspace(.8,.3,nints))
+        if b==1:
+            cmfull = cm0
+        else:
+            cmfull = np.vstack((cmfull,cm0))
+        b+=1
+    
+    newmap = mcolors.LinearSegmentedColormap.from_list('my_colormap',cmfull)
+    return newmap
+
+def cm_partialturbo_r(vmin=2, vmax=9, breaks=[3,5,7.5], red=True, orange=True, gray=True, blue=True):
+    """
+    Modify cmocean oxy colormap with defined break points
+    red (dark to bright) to orange/yellow (orangey orange to orangey yellow)to gray (dark to light) to blue (cyan to bright)
+    
+    Args:
+        vmin (float): colormap minimum
+        vmax (float): colormap maximum
+        breaks (list of floats): break-points within colormap to transition to new color
+        red, orange, gray, blue (bool): whether to include each color in new colormap (default True for all)
+        
+    Raises:
+        ValueError: nans or not-numbers provided for vmin, vmax, or breaks
+        ValueError: values do not increase from vmin -> breaks -> vmax
+        ValueError: number of breakpoints does not match with number of colors in new map
+
+    Returns:
+        object: matplotlib colormap
+    """
+    all_nums = np.append(np.append(vmin,breaks),vmax)
+    # check that all values are real numbers
+    if not np.all(np.isreal(all_nums)) or np.any(np.isnan(all_nums)):
+        raise ValueError('All values must be real numbers.')
+    # check that all values increase in the right order
+    if np.any(np.diff(all_nums)<0):
+        raise ValueError('All values must be increasing from vmin -> breaks -> vmax')
+    # check that number of break values works with number of colors
+    if red+orange+gray+blue != len(breaks)+1:
+        raise ValueError(f'Number of breakpoints ({len(breaks)}) must be one less than number of colors ({red+orange+gray+blue})')
+    # define interval
+    ni = (vmax-vmin)/100
+    
+    b=1
+    if red:
+        nints = int(np.floor((all_nums[b]-all_nums[b-1])/ni))
+        cm0 = plt.cm.turbo(np.linspace(1,.85,nints))
+        if b==1:
+            cmfull = cm0
+        else:
+            cmfull = np.vstack((cmfull,cm0))
+        b+=1
+    if orange:
+        nints = int(np.floor((all_nums[b]-all_nums[b-1])/ni))
+        cm0 = plt.cm.turbo(np.linspace(.75,.6,nints))
+        if b==1:
+            cmfull = cm0
+        else:
+            cmfull = np.vstack((cmfull,cm0))
+        b+=1
+    if gray:
+        nints = int(np.floor((all_nums[b]-all_nums[b-1])/ni))
+        cm0 = cmocean.cm.oxy(np.linspace(.3,.79,nints))
+        if b==1:
+            cmfull = cm0
+        else:
+            cmfull = np.vstack((cmfull,cm0))
+        b+=1
+    if blue:
+        nints = int(np.floor((all_nums[b]-all_nums[b-1])/ni))
+        cm0 = plt.cm.turbo(np.linspace(.3,.1,nints))
+        if b==1:
+            cmfull = cm0
+        else:
+            cmfull = np.vstack((cmfull,cm0))
+        b+=1
+    
+    newmap = mcolors.LinearSegmentedColormap.from_list('my_colormap',cmfull)
+    return newmap
+
+
+def cm_categorical(N, listvals=True):
+    """
+    Define categorical colormap.
+    For 8 or fewer categories uses colorblind-friendly colormap defined by mpetroff (https://mpetroff.net/2018/03/color-cycle-picker/)
+    For 9-10 categories uses tab10
+    For 10-100 categories uses tab10 modified with multiple shades from ImportanceOfBeingErnest (https://stackoverflow.com/a/47232942/2643708)
+    For >100 categories no colormap defined
+    
+    Args:
+        N (int): number of categories to include in colorbar
+        listvals (bool): return array of rgb or hex values instead of colormap object
+        
+    Raises:
+        ValueError: Too many categories for colormap
+
+    Returns:
+        object: matplotlib colormap or list of colors
+    """
+    if N > 100:
+        raise ValueError("Too many categories for colormap.")
+    
+    if N <= 8:
+        # unicorn colormap from mpetroff
+        newmap = np.array(['#7b85d4','#f37738','#83c995','#d7369e','#c4c9d8','#859795','#e9d043','#ad5b50'])
+    else:
+        # modified tab10 from ImportanceOfBeingErnest
+        # if N<=10, equivalent to original tab10
+        nc = 10
+        nsc = int(np.ceil(N/nc))
+        ccolors = plt.cm.tab10(np.arange(nc, dtype=int))
+        for i, c in enumerate(ccolors):
+            chsv = mcolors.rgb_to_hsv(c[:3])
+            arhsv = np.tile(chsv,nsc).reshape(nsc,3)
+            arhsv[:,1] = np.linspace(chsv[1],0.25,nsc)
+            arhsv[:,2] = np.linspace(chsv[2],1,nsc)
+            rgb = mcolors.hsv_to_rgb(arhsv)
+            if i==0:
+                newmap = rgb
+            else:
+                newmap = np.vstack((newmap,rgb))
+            if i==np.mod(N,nc) and i>nc*(nsc-1) and np.mod(N,nc)!=0:
+                nsc-=1
+
+    newmap = newmap[:N]
+
+    if not listvals:
+        newmap = mcolors.ListedColormap(newmap)
+    
+    return newmap
+
+
+def cm_annualcycle(step='monthly', listvals=False):
+    """
+    Define annual cycle colormap.
+    HSV modified with blues at beginning and end of colormap (winter), reds in middle (summer)
+    
+    Args:
+        step (str): 'monthly' (default, 12 steps in colormap) or 'daily' (365 steps in colormap)
+        listvals (bool): return array of rgb or hex values instead of colormap object
+        
+    Raises:
+        ValueError: value for 'step' not a valid option
+
+    Returns:
+        object: matplotlib colormap or list of colors
+    """
+    if step not in ['monthly', 'daily']:
+        raise ValueError("Too many categories for colormap.")
+    
+    if step == 'monthly':
+        nints=12
+        cm_jan_apr = plt.cm.hsv(np.linspace(.7,.2,4))
+        cm_may_aug = plt.cm.hsv(np.linspace(.19,0,4))
+        cm_sep_dec = plt.cm.hsv(np.linspace(.9,.71,4))
+    elif step=='daily':
+        nints=365
+        cm_jan_apr = plt.cm.hsv(np.linspace(.7,.2,120))
+        cm_may_aug = plt.cm.hsv(np.linspace(.19,0,123))
+        cm_sep_dec = plt.cm.hsv(np.linspace(.9,.71,122))
+    newmap = np.vstack((cm_jan_apr, cm_may_aug, cm_sep_dec))
+
+    if not listvals:
+        newmap = mcolors.ListedColormap(newmap)
+    
+    return newmap
diff --git a/cool_maps/download.py b/cool_maps/download.py
index 6d9658e..41f1d03 100644
--- a/cool_maps/download.py
+++ b/cool_maps/download.py
@@ -1,15 +1,23 @@
 from erddapy import ERDDAP
+import os
+import xarray as xr
+import pandas as pd
+import numpy as np
 
 def get_bathymetry(extent=(-100, -45, 5, 46),
                    server="https://hfr.marine.rutgers.edu/erddap/",
-                   dataset_id="bathymetry_srtm15_v24"
+                   dataset_id="bathymetry_srtm15_v24",
+                   file=None
                    ):
     """
     Function to select bathymetry within a bounding box.
     This function pulls GEBCO 2014 bathy data from hfr.marine.rutgers.edu 
+    OR from multi-source elevation/bathy data from CF compliant NetCDF file downloaded from https://www.gmrt.org/GMRTMapTool/
 
     Args:
         extent (tuple, optional): Cartopy bounding box. Defaults to (-100, -45, 5, 46).
+        file (str filename, optional): CF Compliant NetCDF file containing GMRT bathymetry
+                                       if None (default) or if file is not found, will default to ERDDAP
 
     Returns:
         xarray.Dataset: xarray Dataset containing bathymetry data
@@ -17,6 +25,13 @@ def get_bathymetry(extent=(-100, -45, 5, 46),
     lons = extent[:2]
     lats = extent[2:]
 
+    if file and os.path.isfile(file):
+        bathy = xr.open_dataset(file)
+        bathy = bathy.sel(lon=slice(min(lons), max(lons)))
+        bathy = bathy.sel(lat=slice(min(lats), max(lats)))
+        bathy = bathy.rename({'lon':'longitude','lat':'latitude','altitude':'z'})
+        return bathy
+
     e = ERDDAP(
         server=server,
         protocol="griddap"
@@ -41,14 +56,16 @@ def get_totals_from_erddap(time_start,
                            server="https://hfr.marine.rutgers.edu/erddap/",
                            dataset_id='5MHz_6km_realtime-agg_archive_v00'):
     """
-    Function to select bathymetry within a bounding box.
+    Function to select HFR data within a bounding box.
     This function pulls GEBCO 2014 bathy data from hfr.marine.rutgers.edu 
 
     Args:
-        bbox (list, optional): Cartopy bounding box. Defaults to None.
+        time_start (str): Start time to read.
+        time_end (str, optional): End time to read. Defaults to None/latest.
+        extent (tuple, optional): Cartopy bounding box. Defaults to None.
 
     Returns:
-        xarray.Dataset: xarray Dataset containing bathymetry data
+        xarray.Dataset: xarray Dataset containing HFR data
     """
     
     if not time_end:
@@ -77,4 +94,48 @@ def get_totals_from_erddap(time_start,
     e.constraints["time>="] = time_end
         
     # return xarray dataset
-    return e.to_xarray()
\ No newline at end of file
+    return e.to_xarray()
+
+def get_glider_bathymetry(deployment,
+                          time_start=None,
+                          time_end=None):
+    """
+    Function to select bathymetry associated with a glider deployment.
+    This function pulls glider-measured bathymetry data from slocum-data.marine.rutgers.edu and removes spikes
+
+    Args:
+        deployment (str): Name of deployment to get bathymetry for (glider-yyyymmddTHHMM)
+        time_start (str, optional): Start time. Defaults to None/beginning of deployment
+        time_end (str, optional): End time. Defaults to None/end of deployment
+
+    Returns:
+        pandas.DataFrame: pandas DataFrame containing bathymetry data
+    """
+    
+    ru_erddap = ERDDAP(server='http://slocum-data.marine.rutgers.edu/erddap', protocol='tabledap')
+    glider_url = ru_erddap.get_search_url(search_for=f'{deployment}-trajectory', response='csv')
+    try:
+        glider_datasets = pd.read_csv(glider_url)['Dataset ID']
+    except:
+        print(f'Unable to access deployment {deployment}. Check for typos.')
+        return None
+
+    ru_erddap.dataset_id = glider_datasets[0]
+    ru_erddap.constraints = {'m_water_depth>': 2}
+    if time_start:
+        ru_erddap.constraints["time<="] = time_start
+    if time_end:
+        ru_erddap.constraints["time>="] = time_end
+    ru_erddap.variables = ['time', 'm_water_depth']
+    glider_bathy = ru_erddap.to_pandas()
+    glider_bathy.columns=['time', 'water_depth']
+    glider_bathy['time']=pd.to_datetime(glider_bathy['time'])
+    glider_bathy = glider_bathy.sort_values(by='time')
+    glider_bathy['d1']=np.nan
+    glider_bathy['d2']=np.nan
+    glider_bathy.loc[:,'d1']=np.append(np.divide((np.array(glider_bathy['water_depth'][1:]) - np.array(glider_bathy['water_depth'][:len(glider_bathy)-1])), (np.array(glider_bathy['time'][1:].astype('int64')//1e9) - np.array(glider_bathy['time'][:len(glider_bathy)-1].astype('int64')//1e9))), np.nan)
+    glider_bathy.loc[1:,'d2']=np.divide((np.array(glider_bathy['water_depth'][1:]) - np.array(glider_bathy['water_depth'][:len(glider_bathy)-1])), np.array(glider_bathy['time'][1:].astype('int64')//1e9) - np.array(glider_bathy['time'][:len(glider_bathy)-1].astype('int64')//1e9))
+    glider_bathy['d']=np.max(np.abs(glider_bathy[['d1','d2']]), axis=1)
+    glider_bathy=glider_bathy[glider_bathy['d']<.2][['time','water_depth']]
+    
+    return glider_bathy
diff --git a/cool_maps/plot.py b/cool_maps/plot.py
index a5cc891..95f0773 100644
--- a/cool_maps/plot.py
+++ b/cool_maps/plot.py
@@ -11,10 +11,11 @@
 import matplotlib.ticker as mticker
 import numpy as np
 from oceans.ocfis import spdir2uv, uv2spdir
+import cmocean as cmo
 
 # Imports from cool_maps
-from cool_maps.calc import calculate_ticks, dd2dms, fmt
-from cool_maps.download import get_bathymetry
+from cool_maps.calc import calculate_ticks, dd2dms, fmt, calculate_colorbar_ticks
+from cool_maps.download import get_bathymetry, get_glider_bathymetry
 
 
 # Suppresing warnings for a "pretty output."
@@ -29,6 +30,9 @@
 
 def add_bathymetry(ax, lon, lat, elevation,
                    levels=(-1000),
+                   method='contour',
+                   legend_scale='both',
+                   fontsize=13,
                    zorder=5,
                    transform=proj['data'], 
                    transform_first=False):
@@ -41,17 +45,116 @@ def add_bathymetry(ax, lon, lat, elevation,
         lat (array-like): Latitudes of bathymetry
         elevation (array-like): Elevation of bathymetry
         levels (tuple, optional): Determines the number and positions of the contour lines. Defaults to (-1000).
+        method (string): Method for plotting bathymetry. Defaults to contour. Options:
+            - contour: standard black contour at :levels:
+            - shadedcontour: contours in shades of gray varying with depth
+            - blues: pcolormesh using Blues colormap; excludes land and ignores :levels:
+            - blues_log: pcolormesh of log-transformed bathymetry using Blues colormap; excludes land and ignores :levels:
+            - topo: pcolormesh using cmocean topo colormap; excludes land and ignores :levels:
+            - topo_log: pcolormesh of log-transformed bathymetry using cmocean topo colormap; excludes land and ignores :levels:
+            - topofull: pcolormesh using cmocean topo colormap; includes land and ignores :levels:
+            - topofull_log: pcolormesh of log-transformed altitude/bathymetry using cmocean topo colormap; includes land and ignores :levels:
+        legend_scale (string, optional): Measurement system to use for legend. Currently only supported for shadedcontour; no legend for other options. ASSUMES LEVELS PROVIDED ARE IN METERS.
+            - metric: meters
+            - imperial: fathoms
+            - both: meters and fathoms
+            - off: no legend
+        fontsize (int, optional): Font size for legend
         zorder (int, optional): Drawing order for this function on the axes. Defaults to 5.
         transform (_type_, optional): Coordinate system data is defined in. Defaults to crs.PlateCarree.
         transform_first (bool, optional): Indicate that Cartopy points should be transformed before calling the contouring algorithm, which can have a significant impact on speed (it is much faster to transform points than it is to transform patches). Defaults to False.
+
+    Raises:
+        NameError: provided method is not recognized
+
+    Returns:
+        object: plotted bathymetry handle
     """
+    recognized_methods=['contour', 'shadedcontour', 'blues', 'blues_log', 'topo', 'topo_log', 'topofull', 'topofull_log']
+    recognized_legends=['metric', 'imperial', 'both', 'off']
+    if method not in recognized_methods:
+        raise NameError(f'{method} is not a currently supported option for bathymetry plotting. Please choose from {", ".join(recognized_methods)}')
+    if legend_scale not in recognized_legends:
+        raise NameError(f'{legend_scale} is not a currently supported option for the legend. Please choose from {", ".join(recognized_legends)}')
+
     lons, lats = np.meshgrid(lon, lat)
-    h = ax.contour(lons, lats, elevation, levels, 
-                    linewidths=.75, alpha=.5, colors='k', 
-                    transform=transform, 
-                    transform_first=transform_first, # May speed plotting up
-                    zorder=zorder)
-    ax.clabel(h, levels, inline=True, fontsize=6, fmt=fmt)
+    elevation=elevation.copy()
+
+    if method=='contour':
+        h = ax.contour(lons, lats, elevation, levels, 
+                        linewidths=.75, alpha=.5, colors='k', 
+                        transform=transform, 
+                        transform_first=transform_first, # May speed plotting up
+                        zorder=zorder)
+        ax.clabel(h, levels, inline=True, fontsize=6, fmt=fmt)
+    
+    if method=='shadedcontour':
+        ci=1.0/float(len(levels))
+        bathyLabels = {-x: '{}m'.format(int(x)) for x in np.sort(np.abs(levels))[::-1]}
+        if legend_scale == 'imperial':
+            bathyLabels = {-x: '{}fth'.format(int(x * 0.546807)) for x in np.sort(np.abs(levels))[::-1]}
+        if legend_scale == 'both':
+            bathyLabels = {-x: '{}fth, {}m'.format(int(x * 0.546807), int(x)) for x in np.sort(np.abs(levels))[::-1]}
+        rgb_col=[(ci*cs,ci*cs,ci*cs) for cs in range(len(bathyLabels))]
+        h=[]
+        for cs in range(len(bathyLabels)):
+            h = np.append(h, ax.contour(lons,lats,elevation,
+                       levels=[list(bathyLabels.keys())[cs]],
+                       linestyles='solid',linewidths=.75,
+                       colors=[(ci*cs,ci*cs,ci*cs)],
+                       zorder=zorder,transform=transform))
+        if legend_scale != 'off':
+            cs_cols=[plt.Line2D((0,1),(0,1),color=pc) for pc in rgb_col]
+            ax.legend(cs_cols[::-1],list(bathyLabels.values())[::-1],
+                      loc='upper center',
+                      bbox_to_anchor=(0.5,-0.05),
+                      fancybox=True,ncol=len(levels),
+                      fontsize=fontsize)
+    
+    if method in ['blues_log', 'topo_log', 'topofull_log']:
+        elevation[np.abs(elevation)<1] = 0
+        elevation[elevation>0] = np.log10(elevation[elevation>0])
+        elevation[elevation<0] = -np.log10(np.abs(elevation[elevation<0]))
+    if method in ['blues', 'topo', 'blues_log', 'topo_log']:
+        elevation[elevation>0] = np.nan
+    if method in ['blues', 'blues_log']:
+        cmap = plt.cm.Blues_r
+        vmin = np.nanquantile(elevation, 0.05)
+        vmax = 0
+    if method in ['topo', 'topo_log', 'topofull', 'topofull_log']:
+        cmap = cmo.cm.topo
+        vmin = -np.nanquantile(np.abs(elevation), 0.95)
+        vmax = np.nanquantile(np.abs(elevation), 0.95)
+    if method in ['blues', 'blues_log', 'topo', 'topo_log', 'topofull', 'topofull_log']:
+        h = plt.pcolormesh(lons, lats, elevation, cmap=cmap, vmin=vmin, vmax=vmax, 
+                           transform=transform, zorder=zorder)
+
+    return h
+
+
+def add_glider_bathymetry(ax, deployment,
+                          time_start=None,
+                          time_end=None,
+                          color='black'):
+    """
+    Download and plot bathymetry measured by glider during a given deployment
+
+    Args:
+        ax (matplotlib.axes): matplotlib axes
+        deployment (str): name of deployment to grab and plot bathymetry for
+        time_start (str, optional): Start time. Defaults to None/beginning of deployment
+        time_end (str, optional): End time. Defaults to None/end of deployment
+        color (str, optional): name of color to plot bathymetry
+    Returns:
+        object: Patch
+    """
+
+    glider_bathy = get_glider_bathymetry(deployment, time_start=time_start, time_end=time_end)
+    floor_depth = np.max(glider_bathy['water_depth'])*1.05
+    h = plt.fill(np.append(glider_bathy['time'], [max(glider_bathy['time']), min(glider_bathy['time'])]), 
+               np.append(glider_bathy['water_depth'], [floor_depth, floor_depth]), 
+               color=color)
+
     return h
 
 
@@ -157,6 +260,61 @@ def add_features(ax,
     ax.add_feature(cfeature.BORDERS, linestyle='--', zorder=zorder+10.3)
 
 
+def add_double_temp_colorbar(ax, h, vmin, vmax,
+                             anomaly=False, fontsize=13):
+    """
+    Add colorbar with Celsius and Fahrenheit units
+    https://pythonmatplotlibtips.blogspot.com/2019/07/draw-two-axis-to-one-colorbar.html
+
+    Args:
+        ax (matplotlib.axes): matplotlib axes
+        h (matplotlib object handle to match colorbar to): 
+        vmin (float): minimum value of colorbar (should match vmin of object colorbar is mapped to)
+        vmax (float): maximum value of colorbar (should match vmax of object colorbar is mapped to)
+        anomaly (bool): whether product is an anomaly
+        fontsize (int, optional): font size for tick labels
+
+    Returns:
+        object: Colorbar
+        axes: Twin axes for colorbar
+    """
+        
+    cbticks = calculate_colorbar_ticks(vmin, vmax, c0=anomaly)
+    if anomaly:
+        cbticksF = calculate_colorbar_ticks(vmin*1.8, vmax*1.8, c0=anomaly)
+    else:
+        cbticksF = calculate_colorbar_ticks(vmin*1.8+32, vmax*1.8+32, c0=anomaly)
+    cbCLabels=[str(int(cbticks[i]))+u"\N{DEGREE SIGN}"+"C" for i in range(len(cbticks))]
+    cbFLabels = [str(int(cbticksF[i])) + u"\N{DEGREE SIGN}" + "F" for i in range(len(cbticksF))]
+
+    cb = plt.colorbar(h) #, ticks=cbticks)
+    cb.ax.set_yticks(cbticks, labels=cbCLabels, fontsize=fontsize)
+    pcb=cb.ax.get_position()
+    pax=ax.get_position()
+
+    # set up axis overlapping colorbar to have degree F and degree C labels
+    cb.ax.set_aspect('auto')
+    pcb.x0 = pax.x1+.055
+    pcb.x1 = pax.x1+.085
+    pcb.y0 = pax.y0
+    pcb.y1 = pax.y1
+    cb2 = cb.ax.twinx()
+    ax.set_position(pax)
+    cb2.set_ylim(np.array([vmin, vmax])*1.8+(32*(not anomaly)))
+    cb2.yaxis.set_label_position('left')
+    cb2.yaxis.set_ticks_position('left')
+    cb.ax.yaxis.set_label_position('right')
+    cb.ax.yaxis.set_ticks_position('right')
+    cb2.set_yticks(cbticksF, labels=cbFLabels, fontsize=fontsize)
+    cb.ax.set_position(pcb)
+    cb2.set_position(pcb)
+    cb2.spines['right'].set_visible(False)
+    cb2.spines['top'].set_visible(False)
+    cb2.spines['bottom'].set_visible(False)
+
+    return cb, cb2
+
+
 def add_ticks(ax, extent,
               proj=proj['data'], 
               fontsize=13, 
@@ -164,7 +322,8 @@ def add_ticks(ax, extent,
               label_right=False,
               label_bottom=True, 
               label_top=False, 
-              gridlines=False):
+              gridlines=False,
+              decimal_degrees=False):
     """
     Calculate and add nicely formatted ticks to your map
 
@@ -174,15 +333,16 @@ def add_ticks(ax, extent,
         proj (cartopy.crs class, optional): Define a projected coordinate system for ticks. Defaults to ccrs.PlateCarree().
         fontsize (int, optional): Font size of tick labels. Defaults to 13.
         gridlines (bool, optional): Add gridlines to map. Defaults to False.
+        decimal_degrees (bool, optional): Label axes with decimal degrees instead of degree-minute-second. Defaults to False.
     """
     # Calculate Longitude ticks
-    tick0x, tick1, ticklab = calculate_ticks(extent, 'longitude')
+    tick0x, tick1, ticklab = calculate_ticks(extent, 'longitude', decimal_degrees=decimal_degrees)
     ax.set_xticks(tick0x, minor=True, crs=proj)
     ax.set_xticks(tick1, crs=proj)
     ax.set_xticklabels(ticklab, fontsize=fontsize)
 
     # Calculate Latitude Ticks
-    tick0y, tick1, ticklab = calculate_ticks(extent, 'latitude')
+    tick0y, tick1, ticklab = calculate_ticks(extent, 'latitude', decimal_degrees=decimal_degrees)
     ax.set_yticks(tick0y, minor=True, crs=proj)
     ax.set_yticks(tick1, crs=proj)
     ax.set_yticklabels(ticklab, fontsize=fontsize)
@@ -224,12 +384,15 @@ def create(extent,
            gridlines=False,
            bathymetry=False,
            isobaths=(-1000, -100),
+           bathymetry_method='contour',
+           bathymetry_file=None,
            xlabel=None,
            ylabel=None,
            tick_label_left=True,
            tick_label_right=False,
            tick_label_bottom=True,
            tick_label_top=False,
+           decimal_degrees=False,
            labelsize=14,
            ax=None,
            figsize=(11,8),
@@ -249,12 +412,15 @@ def create(extent,
         gridlines (bool, optional): Add gridlines. Defaults to False
         bathymetry (bool or tuple, optional): Download and plot bathymetry on map. Defaults to False.
         isobaths (tuple or list, optional): Elevation at which to create bathymetric contour lines. Defaults to (-1000, -100)
+        bathymetry_method (str, optional): Method for plotting bathymetry (see cool_maps.plot.plot_bathymetry for options). Defaults to contour.
+        bathymetry_file (str filename or None): Name of CF compliant GMRT nc file, or None (default) to use ERDDAP.
         xlabel (str, optional): X Axis Label. Defaults to None
         ylabel (str, optional): Y Axis Label. Defaults to None
         tick_label_left (bool, optional): Add tick labels to left side of plot. Defaults to True.
         tick_label_right (bool, optional): Add tick labels to right side of plot. Defaults to False.
         tick_label_bottom (bool, optional): Add tick labels to bottom side of plot. Defaults to True.
         tick_label_top (bool, optional): Add tick labels to top side of plot. Defaults to False.
+        decimal_degrees (bool, optional): Label axes with decimal degrees instead of degree-minute-second. Defaults to False.
         labelsize (int, optional): Font size for axis labels. Defaults to 14.
         ax (matplotlib.Axis, optional): Pass matplotlib axis to function. Not necessary if plotting to subplot. Defaults to None.
         figsize (tuple, optional): (width, height) of the figure. Defaults to (11,8).
@@ -295,12 +461,19 @@ def create(extent,
 
     # Add bathymetry
     if bathymetry:
+        if 'contour' in bathymetry_method:
+            bathy_zorder_add=99
+        elif 'topofull' in bathymetry_method:
+            bathy_zorder_add=20
+        else:
+            bathy_zorder_add=5
         bargs = {
             "levels": isobaths,
-            "zorder": zorder+99,
+            "zorder": zorder+bathy_zorder_add,
+            "method": bathymetry_method,
         }
-        bathy = get_bathymetry(extent)
-        add_bathymetry(ax, bathy['longitude'], bathy['latitude'], bathy['z'], **bargs)
+        bathy = get_bathymetry(extent, file=bathymetry_file)
+        add_bathymetry(ax, bathy['longitude'].data, bathy['latitude'].data, bathy['z'].data, **bargs)
 
     # Add ticks using custom functions
     if ticks:
@@ -310,6 +483,7 @@ def create(extent,
         tick_dict['label_bottom'] = tick_label_bottom
         tick_dict['label_top'] = tick_label_top
         tick_dict['gridlines'] = gridlines
+        tick_dict['decimal_degrees'] = decimal_degrees
         add_ticks(ax, extent, **tick_dict)
     else:
         # Add gridlines using built-in cartopy gridliner, provided we didn't