#!/usr/bin/env python # coding: utf-8 # ### Oxygen Evals in PugetSound (data from Dakota - UW) # In[16]: import xarray as xr import netCDF4 as nc import numpy as np import pandas as pd import matplotlib.pyplot as plt from pathlib import Path import cmocean.cm as cm import copy from salishsea_tools import visualisations as vis from salishsea_tools import viz_tools import datetime as dt import arrow import copy import math import matplotlib.dates as mdates import gsw from mpl_toolkits.basemap import Basemap import matplotlib.patches as mpatches import scipy.stats as stat # In[17]: mesh = nc.Dataset('/ocean/atall/MOAD/grid/mesh_mask202108.nc') bathy = nc.Dataset('/ocean/atall/MOAD/grid/bathymetry_202108.nc') meshb = nc.Dataset('/ocean/atall/MOAD/grid/mesh_mask_202310b.nc') bathyb = nc.Dataset('/ocean/atall/MOAD/grid/bathymetry_202310b.nc') grid = xr.open_dataset('/ocean/atall/MOAD/grid/bathymetry_202310b.nc', mask_and_scale=False) depthb = meshb.variables['gdept_0'][:] with xr.open_dataset('/ocean/atall/MOAD/grid/mesh_mask_202310b.nc') as mesh: tmask = mesh.tmask mbathy = mesh.mbathy long = mesh.nav_lon latg = mesh.nav_lat grid_dir = Path("/ocean/atall/MOAD/grid/") grid_map = Path("grid_from_lat_lon_mask999.nc") grid_lons_lats = xr.open_dataset(grid_dir / grid_map) # In[18]: ctdd = pd.read_pickle('/ocean/atall/MOAD/Obs/ecology_fromDakota_UW/ctd_1999_2016_fixed.p') ctd17 = pd.read_pickle('/ocean/atall/MOAD/Obs/ecology_fromDakota_UW/ctd_2017_fixed.p') ctd18 = pd.read_pickle('/ocean/atall/MOAD/Obs/ecology_fromDakota_UW/ctd_2018_fixed.p') ctd19 = pd.read_pickle('/ocean/atall/MOAD/Obs/ecology_fromDakota_UW/ctd_2019_fixed.p') bot = pd.read_pickle('/ocean/atall/MOAD/Obs/ecology_fromDakota_UW/bottle_2006_2017_fixed.p') sta_df = pd.read_pickle('/ocean/atall/MOAD/Obs/ecology_fromDakota_UW/sta_df.p') sta_df # In[19]: ctd17 # ### https://www.ices.dk/data/tools/Pages/Unit-conversions.aspx # O2: 1 mg/l = 0.7 ml/l = 0.7*44.661 umol/l # In[20]: # Put coordonates in dataframe --- matching with model Lat=dict() Lon=dict() for alpha in range(0,5): if alpha==0: ctd=ctdd filename = 'ctd_1999-2016' ctd['DO_adjusted']=ctd.DO_adjusted*0.7*44.661 elif alpha==1: ctd=ctd17 filename = 'ctd_2017' ctd['DO_raw']=ctd.DO_raw*0.7*44.661 elif alpha==2: ctd=ctd18 filename = 'ctd_2018' ctd['DO_adjusted']=ctd.DO_adjusted*0.7*44.661 elif alpha==3: ctd=ctd19 filename = 'ctd_2019' ctd['DO_raw']=ctd.DO_raw*0.7*44.661 else: ctd=bot filename = 'bot_2006-2017' for i in range(0,len(ctd.Station)): if ctd.Station[i]=='ADM001': Lat[i] = 48.0298133333 Lon[i] = -122.6179333333 elif ctd.Station[i]=='ADM002': Lat[i] = 48.1873183333 Lon[i] = -122.84295 elif ctd.Station[i]=='ADM003': Lat[i] = 47.8789833333 Lon[i] = -122.483195 elif ctd.Station[i]=='BLL009': Lat[i] = 48.68594 Lon[i] = -122.5996183333 elif ctd.Station[i]=='BUD005': Lat[i] = 47.09204 Lon[i] = -122.9181966667 elif ctd.Station[i]=='CMB003': Lat[i] = 47.2903766667 Lon[i] = -122.4501233333 elif ctd.Station[i]=='CRR001': Lat[i] = 47.276485 Lon[i] = -122.709575 elif ctd.Station[i]=='CSE001': Lat[i] = 47.26454 Lon[i] = -122.844305 elif ctd.Station[i]=='DNA001': Lat[i] = 47.1614833333 Lon[i] = -122.871805 elif ctd.Station[i]=='EAP001': Lat[i] = 47.4170433333 Lon[i] = -122.3804016667 elif ctd.Station[i]=='ELB015': Lat[i] = 47.5964866667 Lon[i] = -122.3695716667 elif ctd.Station[i]=='GOR001': Lat[i] = 47.1831516667 Lon[i] = -122.6345716667 elif ctd.Station[i]=='GRG002': Lat[i] = 48.80816 Lon[i] = -122.9540766667 elif ctd.Station[i]=='GYS004': Lat[i] = 46.97787 Lon[i] = -123.78461 elif ctd.Station[i]=='GYS008': Lat[i] = 46.9373133333 Lon[i] = -123.9132233333 elif ctd.Station[i]=='GYS016': Lat[i] = 46.9534216667 Lon[i] = -124.09295 elif ctd.Station[i]=='HCB003': Lat[i] = 47.53787 Lon[i] = -123.0096 elif ctd.Station[i]=='HCB004': Lat[i] = 47.356205 Lon[i] = -123.0248733333 elif ctd.Station[i]=='HCB007': Lat[i] = 47.3981483333 Lon[i] = -122.9295916667 elif ctd.Station[i]=='HCB010': Lat[i] = 47.67 Lon[i] = -122.82 elif ctd.Station[i]=='NSQ002': Lat[i] = 47.1673166667 Lon[i] = -122.78819 elif ctd.Station[i]=='OAK004': Lat[i] = 47.2134266667 Lon[i] = -123.07765 elif ctd.Station[i]=='PSB003': Lat[i] = 47.6598183333 Lon[i] = -122.4429083333 elif ctd.Station[i]=='PSS019': Lat[i] = 48.0109266667 Lon[i] = -122.30125 elif ctd.Station[i]=='PTH005': Lat[i] = 48.0831483333 Lon[i] = -122.7646116667 elif ctd.Station[i]=='RSR837': Lat[i] = 48.6164916667 Lon[i] = -122.7629583333 elif ctd.Station[i]=='SAR003': Lat[i] = 48.107595 Lon[i] = -122.4915416667 elif ctd.Station[i]=='SIN001': Lat[i] = 47.5492616667 Lon[i] = -122.6434716667 elif ctd.Station[i]=='SJF000': Lat[i] = 48.4166666667 Lon[i] = -123.025 elif ctd.Station[i]=='SJF001': Lat[i] = 48.3333333333 Lon[i] = -123.025 elif ctd.Station[i]=='SJF002': Lat[i] = 48.25 Lon[i] = -123.025 elif ctd.Station[i]=='SKG003': Lat[i] = 48.2964883333 Lon[i] = -122.489605 elif ctd.Station[i]=='WPA001': Lat[i] = 46.6873216667 Lon[i] = -123.7498816667 elif ctd.Station[i]=='WPA003': Lat[i] = 46.7039866667 Lon[i] = -123.837385 elif ctd.Station[i]=='WPA004': Lat[i] = 46.6867633333 Lon[i] = -123.9735 elif ctd.Station[i]=='WPA006': Lat[i] = 46.5453766667 Lon[i] = -123.9801616667 elif ctd.Station[i]=='WPA007': Lat[i] = 46.453155 Lon[i] = -124.0096033333 elif ctd.Station[i]=='WPA008': Lat[i] = 46.463155 Lon[i] = -123.9412683333 elif ctd.Station[i]=='WPA113': Lat[i] = 46.644 Lon[i] = -123.993 ctd['Lat']=Lat ctd['Lon']=Lon ctd.to_pickle('/ocean/atall/MOAD/Obs/ecology_fromDakota_UW/fixedTall/{}.p'.format(f'{filename}')) # In[46]: # matched data ctd202111_12 = pd.read_csv('/ocean/atall/MOAD/ObsModel/202111/ObsModel_202111_ctd_from_UW_20120101_20121231.csv') ctd202410_12 = pd.read_csv('/ocean/atall/MOAD/ObsModel/202410/ObsModel_202410_ctd_from_UW_20120101_20121231.csv') ctd202111_13 = pd.read_csv('/ocean/atall/MOAD/ObsModel/202111/ObsModel_202111_ctd_from_UW_20130101_20131231.csv') ctd202410_13 = pd.read_csv('/ocean/atall/MOAD/ObsModel/202410/ObsModel_202410_ctd_from_UW_20130101_20131231.csv') ctd202111_14 = pd.read_csv('/ocean/atall/MOAD/ObsModel/202111/ObsModel_202111_ctd_from_UW_20140101_20141231.csv') ctd202410_14 = pd.read_csv('/ocean/atall/MOAD/ObsModel/202410/ObsModel_202410_ctd_from_UW_20140101_20141231.csv') ctd202111 = pd.concat([ctd202111_12, ctd202111_13, ctd202111_14 ], ignore_index=True) ctd202410 = pd.concat([ctd202410_12, ctd202410_13, ctd202410_14 ], ignore_index=True) # In[47]: fig, ax = plt.subplots(1,1,figsize=(5, 9)) lon1,lon2,lat1,lat2 = (-123.2,-122.2,47.05,48.13) ax.contourf(long, latg, mbathy[0,:,:], linewidths=1, levels=[-0.01, 0.01], colors='whitesmoke') ax.contour(long, latg, mbathy[0,:,:], linewidths=1, levels=[-0.01, 0.01], colors='dimgray') ax.set_ylabel('Latitude') ax.set_xlabel('Longitude') ax.set_ylim([47,49]) ax.set_xlim([-123.5,-122]) # plot the location of observations ax.scatter(ctd202111.Lon, ctd202111.Lat, s=5, label="ctd_1999-2016") ax.legend() left, bottom, width, height = (lon2, lat1, lon1-lon2, lat2-lat1) rect=mpatches.Rectangle((left,bottom),width,height, fill=False, #alpha=0.1 color="purple", linewidth=2, label="Puget Sound") plt.gca().add_patch(rect) # In[48]: df_21 = ctd202111[ ctd202111['Lon'].between(lon1, lon2) & ctd202111['Lat'].between(lat1, lat2) ] df_24 = ctd202410[ ctd202410['Lon'].between(lon1, lon2) & ctd202410['Lat'].between(lat1, lat2) ] # In[49]: def profiles(tracer,colour,ax): if tracer == 'Salinity': t_obs = 'Salinity' t_mod = 'mod_vosaline' unit = 'g/kg' unity ='meter' elif tracer == 'DO': t_obs = 'Oxygen_Dissolved' t_mod = 'mod_dissolved_oxygen' unit = 'uM' unity ='meter' avg_obs, binsa, _ = stat.binned_statistic(-df_21['Z'][(np.isfinite(df_21[t_obs]))],df_21[t_obs][(np.isfinite(df_21[t_obs]))],statistic='mean',bins=100) avg_21, bins, _ = stat.binned_statistic(-df_21['Z'][(np.isfinite(df_21[t_mod]))],df_21[t_mod][(np.isfinite(df_21[t_mod]))],statistic='mean',bins=100) avg_24, binsa, _ = stat.binned_statistic(-df_24['Z'][(np.isfinite(df_24[t_mod]))],df_24[t_mod][(np.isfinite(df_24[t_mod]))],statistic='mean',bins=100) ax.plot(avg_obs, binsa[:-1], lw=2,label='obs') ax.plot(avg_21, bins[:-1], lw=2,label='202111') ax.plot(avg_24, binsa[:-1], lw=2,label='202410') title = tracer #ax.set_title(title) ax.set_xlabel(unit) ax.set_ylabel(unity) # In[51]: print('2012-2014') i, j = (0, 1) fig, ax = plt.subplots(1, 2, figsize=(8, 4)) title = 'Salinity, PugetSound' ax[i].set_title(title,fontsize=12) title = 'Oxygen, PugetSound' ax[j].set_title(title,fontsize=12) ax[i].plot(df_24.Salinity, -df_24.Z, '.',label='obs') ax[i].plot(df_21.mod_vosaline, -df_21.Z, '.', label='202111') ax[i].plot(df_24.mod_vosaline, -df_24.Z, '.', label='202410') ax[i].legend() ax[j].plot(df_24.Oxygen_Dissolved, -df_24.Z, '.',label='obs') ax[j].plot(df_21.mod_dissolved_oxygen, -df_21.Z, '.', label='202111') ax[j].plot(df_24.mod_dissolved_oxygen, -df_24.Z, '.', label='202410') fig, ax = plt.subplots(1, 2, figsize=(8,4)) title = 'Salinity, PugetSound' ax[i].set_title(title,fontsize=12) title = 'Oxygen, PugetSound' ax[j].set_title(title,fontsize=12) # plot profiles profiles('Salinity','k',ax[0]) profiles('DO','k',ax[1]) ax[0].legend() plt.tight_layout() # In[ ]: