test plotting diurnal with adcp data.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import xarray as xr\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib\n",
    "import numpy.ma as ma\n",
    "import numpy\n",
    "import datetime as dt\n",
    "from datetime import datetime, timedelta\n",
    "import matplotlib.cm as cm\n",
    "from mpl_toolkits.basemap import Basemap\n",
    "import matplotlib.pyplot as plt\n",
    "from palettable import colorbrewer\n",
    "import numpy as np\n",
    "\n",
    "\n",
    "def smooth(x,window_len=11,window='hanning'):\n",
    "    \"\"\"smooth the data using a window with requested size.\n",
    "    \n",
    "    This method is based on the convolution of a scaled window with the signal.\n",
    "    The signal is prepared by introducing reflected copies of the signal \n",
    "    (with the window size) in both ends so that transient parts are minimized\n",
    "    in the begining and end part of the output signal.\n",
    "    \n",
    "    input:\n",
    "        x: the input signal \n",
    "        window_len: the dimension of the smoothing window; should be an odd integer\n",
    "        window: the type of window from 'flat', 'hanning', 'hamming', 'bartlett', 'blackman'\n",
    "            flat window will produce a moving average smoothing.\n",
    "\n",
    "    output:\n",
    "        the smoothed signal\n",
    "        \n",
    "    example:\n",
    "\n",
    "    t=linspace(-2,2,0.1)\n",
    "    x=sin(t)+randn(len(t))*0.1\n",
    "    y=smooth(x)\n",
    "    \n",
    "    see also: \n",
    "    \n",
    "    numpy.hanning, numpy.hamming, numpy.bartlett, numpy.blackman, numpy.convolve\n",
    "    scipy.signal.lfilter\n",
    " \n",
    "    TODO: the window parameter could be the window itself if an array instead of a string\n",
    "    NOTE: length(output) != length(input), to correct this: return y[(window_len/2-1):-(window_len/2)] instead of just y.\n",
    "    \"\"\"\n",
    "\n",
    "    if x.ndim != 1:\n",
    "        raise ValueError( \"smooth only accepts 1 dimension arrays.\")\n",
    "\n",
    "    if x.size < window_len:\n",
    "        raise ValueError(\"Input vector needs to be bigger than window size.\")\n",
    "\n",
    "\n",
    "    if window_len<3:\n",
    "        return x\n",
    "\n",
    "\n",
    "    if not window in ['flat', 'hanning', 'hamming', 'bartlett', 'blackman']:\n",
    "        raise ValueError(\"Window is on of 'flat', 'hanning', 'hamming', 'bartlett', 'blackman'\")\n",
    "\n",
    "\n",
    "    s=numpy.r_[x[window_len-1:0:-1],x,x[-2:-window_len-1:-1]]\n",
    "    #print(len(s))\n",
    "    if window == 'flat': #moving average\n",
    "        w=numpy.ones(window_len,'d')\n",
    "    else:\n",
    "        w=eval('numpy.'+window+'(window_len)')\n",
    "\n",
    "    y=numpy.convolve(w/w.sum(),s,mode='valid')\n",
    "    return y\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "adir_adcp = 'F:/data/cruise_data/saildrone/baja-2018/adcp_data/'\n",
    "filename_adcp = adir_adcp + 'saildrone-gen_4-baja_2018-sd1002-adcp-20180411T175730-20180611T235730_300s-v1.nc'\n",
    "ds_adcp=xr.open_dataset(filename_adcp)\n",
    "ds_adcp.close()\n",
    "ds_adcp = ds_adcp.isel(trajectory=0)\n",
    "ds_adcp=ds_adcp.swap_dims({'obs':'time'})\n",
    "print(ds_adcp)\n",
    "time_offset_to_lmt=(ds_adcp.longitude/360.)*24.*60\n",
    "tlmt= ds_adcp.time.data+time_offset_to_lmt*np.timedelta64(1,'m')# dt.timedelta(seconds=1)\n",
    "ds_adcp['tlmt']=tlmt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "adir='F:/data/cruise_data/saildrone/baja-2018/'\n",
    "#adir='C:/Users/gentemann/Desktop/data/'\n",
    "fname=adir+'saildrone-gen_4-baja_2018-sd1002-20180411T180000-20180611T055959-10_minutes-v1_withloggers_matchedadcp_v3.nc'\n",
    "#fname=adir+'saildrone-gen_4-baja_2018-EP-sd1002-ALL-1_min-v1_with_loggers.nc'\n",
    "ds=xr.open_dataset(fname)\n",
    "ds.close()\n",
    "ds = ds.isel(obs=slice(0,8674))\n",
    "ds=ds.swap_dims({'obs':'time'})\n",
    "lons_usv=ds.longitude\n",
    "lats_usv=ds.latitude\n",
    "time_usv=ds.time\n",
    "xlon = ds.longitude\n",
    "ds.sea_water_temperature_logger1_mean\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#calculate local time\n",
    "time_offset_to_lmt=(xlon/360.)*24.*60\n",
    "tlmt= ds.time.data+time_offset_to_lmt*np.timedelta64(1,'m')# dt.timedelta(seconds=1)\n",
    "ds['tlmt']=tlmt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "z=[-0.295,-0.6,-0.985,-1.42,-1.785]\n",
    "i1save=[12800,35800,46000,61400,62080]\n",
    "i2save=[13300,36300,46700,62080,63400]\n",
    "file_str=['p1','p2','p3','p4','p5','p6','p7']\n",
    "ilen=len(i1save)\n",
    "ilen_ave=ds.sea_water_temperature_logger1_mean.size\n",
    "tdata=np.ones([ilen_ave,5])\n",
    "for i in range(0,ilen_ave):\n",
    "    tdata[i,0]=ds.sea_water_temperature_logger1_mean[i]-ds.sea_water_temperature_logger4_mean[i]\n",
    "    tdata[i,1]=ds.TEMP_CTD_MEAN[i]-(ds.sea_water_temperature_logger4_mean[i]) #-273.15)\n",
    "    tdata[i,2]=ds.sea_water_temperature_logger2_mean[i]-ds.sea_water_temperature_logger4_mean[i]\n",
    "    tdata[i,3]=ds.sea_water_temperature_logger3_mean[i]-ds.sea_water_temperature_logger4_mean[i]\n",
    "    tdata[i,4]=ds.sea_water_temperature_logger4_mean[i]-ds.sea_water_temperature_logger4_mean[i]\n",
    "#create xarray structure\n",
    "profiles = xr.DataArray(tdata, coords={'time': ds.time,'depth':z}, dims=('time','depth'))\n",
    "profiles_lmt = xr.DataArray(ds.tlmt, name='tlmt',coords={'time': ds.time}, dims=('time'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#print(profiles_lmt)\n",
    "print(profiles_lmt[0].data,profiles.time[0].data)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib as mpl\n",
    "from palettable import colorbrewer\n",
    "import matplotlib.cm as cm\n",
    "#make plots of times with diurnal warming\n",
    "#add wind speed to figure\n",
    "#i1save=[12800,35800,46000,61400,62000]\n",
    "#i2save=[13300,36300,46700,62080,63400]\n",
    "#i1save=[12800,35800,46000,61600,63080]\n",
    "#i2save=[13300,36300,46700,63080,63800]\n",
    "i1save=[np.datetime64('2018-04-20T06:00'),np.datetime64('2018-05-06T06:00'),np.datetime64('2018-05-13T06:00'),np.datetime64('2018-05-24T06:00'),np.datetime64('2018-05-25T06:00')]\n",
    "i2save=[np.datetime64('2018-04-21T06:00'),np.datetime64('2018-05-07T06:00'),np.datetime64('2018-05-14T06:00'),np.datetime64('2018-05-25T06:00'),np.datetime64('2018-05-26T06:00')]\n",
    "#i2save=[np.datetime64('2018-04-21'),np.datetime64('2018-05-07')]\n",
    "                                                 \n",
    "\n",
    "file_str=['p1','p2','p3','p4','p5','p6','p7']\n",
    "ilen=len(i1save)\n",
    "print(ilen)\n",
    "for ii in range(0,ilen):\n",
    "    i1=i1save[ii]\n",
    "    i2=i2save[ii]\n",
    "\n",
    "  #  t=tlmt.sel(time=slice(i1,i2)) #local time\n",
    "  #  tem2=(t[0].values-np.datetime64('1970-01-01'))/np.timedelta64(1,'s')\n",
    "  #  tem3=datetime.utcfromtimestamp(tem2)\n",
    "  #  day_of_year = tem3.timetuple().tm_yday\n",
    "\n",
    "    subset_adcp = ds_adcp.sel(time=slice(i1,i2))\n",
    "    subset_adcp['vel_horiz']=np.sqrt(subset_adcp.vel_east**2+subset_adcp.vel_north**2)\n",
    "    subset_adcp_masked = subset_adcp.where(subset_adcp.percent_good>10)      \n",
    "    subset_usv = ds.sel(time=slice(i1,i2))\n",
    "    \n",
    "#    t=subset_usv.time\n",
    "    t=subset_usv.tlmt\n",
    "    s0=subset_usv.TEMP_CTD_MEAN.sel(time=slice(i1,i2))\n",
    "    s1=subset_usv.sea_water_temperature_logger1_mean\n",
    "    s2=subset_usv.sea_water_temperature_logger2_mean\n",
    "    s3=subset_usv.sea_water_temperature_logger3_mean\n",
    "    s4=subset_usv.sea_water_temperature_logger4_mean\n",
    "    wnd=(subset_usv.UWND_MEAN**2+subset_usv.VWND_MEAN**2)**.5\n",
    "    s0mask = np.isfinite(s0)\n",
    "    s1mask = np.isfinite(s1)\n",
    "    s2mask = np.isfinite(s2)\n",
    "    s3mask = np.isfinite(s3)\n",
    "    s4mask = np.isfinite(s4)\n",
    "    wmask = np.isfinite(wnd)\n",
    "    \n",
    "    fig = plt.figure(figsize=(8,10))\n",
    "    ax = fig.add_axes([0.1, 0.8, 0.8, 0.2])\n",
    "    ax2 = ax.twinx()\n",
    "    ax.plot_date(t[s1mask], s1[s1mask],'-',markersize=.5,lw=1)\n",
    "    ax.plot_date(t[s0mask], s0[s0mask],'-',markersize=.5,lw=1)\n",
    "    ax.plot_date(t[s2mask], s2[s2mask],'-',markersize=.5,lw=1)\n",
    "    ax.plot_date(t[s3mask], s3[s3mask],'-',markersize=.5,lw=1)\n",
    "    ax.plot_date(t[s4mask], s4[s4mask],'-',markersize=.5,lw=1)       \n",
    "    ax.set_prop_cycle(None)\n",
    "#    fig.gca().set_prop_cycle(None)   \n",
    "    [y1,y2]=ax.get_ylim()\n",
    "    [x1,x2]=ax.get_xlim()\n",
    "    for ihr in range(3,24,3):\n",
    "        tt=t[0]+np.timedelta64(ihr,'h')\n",
    "        ax.plot([tt.data,tt.data],[y1,y2])\n",
    "    ax2.plot_date(t[wmask], wnd[wmask],'k-',markersize=.5,lw=2)\n",
    "    ax.set(xlabel='Date (LMT)', ylabel='SST (C)')\n",
    "    ax2.set(ylabel='Wind Speed (ms$^{-1}$)')\n",
    "    # rotate and align the tick labels so they look better\n",
    "    fig.autofmt_xdate()\n",
    "    ax.set_xlim(t[0].data,t[-1].data)\n",
    "    #ax.set_xlim(subset)\n",
    "    #ax.fmt_xdata = t.DateFormatter('%Y-%m-%d')\n",
    "    tem=str(ds.TEMP_CTD_MEAN.installed_height)+' m'\n",
    "    ax.legend([ds.sea_water_temperature_logger1_mean.installed_height,tem,ds.sea_water_temperature_logger2_mean.installed_height,ds.sea_water_temperature_logger3_mean.installed_height,ds.sea_water_temperature_logger4_mean.installed_height])\n",
    "    ax.text(t[0].data+np.timedelta64(10,'m'),s1[s1mask].max().data,'A',size=12,weight='bold')\n",
    "\n",
    "    #fig, ax = plt.subplots()\n",
    "    \n",
    "    profile_subset = profiles.sel(time=slice(i1,i2))\n",
    "    profile_subset_lmt = profiles_lmt.sel(time=slice(i1,i2))\n",
    "    fig.gca().set_prop_cycle(None)\n",
    "    ax = fig.add_axes([0.1, 0.55, 0.8, 0.2])   \n",
    "    x = profile_subset_lmt.data\n",
    "    y = profile_subset.depth\n",
    "    X, Y = np.meshgrid(x, y)\n",
    "    c = profile_subset.values.T\n",
    "    s = ax.contourf(X, Y, c, 8, cmap=cm.YlOrBr, origin='lower', vmax=1, vmin=0)\n",
    "        # rotate and align the tick labels so they look better\n",
    "    ax.set_ylabel('Depth (m)')\n",
    "    [y1,y2]=ax.get_ylim()\n",
    "    for ihr in range(3,24,3):\n",
    "        tt=t[0]+np.timedelta64(ihr,'h')\n",
    "        ax.plot([tt.data,tt.data],[y1,y2])\n",
    "    ax.text(t[0].data+np.timedelta64(10,'m'),-0.4,'B',size=12,weight='bold')\n",
    "    fig.autofmt_xdate()\n",
    "    m = plt.cm.ScalarMappable(cmap=cm.YlOrBr)\n",
    "    m.set_array(c)\n",
    "    m.set_clim(0., 1.)\n",
    "    colorbar_ax = fig.add_axes([0.905, 0.55, 0.015, 0.2])\n",
    "    fig.colorbar(m, cax=colorbar_ax,orientation='vertical',label='$\\Delta$SST [K]',ticks=[0,.25,.5,.75,1])\n",
    "    #fig.set_label('$\\Delta$SST [K]')\n",
    "    ax = fig.add_axes([0.1, 0.3, 0.8, 0.2])\n",
    "    ax0 = ax.pcolormesh(subset_adcp.tlmt,subset_adcp.cell_depth,subset_adcp_masked.vel_horiz.T,vmin=0,vmax=.5,cmap='seismic')\n",
    "    #ax0.set_xlabel('Time')\n",
    "    [y1,y2]=ax.get_ylim()\n",
    "    for ihr in range(3,24,3):\n",
    "        tt=t[0]+np.timedelta64(ihr,'h')\n",
    "        ax.plot([tt.data,tt.data],[y1,y2])\n",
    "    ax.set_ylabel('Depth (m)')\n",
    "    ax.invert_yaxis()    \n",
    "    ax.text(t[0].data+np.timedelta64(10,'m'),90,'C',size=12,weight='bold')\n",
    "    \n",
    "    colorbar_ax = fig.add_axes([0.905, 0.3, 0.015, 0.2])\n",
    "    fig.colorbar(ax0, cax=colorbar_ax,orientation='vertical',label='Cur. Velocity (ms$^{-1}$)',ticks=[0,.5,1])\n",
    "\n",
    "    ax2 = fig.add_axes([.1, 0.05, 0.4, 0.2])\n",
    "    ilen_profile = profile_subset.shape[0]\n",
    "    colors = plt.cm.jet(np.linspace(0,1,ilen_profile))\n",
    "    for i in range(0,ilen_profile):\n",
    "        ax0= plt.plot(profile_subset[i,:].values,z,color=colors[i])\n",
    "    norm = mpl.colors.Normalize(vmin=6,vmax=30)\n",
    "    cmap = colorbrewer.get_map('Spectral', 'diverging', 11, reverse=True).mpl_colormap\n",
    "    sm = plt.cm.ScalarMappable(cmap=cmap, norm=norm)\n",
    "    sm.set_array([])\n",
    "    plt.colorbar(sm,boundaries=np.arange(6,30,.1),label='Hour')\n",
    "   \n",
    "#    for i in range(0,len(profile_subset.time)):\n",
    "#        ax2.plot(profile_subset[i,:].values,z,'-',markersize=.5,lw=1)\n",
    "#        ax2.plot(profile_subset[i,:].values,z,'.',markersize=.5,lw=1)\n",
    "    ax2.set_xlabel('$\\Delta$SST (K)')\n",
    "    ax2.set_ylabel('Depth (m)')\n",
    "    [x1,x2] = ax2.get_xlim()\n",
    "    [y1,y2] = ax2.get_xlim()\n",
    "    ax2.text(x1,-.35,'D',size=12,weight='bold')\n",
    "\n",
    "    \n",
    "    tem2=(subset_usv.time[0].values-np.datetime64('1970-01-01'))/np.timedelta64(1,'s')\n",
    "    tem3=datetime.utcfromtimestamp(tem2)\n",
    "    day_of_year = tem3.timetuple().tm_yday\n",
    "   \n",
    "    mn = xr.open_dataset(adir+'sat_data/virs_sst/all_sst'+str(day_of_year)+'.nc')\n",
    "    mn.close()\n",
    "#calculate overpass with most values\n",
    "    subset_mn = mn.sel(lon=slice(subset_usv.longitude.min().data,subset_usv.longitude.max().data),\n",
    "                       lat=slice(subset_usv.latitude.max().data,subset_usv.latitude.min().data))\n",
    "    s=np.zeros(subset_mn.dims['time'])\n",
    "    for i in range(0,subset_mn.dims['time']):\n",
    "        s[i]=(subset_mn.sea_surface_temperature[i,:,:]/subset_mn.sea_surface_temperature[i,:,:]).sum()\n",
    "    imax = np.argmax(s)\n",
    "    subset_mn = mn.sel(lon=slice(subset_usv.longitude.min().data-.1,subset_usv.longitude.max().data+.1),\n",
    "                       lat=slice(subset_usv.latitude.max().data+.1,subset_usv.latitude.min().data-.1))\n",
    "    subset_mn = subset_mn.isel(time=imax)\n",
    "#end calculate overpass with most values of sst\n",
    "\n",
    "#    subset_mn = mn.sel(lon=slice(subset_usv.longitude.min().data-1,subset_usv.longitude.max().data+1),\n",
    "#                       lat=slice(subset_usv.latitude.max().data+1,subset_usv.latitude.min().data-1))\n",
    "    minsst = subset_mn.sea_surface_temperature.min().data-273.15\n",
    "    maxsst = subset_mn.sea_surface_temperature.max().data-273.15\n",
    "    print(minsst,maxsst)\n",
    "    #blue marble image of ssts for cruise\n",
    "    ax4 = fig.add_axes([.58, 0.05, 0.29, 0.2])    \n",
    "    cmap = colorbrewer.get_map('Spectral', 'diverging', 11, reverse=True).mpl_colormap\n",
    "    ax0 = ax4.pcolormesh(subset_mn.lon,subset_mn.lat, subset_mn.sea_surface_temperature - 273.15, cmap=cmap, vmin = minsst, vmax = maxsst)\n",
    "    ax4.plot(subset_usv.longitude.data,subset_usv.latitude.data,'w', linewidth=3)\n",
    "    ax4.plot(subset_usv.longitude.data,subset_usv.latitude.data,color='magenta', linewidth=1)\n",
    "    for ihr in range(3,24,3):\n",
    "        tt=subset_usv.time[0]+np.timedelta64(ihr,'h')\n",
    "        sub_lons = subset_usv.longitude.sel(time=tt,method='nearest')\n",
    "        sub_lats = subset_usv.latitude.sel(time=tt,method='nearest')\n",
    "        ax4.plot(sub_lons,sub_lats,'o')\n",
    "    #print('sub',sub_lons[0].data,sub_lats[0].data)\n",
    "    ax4.text(subset_mn.lon.max().data-.04,subset_mn.lat.min().data+.02,'E',size=12,weight='bold')\n",
    "    colorbar_ax = fig.add_axes([0.9, 0.05, 0.015, 0.2])\n",
    "    fig.colorbar(ax0, cax=colorbar_ax,orientation='vertical',label='SST (C)')\n",
    "#    ax0.colorbar(location=\"bottom\",label='SST [C]') # draw colorbar\n",
    "    adir_fig='F:/data/cruise_data/saildrone/baja-2018/figs/diurnal/'\n",
    "    fig.savefig(adir_fig+'temp_loggers__diurnal_adcp'+str(day_of_year)+'3.png')\n",
    "    \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sub_lons"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#old sst plotting code that used basemap which was nice because it has nice land\n",
    "#but lacks lat / lon so the zoomed in images didn't tell a lot without location\n",
    "m = Basemap(projection='merc', lat_0 = 27, lon_0 = -120,    resolution = 'l', area_thresh = 0.1,    \n",
    "                                  llcrnrlon=min(subset_usv.longitude.data)-.05, llcrnrlat=min(subset_usv.latitude.data)-.05,\n",
    "                urcrnrlon=max(subset_usv.longitude.data)+.05, urcrnrlat=max(subset_usv.latitude.data)+.05)\n",
    "    m.bluemarble()\n",
    "    #    m.fillcontinents(color='grey',lake_color='white')\n",
    "    lon_grid, lat_grid = np.meshgrid(subset_mn.lon, subset_mn.lat)\n",
    "    x,y = m(lon_grid,lat_grid)\n",
    "    cs = m.pcolormesh(x,y, subset_mn.sea_surface_temperature - 273.15, cmap=cmap, vmin = minsst, vmax = maxsst)\n",
    "    x,y = m(subset_usv.longitude.data,subset_usv.latitude.data)\n",
    "    m.plot(x,y,'w', linewidth=3)\n",
    "    m.plot(x,y,color='magenta', linewidth=1)\n",
    "    for ihr in range(3,24,3):\n",
    "        tt=subset_usv.time[0]+np.timedelta64(ihr,'h')\n",
    "        sub_lons = subset_usv.longitude.sel(time=tt,method='nearest')\n",
    "        sub_lats = subset_usv.latitude.sel(time=tt,method='nearest')\n",
    "        x,y = m(sub_lons,sub_lats)\n",
    "        m.plot(x,y,'o')\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}