diff --git a/biascorr/notebooks/plot_series.ipynb b/biascorr/notebooks/plot_series.ipynb new file mode 100644 index 0000000..804d829 --- /dev/null +++ b/biascorr/notebooks/plot_series.ipynb @@ -0,0 +1,273 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import xarray as xr\n", + "import matplotlib as mpl\n", + "from matplotlib.gridspec import GridSpec, GridSpecFromSubplotSpec\n", + "import matplotlib.pyplot as plt\n", + "import numpy as np\n", + "import cmocean as cmo\n", + "\n", + "import sys\n", + "sys.path.append('../src')\n", + "\n", + "from run import Run\n", + "from region import Region\n", + "\n", + "import warnings\n", + "warnings.filterwarnings('ignore')\n", + "\n", + "%load_ext autoreload\n", + "%autoreload 2" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Plot regions\n", + "\n", + "fig,ax = plt.subplots(1,4,figsize=(16,4))\n", + "\n", + "for m,mod in enumerate(['CESM2','EC-Earth3','UKESM1-0-LL','WOA23']):\n", + " self = Run(mod,'Ross','historical',1995,1996)\n", + " ax[m].pcolormesh(self.lon,self.lat,np.nansum(self.V,axis=0))\n", + " #ax[m].pcolormesh(np.nansum(self.V,axis=0))\n", + " ax[m].set_title(mod)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "def get_series(region,years):\n", + " Traw_av = np.zeros((len(years),len(region.prd.depth)))\n", + " Tcor_av = np.zeros((len(years),len(region.ref.depth)))\n", + " dTraw_av = np.zeros((len(years),len(region.prd.depth)))\n", + " dTcor_av = np.zeros((len(years),len(region.ref.depth)))\n", + "\n", + " for y,year in enumerate(years):\n", + " if year>2014:\n", + " run = 'ssp585'\n", + " else:\n", + " run = 'historical'\n", + " region.get_future(run,year,year+1)\n", + "\n", + " dTraw_av[y,:] = region.calc_horav(region.fut.dTraw,region.prd)\n", + " Traw_av[y,:] = region.calc_horav(region.fut.T,region.fut)\n", + " dTcor_av[y,:] = region.calc_horav(region.fut.dTcor,region.ref)\n", + " Tcor_av[y,:] = region.calc_horav(region.ref.T+region.fut.dTcor,region.ref)\n", + " return Traw_av,Tcor_av,dTraw_av,dTcor_av\n", + "\n", + "def add_Tvar(fig,gs,n,m,zvar,var,years,title='',ylabel=''):\n", + " ax = fig.add_subplot(gs[n,m])\n", + " im = ax.pcolormesh(years,zvar,var.T,cmap='cmo.thermal',vmin=-2,vmax=3)\n", + " ax.set_ylim([-1500,0])\n", + " if m>0:\n", + " ax.set_yticklabels([])\n", + " ax.set_ylabel(ylabel)\n", + " if n==0:\n", + " ax.set_title(title)\n", + " if n<2:\n", + " ax.set_xticklabels([])\n", + " return im\n", + "\n", + "def add_dTvar(fig,gs,n,m,zvar,var,years,title=''):\n", + " ax = fig.add_subplot(gs[n,m])\n", + " im = ax.pcolormesh(years,zvar,var.T,cmap='cmo.balance',vmin=-4,vmax=4)\n", + " ax.set_ylim([-1500,0])\n", + " ax.set_yticklabels([])\n", + " if n==0:\n", + " ax.set_title(title)\n", + " if n<2:\n", + " ax.set_xticklabels([])\n", + " return im\n", + "\n", + "def add_row(fig,gs,n,regionname,years,model,refmodel):\n", + " region = Region(regionname,model,refmodel,k0=50,k1=2000)\n", + " Traw_av,Tcor_av,dTraw_av,dTcor_av = get_series(region,years)\n", + "\n", + " im = add_Tvar(fig,gs,n,0,region.prd.depth,Traw_av,years,title='Raw model',ylabel=region.model)\n", + " im = add_Tvar(fig,gs,n,1,region.ref.depth,Tcor_av,years,title=f'Corrected (ref: {region.refmodel})')\n", + "\n", + " if n==0:\n", + " ax = fig.add_subplot(gs[:,2])\n", + " cb = plt.colorbar(im,cax=ax,extend='both')\n", + "\n", + " im = add_dTvar(fig,gs,n,4,region.prd.depth,dTraw_av,years,title='Raw model')\n", + " im = add_dTvar(fig,gs,n,5,region.ref.depth,dTcor_av,years,title=f'Corrected')\n", + "\n", + " if n==0:\n", + " ax = fig.add_subplot(gs[:,6])\n", + " cb = plt.colorbar(im,cax=ax,extend='both')\n", + "\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "#years = np.arange(1850,2110,10)\n", + "years = np.arange(1850,2101,10)\n", + "refmodel = 'CESM2'\n", + "models = ['CESM2','UKESM1-0-LL','EC-Earth3']\n", + "#models = ['EC-Earth3']\n", + "regionname = 'Amundsen'\n", + "\n", + "fig = plt.figure(figsize=(15,5))#,constrained_layout=True)\n", + "gs = GridSpec(3,8, figure=fig,\n", + " width_ratios=[1,1,.05,.05,1,1,.05,.05],\n", + " height_ratios=[1,1,1],wspace=.05,hspace=.06)\n", + "\n", + "for n,model in enumerate(models):\n", + " add_row(fig,gs,n,regionname,years,model,refmodel)\n", + "\n", + "fig.suptitle(regionname)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "years = np.arange(1850,2110,10)\n", + "#years = np.arange(1850,2101,1)\n", + "refmodel = 'WOA23'\n", + "models = ['CESM2','UKESM1-0-LL','EC-Earth3']\n", + "#models = ['EC-Earth3']\n", + "regionname = 'Ross'\n", + "\n", + "fig = plt.figure(figsize=(15,5))#,constrained_layout=True)\n", + "gs = GridSpec(3,8, figure=fig,\n", + " width_ratios=[1,1,.05,.05,1,1,.05,.05],\n", + " height_ratios=[1,1,1],wspace=.05,hspace=.06)\n", + "\n", + "for n,model in enumerate(models):\n", + " add_row(fig,gs,n,regionname,years,model,refmodel)\n", + "\n", + "fig.suptitle(regionname)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "#years = np.arange(1850,2110,10)\n", + "years = np.arange(1850,2101,1)\n", + "refmodel = 'WOA23'\n", + "models = ['CESM2','UKESM1-0-LL','EC-Earth3']\n", + "#models = ['EC-Earth3']\n", + "regionname = 'Weddell'\n", + "\n", + "fig = plt.figure(figsize=(15,5))#,constrained_layout=True)\n", + "gs = GridSpec(3,8, figure=fig,\n", + " width_ratios=[1,1,.05,.05,1,1,.05,.05],\n", + " height_ratios=[1,1,1],wspace=.05,hspace=.06)\n", + "\n", + "for n,model in enumerate(models):\n", + " add_row(fig,gs,n,regionname,years,model,refmodel)\n", + "\n", + "fig.suptitle(regionname)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "#years = np.arange(1850,2110,10)\n", + "years = np.arange(1850,2101,1)\n", + "refmodel = 'WOA23'\n", + "models = ['CESM2','UKESM1-0-LL','EC-Earth3']\n", + "#models = ['EC-Earth3']\n", + "regionname = 'Totten'\n", + "\n", + "fig = plt.figure(figsize=(15,5))#,constrained_layout=True)\n", + "gs = GridSpec(3,8, figure=fig,\n", + " width_ratios=[1,1,.05,.05,1,1,.05,.05],\n", + " height_ratios=[1,1,1],wspace=.05,hspace=.06)\n", + "\n", + "for n,model in enumerate(models):\n", + " add_row(fig,gs,n,regionname,years,model,refmodel)\n", + "\n", + "fig.suptitle(regionname)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "region = Region('Amundsen','EC-Earth3',refmodel='CESM2',k0=50,k1=2000)\n", + "region.get_future('ssp585',2050,2060)\n", + "region.plot_delta()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "region = Region('Totten','CESM2',refmodel='WOA23',k0=50,k1=2000)\n", + "region.get_future('ssp585',2050,2060)\n", + "region.plot_delta()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "region = Region('Ross','EC-Earth3',refmodel='WOA23',k0=50,k1=2000)\n", + "region.get_future('ssp585',2050,2060)\n", + "region.plot_delta()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3 (ipykernel)", + "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.8.12" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +} diff --git a/biascorr/src/region.py b/biascorr/src/region.py new file mode 100644 index 0000000..95662fe --- /dev/null +++ b/biascorr/src/region.py @@ -0,0 +1,252 @@ +import xarray as xr +import matplotlib as mpl +import matplotlib.pyplot as plt +import numpy as np +import cmocean as cmo + +from run import Run +from run import dpm + +# Region class with functions for processing and plotting +class Region: + def __init__(self, region,model,refmodel='UKESM1-0-LL',k0=0,k1=6000): + # Get volume CMIP + self.region = region + self.model = model # Reference model considered 'truth' + self.refmodel = refmodel # Model to be bias corrected + self.k0 = k0 # Start depth (top) + self.k1 = k1 # End depth (bottom) + + # Read historical period of both refmodel and model to get started + self.get_historical() + + def get_historical(self): + # Get data for refmodel over historical period + self.ref = Run(self.refmodel,self.region,'historical',1995,2015,k0=self.k0,k1=self.k1) + + # Get data for model over historical period + self.prd = Run(self.model,self.region,'historical',1995,2015,k0=self.k0,k1=self.k1) + + # Define the bias in the model with respect to the reference model + self.get_bias() + + def get_future(self,run,y0,y1): + # Get data for model over some future period y0 - y1 + self.fut = Run(self.model,self.region,run,y0,y1,k0=self.k0,k1=self.k1) + + # Compute change in T and S between future and historical in T,S space of model + # Example: between y0 and y1, the water mass that historically had a temperature of X degrees and a salinity of Y psu + # has warmed by an amount of Z degrees. + self.get_delta() + + # Construct deltaT and deltaS on the grid of the reference model for this future period + self.get_anom() + + def plot_region_cmip(self,ii,jj): + plt.pcolormesh(self.lonc,self.latc,np.nansum(self.thkcello,axis=0),cmap='Blues',norm=mpl.colors.LogNorm(vmin=10,vmax=1e4)) + plt.scatter(self.lonc[jj,ii],self.latc[jj,ii],50,c='tab:red') + + def plot_volumes(self): + fig,ax = plt.subplots(1,3,figsize=(3*5,6),sharex=True,sharey=True) + + ax[0].pcolormesh(self.ref.Sb,self.ref.Tb,self.prd.Vb.T,norm=mpl.colors.LogNorm(vmin=1e9,vmax=1e13),cmap='Greys') + plt.colorbar(im,ax=ax[0],orientation='horizontal') + ax[0].set_title(f'Reference ({self.ref.model} {self.ref.y0}-{self.ref.y1})') + ax[0].set_ylabel(self.region) + + im = ax[1].pcolormesh(self.prd.Sb,self.prd.Tb,self.prd.Vb.T,norm=mpl.colors.LogNorm(vmin=1e9,vmax=1e13)) + plt.colorbar(im,ax=ax[1],orientation='horizontal') + ax[1].set_title(f'Present day {self.model} ({self.prd.y0}-{self.prd.y1})') + #ax[1].plot(self.prd.S[:,jj,ii],self.prd.T[:,jj,ii],c='tab:red') + + im = ax[2].pcolormesh(self.fut.Sb,self.fut.Tb,self.fut.Vb.T,norm=mpl.colors.LogNorm(vmin=1e9,vmax=1e13),cmap='Oranges') + plt.colorbar(im,ax=ax[2],orientation='horizontal') + ax[2].set_title(f'{self.fut.run} ({self.fut.y0}-{self.fut.y1})') + #ax[2].plot(self.fut.S[:,jj,ii],self.fut.T[:,jj,ii],c='tab:red') + + def plot_delta(self): + fig,ax = plt.subplots(1,5,figsize=(5*3.5,4.5),sharex=True,sharey=True) + + im = ax[0].pcolormesh(self.ref.Sb,self.ref.Tb,self.ref.Vb.T,norm=mpl.colors.LogNorm(vmin=1e9,vmax=1e13)) + plt.colorbar(im,ax=ax[0],orientation='horizontal') + ax[0].set_title(f'Reference ({self.ref.model} {self.ref.y0}-{self.ref.y1})') + ax[0].set_ylabel(self.region) + + im = ax[1].pcolormesh(self.prd.Sb,self.prd.Tb,self.prd.Vb.T,norm=mpl.colors.LogNorm(vmin=1e9,vmax=1e13)) + plt.colorbar(im,ax=ax[1],orientation='horizontal') + ax[1].set_title(f'present day {self.model} ({self.prd.y0}-{self.prd.y1})') + #ax[0].plot(self.prd.S[:,jj,ii],self.prd.T[:,jj,ii],c='tab:red') + + #im = ax[2].pcolormesh(self.prd.Sc,self.prd.Tc,self.fut.Vb.T,norm=mpl.colors.LogNorm(vmin=1e9,vmax=1e13),cmap='Greys') + im = ax[2].pcolormesh(self.prd.Sb,self.prd.Tb,self.fut.deltaT.T,cmap='cmo.balance',vmin=-3,vmax=3) + plt.colorbar(im,ax=ax[2],orientation='horizontal') + ax[2].set_title(f'{self.fut.run} ({self.fut.y0}-{self.fut.y1})') + + im = ax[3].pcolormesh(self.prd.Sb,self.prd.Tb,self.fut.deltaTf.T,cmap='cmo.balance',vmin=-3,vmax=3) + + plt.colorbar(im,ax=ax[3],orientation='horizontal') + ax[3].set_title(f'Filled') + #ax[1].plot(self.fut.S[:,jj,ii],self.fut.T[:,jj,ii],c='tab:red') + + im = ax[4].pcolormesh(self.ref.Sb,self.ref.Tb,self.fut.dTcorb.T,cmap='cmo.balance',vmin=-3,vmax=3) + plt.colorbar(im,ax=ax[4],orientation='horizontal') + #ax[3].set_title(f'{self.fut.run} ({self.fut.y0}-{self.fut.y1})') + + def plot_profiles(self,vlim=-1): + fig,ax = plt.subplots(1,3,figsize=(10,5),sharey=True) + + Tav = np.average(np.where(np.isnan(self.ref.T[:vlim,:,:]),0,self.ref.T[:vlim,:,:]),axis=(1,2),weights=self.ref.V[:vlim,:,:]) + ax[1].plot(Tav,self.ref.depth[:vlim],label='optimised',c='tab:red',ls='--') + + Tav = np.average(np.where(np.isnan(self.prd.T[:vlim,:,:]),0,self.prd.T[:vlim,:,:]),axis=(1,2),weights=self.prd.V[:vlim,:,:]) + ax[0].plot(Tav,self.prd.depth[:vlim],label='raw',c='.5',ls='--') + + Tav = np.average(np.where(np.isnan(self.ref.T[:vlim,:,:]),0,self.fut.dTcor[:vlim,:,:]),axis=(1,2),weights=self.ref.V[:vlim,:,:]) + ax[2].plot(Tav,self.ref.depth[:vlim],c='tab:red') + + Tav = np.average(np.where(np.isnan(self.prd.T[:vlim,:,:]),0,self.fut.dTraw[:vlim,:,:]),axis=(1,2),weights=self.prd.V[:vlim,:,:]) + ax[2].plot(Tav,self.prd.depth[:vlim],c='.5') + + Tav = np.average(np.where(np.isnan(self.ref.T[:vlim,:,:]),0,(self.ref.T+self.fut.dTcor)[:vlim,:,:]),axis=(1,2),weights=self.ref.V[:vlim,:,:]) + ax[1].plot(Tav,self.ref.depth[:vlim],c='tab:red') + + Tav = np.average(np.where(np.isnan(self.prd.T[:vlim,:,:]),0,self.fut.T[:vlim,:,:]),axis=(1,2),weights=self.prd.V[:vlim,:,:]) + ax[0].plot(Tav,self.prd.depth[:vlim],c='.5') + + ax[0].set_title('Raw model') + ax[1].set_title('Optimised') + ax[2].set_xlabel('Warming T') + #ax[0].legend() + + for Ax in ax: + Ax.set_ylim([-1500,0]) + #Ax.invert_yaxis() + + def get_bias(self,Tperc=.99,Sperc=.99): + + # S-bias constrain 95th percentile + A,B = np.histogram((self.ref.S).flatten()[self.ref.V.flatten()>0],weights=self.ref.V.flatten()[self.ref.V.flatten()>0],bins=1000,density=True) + AA = np.cumsum(A)/np.cumsum(A)[-1] + aa = np.argmin((AA-Sperc)**2) + self.ref.S95 = B[aa] + + A,B = np.histogram((self.prd.S).flatten()[self.prd.V.flatten()>0],weights=self.prd.V.flatten()[self.prd.V.flatten()>0],bins=1000,density=True) + AA = np.cumsum(A)/np.cumsum(A)[-1] + aa = np.argmin((AA-Sperc)**2) + self.prd.S95 = B[aa] + + # S-bias + a = np.arange(.7,1.3,.01) + ref,dum = np.histogram(self.ref.S,bins=self.ref.Sb,weights=self.ref.V) + rmse = np.zeros((len(a))) + for A,AA in enumerate(a): + prd,dum = np.histogram(AA*(self.prd.S-self.prd.S95)+self.ref.S95,bins=self.ref.Sb,weights=self.prd.V) + rmse[A] = np.sum(np.where(prd==0,0,(np.log10(prd/np.sum(prd))-np.log10(ref/np.sum(ref))))**2)**.5 / np.sum(np.where(prd==0,0,1)) + #rmse[A] = np.sum(np.where(prd==0,0,(np.log10(prd/np.sum(prd))-np.log10(ref/np.sum(ref))))**2) / np.sum(np.where(prd==0,0,1))**.5 + + #print(AA,rmse[A]) + out = np.unravel_index(np.nanargmin(rmse, axis=None), rmse.shape) + self.dSa = a[out[0]] + self.prd.Sc = self.dSa*(self.prd.Sb-self.prd.S95)+self.ref.S95 #Corrected S + #self.fut.Sc = self.dSa*(self.fut.Sb-self.prd.S95)+self.ref.S95 #Corrected S + + # T-bias (scale 95th percentile of T-Tmin) + A,B = np.histogram((self.ref.T-self.ref.Tb[0]).flatten()[self.ref.V.flatten()>0],weights=self.ref.V.flatten()[self.ref.V.flatten()>0],bins=1000,density=True) + AA = np.cumsum(A)/np.cumsum(A)[-1] + aa = np.argmin((AA-Tperc)**2) + self.ref.TF95 = B[aa] + + A,B = np.histogram((self.prd.T-self.prd.Tb[0]).flatten()[self.prd.V.flatten()>0],weights=self.prd.V.flatten()[self.prd.V.flatten()>0],bins=1000,density=True) + AA = np.cumsum(A)/np.cumsum(A)[-1] + aa = np.argmin((AA-Tperc)**2) + self.prd.TF95 = B[aa] + self.dTa = self.ref.TF95/self.prd.TF95 + self.prd.Tc = self.dTa*(self.prd.Tb-self.prd.Tb[0])+self.ref.Tb[0] + #self.fut.Tc = self.dTa*(self.fut.Tb-self.fut.Tb[0])+self.ref.Tb[0] + + def get_delta(self): + + # Get difference in T and S between future and historical period + self.fut.dTraw = self.fut.T-self.prd.T + self.fut.dSraw = self.fut.S-self.prd.S + + # Compute binned deltaT and deltaS on historical S- and T- bins + VTdel,Sdel,Tdel = np.histogram2d(self.prd.S.flatten()[self.prd.V.flatten()>0],self.prd.T.flatten()[self.prd.V.flatten()>0],bins=100,weights=(self.prd.V*self.fut.dTraw).flatten()[self.prd.V.flatten()>0]) + self.fut.deltaT = VTdel/self.prd.Vb + VSdel,Sdel,Tdel = np.histogram2d(self.prd.S.flatten()[self.prd.V.flatten()>0],self.prd.T.flatten()[self.prd.V.flatten()>0],bins=100,weights=(self.prd.V*self.fut.dSraw).flatten()[self.prd.V.flatten()>0]) + self.fut.deltaS = VSdel/self.prd.Vb + + # Create filled binned deltaT and deltaS through extrapolation + self.fut.deltaTf = self.fill_delta(self.fut.deltaT) + self.fut.deltaSf = self.fill_delta(self.fut.deltaS) + + def fill_delta(self,deltavar): + '''Routine to fill deltaT and deltaS in full normalised T,S space''' + + newval = np.nan*np.ones((len(self.prd.Sc),len(self.prd.Tc))) + + newval[1:,1:] = deltavar + Nleft = np.sum(np.isnan(newval[1:,1:])) + while Nleft>0: + mask = np.where(np.isnan(newval),0,1) + newval2 = np.where(np.isnan(newval),0,newval) + AA = newval2*mask + AAm1 = np.roll(AA,-1,axis=0) + m1 = np.roll(mask,-1,axis=0) + AAp1 = np.roll(AA,1,axis=0) + p1 = np.roll(mask,1,axis=0) + newval3 = (np.roll(AAm1,-1,axis=1)+AAm1+np.roll(AAm1,1,axis=1) + np.roll(AA,-1,axis=1)+np.roll(AA,1,axis=1) + np.roll(AAp1,-1,axis=1)+AAp1+np.roll(AAp1,1,axis=1)) / (np.roll(m1,-1,axis=1)+m1+np.roll(m1,1,axis=1) + np.roll(mask,-1,axis=1)+np.roll(mask,1,axis=1) + np.roll(p1,-1,axis=1)+p1+np.roll(p1,1,axis=1)) + newval4 = np.where(mask,newval,newval3) + newval = newval4 + newval[0,:] = np.nan + newval[:,0] = np.nan + + Nleft = np.sum(np.isnan(newval[1:,1:])) + + return newval[1:,1:] + + def get_anom(self): + '''Get anomalies of reference T,S values''' + + self.fut.dTcor = np.zeros(self.ref.T.shape) + self.fut.dScor = np.zeros(self.ref.S.shape) + + for i in range(self.fut.dTcor.shape[0]): + for j in range(self.fut.dTcor.shape[1]): + for k in range(self.fut.dTcor.shape[2]): + if self.ref.V[i,j,k]==0: + continue + else: + isref = min(99,max(0,int(99*(self.ref.S[i,j,k]-self.prd.Sc[0])/(self.prd.Sc[-1]-self.prd.Sc[0])))) + jtref = min(99,max(0,int(99*(self.ref.T[i,j,k]-self.prd.Tc[0])/(self.prd.Tc[-1]-self.prd.Tc[0])))) + self.fut.dScor[i,j,k] = self.fut.deltaSf[isref,jtref] + self.fut.dTcor[i,j,k] = self.fut.deltaTf[isref,jtref] + + #Distributions of dS and dT, only for plotting purposes + VTdel,Sdel,Tdel = np.histogram2d(self.ref.S.flatten()[self.ref.V.flatten()>0],self.ref.T.flatten()[self.ref.V.flatten()>0],bins=100,range=[[self.ref.Sb[0],self.ref.Sb[-1]],[self.ref.Tb[0],self.ref.Tb[-1]]],weights=(self.ref.V*self.fut.dTcor).flatten()[self.ref.V.flatten()>0]) + self.fut.dTcorb = VTdel/self.ref.Vb + VSdel,Sdel,Tdel = np.histogram2d(self.ref.S.flatten()[self.ref.V.flatten()>0],self.ref.S.flatten()[self.ref.V.flatten()>0],bins=100,range=[[self.ref.Sb[0],self.ref.Sb[-1]],[self.ref.Tb[0],self.ref.Tb[-1]]],weights=(self.ref.V*self.fut.dScor).flatten()[self.ref.V.flatten()>0]) + self.fut.dScorb = VSdel/self.ref.Vb + + def get_profiles(self): + # Construct vertical profiles of reference T and S + + self.ref.Tz = np.nan*np.ones(len(self.ref.depth)) + self.ref.Sz = np.nan*np.ones(len(self.ref.depth)) + + kmax = np.where(np.sum(self.ref.V,axis=(1,2))==0)[0][0] + self.ref.Tz[:kmax] = np.average(np.where(np.isnan(self.ref.T[:kmax,:,:]),0,self.ref.T[:kmax,:,:]),axis=(1,2),weights=self.ref.V[:kmax,:,:]) + self.ref.Sz[:kmax] = np.average(np.where(np.isnan(self.ref.S[:kmax,:,:]),0,self.ref.S[:kmax,:,:]),axis=(1,2),weights=self.ref.V[:kmax,:,:]) + + def calc_horav(self,var,run): + '''Calculate horizontally averaged vertical profile''' + + try: + out = np.average(np.where(np.isnan(var),0,var),axis=(1,2),weights=run.V) + except: + out = np.nan*np.ones(var.shape[0]) + kmax = np.where(np.sum(run.V,axis=(1,2))==0)[0][0] + out[:kmax] = np.average(np.where(np.isnan(var[:kmax,:,:]),0,var[:kmax,:,:]),axis=(1,2),weights=run.V[:kmax,:,:]) + + return out + diff --git a/biascorr/src/run.py b/biascorr/src/run.py new file mode 100644 index 0000000..5b84cbe --- /dev/null +++ b/biascorr/src/run.py @@ -0,0 +1,244 @@ +import xarray as xr +import numpy as np + +#Days per month +def dpm(years): + # Construct an array of days per month for the required amount of years + return np.tile([31,28,31,30,31,30,31,31,30,31,30,31],years) + +# Run class, which can be model output or observations +class Run: + def __init__(self,model,region,run,y0,y1,k0=0,k1=6000): + self.model = model + self.region = region + self.run = run + self.y0 = y0 # Start year + self.y1 = y1 # End year + self.k0 = k0 # Start depth (top) + self.k1 = k1 # End depth (bottom) + + # Specify region-specific coordinate boundaries + self.get_coordinates() + + # Compute volume per grid cell + self.get_volume() + + # Read temperature and salinity, averaged over specified time period + self.get_TS() + self.get_bins() + print(f'Got run {self.model} {self.region} {self.run} {self.y0}-{self.y1}') + + def get_coordinates(self): + # Hard-coded domain boundaries for specific regions. + # Currently in model-specific coordinates, + # Should be generalised for all models using a mask in lon-lat dimensions + if self.region == 'Amundsen': + if self.model == 'CESM2': + self.i0,self.i1,self.j0,self.j1 = 250,270,8,18 + elif self.model == 'UKESM1-0-LL': + self.i0,self.i1,self.j0,self.j1 = 170,190,50,70 + elif self.model == 'EC-Earth3': + self.i0,self.i1,self.j0,self.j1 = 170,190,12,32 + elif self.model == 'WOA23': + self.i0,self.i1,self.j0,self.j1 = -116,-96,-76,-70 + elif self.region == 'Ross': + if self.model == 'CESM2': + self.i0,self.i1,self.j0,self.j1 = 180,220,2,7 + elif self.model == 'UKESM1-0-LL': + self.i0,self.i1,self.j0,self.j1 = 85,140,38,48 + elif self.model == 'EC-Earth3': + self.i0,self.i1,self.j0,self.j1 = 85,140,0,10 + elif self.model == 'WOA23': + self.i0,self.i1,self.j0,self.j1 = -180,-150,-79,-76 + elif self.region == 'Weddell': + if self.model == 'CESM2': + self.i0,self.i1,self.j0,self.j1 = 300,360,3,9 + elif self.model == 'UKESM1-0-LL': + self.i0,self.i1,self.j0,self.j1 = 224,246,38,55 + elif self.model == 'EC-Earth3': + self.i0,self.i1,self.j0,self.j1 = 226,248,0,17 + elif self.model == 'WOA23': + self.i0,self.i1,self.j0,self.j1 = -65,-40,-78.5,-74.5 + elif self.region == 'Totten': + if self.model == 'CESM2': + self.i0,self.i1,self.j0,self.j1 = 137,146,24,28 + elif self.model == 'UKESM1-0-LL': + self.i0,self.i1,self.j0,self.j1 = 40,50,78,84 + elif self.model == 'EC-Earth3': + self.i0,self.i1,self.j0,self.j1 = 42,52,40,46 + elif self.model == 'WOA23': + self.i0,self.i1,self.j0,self.j1 = 113,124,-67.5,-64.5 + + def get_volume(self): + # Read volume and other relevant variables for volume-weighting of T and S + if self.model in ['CESM2']: + # Read all available volcello files + ds = xr.open_mfdataset(f'/usr/people/lambert/work2/data/cmip6/{self.model}/volcello*.nc') + + # Select specified domain (horizontal and vertical). Note: volume is time-independent for this model + ds = ds.sel(nlat=slice(self.j0,self.j1),nlon=slice(self.i0,self.i1),lev=slice(self.k0*100.,self.k1*100.)) + + # Get bunch of variables + self.lev = ds.lev_bnds.values + self.depth = -np.average(self.lev,axis=1) + self.lon = ds.lon.values + self.lat = ds.lat.values + self.V = ds.volcello.values + ds.close() + + # Read thetao files + ds = xr.open_mfdataset(f'/usr/people/lambert/work2/data/cmip6/{self.model}/historical/thetao/thetao*.nc') + ds = ds.sel(nlat=slice(self.j0,self.j1),nlon=slice(self.i0,self.i1),lev=slice(self.k0*100.,self.k1*100.)) + ds = ds.isel(time=0) + # Overwrite volume with zeros where thetao is nan + self.V = np.where(np.isnan(ds.thetao),0,self.V) + ds.close() + + elif self.model in ['EC-Earth3']: + + # Read thickness files + if self.y1-self.y0==1: + ds = xr.open_dataset(f'/usr/people/lambert/work2/data/cmip6/{self.model}/{self.run}/thkcello/thkcello_Omon_{self.model}_{self.run}_r1i1p1f1_gn_{self.y0}01-{self.y0}12.nc') + else: + ds = xr.open_mfdataset(f'/usr/people/lambert/work2/data/cmip6/{self.model}/{self.run}/thkcello/*.nc',combine='by_coords') + + # Select temporal and spatial domains + ds = ds.sel(time=slice(f'{self.y0}-01-01',f'{self.y1}-01-01'),j=slice(self.j0,self.j1),i=slice(self.i0,self.i1),lev=slice(self.k0,self.k1)) + if self.y1-self.y0==1: + self.lev = ds.lev_bnds.values + else: + self.lev = np.average(ds.lev_bnds.values,axis=0,weights=dpm(self.y1-self.y0)) + + # Get variables including thickness, time-averaged + self.depth = -np.average(self.lev,axis=1) + self.lon = ds.longitude.values + self.lat = ds.latitude.values + self.thkcello = np.average(ds.thkcello,axis=0,weights=dpm(self.y1-self.y0)) + ds.close() + + # Read horizontal area of cells + ds = xr.open_dataset(f'/usr/people/lambert/work2/data/cmip6/{self.model}/areacello_Ofx_EC-Earth3_historical_r1i1p1f1_gn.nc') + ds = ds.sel(j=slice(self.j0,self.j1),i=slice(self.i0,self.i1)) + + self.areacello = ds.areacello.values + ds.close() + + # Multiply thickness and area to get volume + self.V = self.thkcello*self.areacello + + # Set invalid cells to zero + self.V = np.where(np.isnan(self.V),0,self.V) + + elif self.model in ['UKESM1-0-LL']: + + # Read thickness files + ds = xr.open_mfdataset(f'/usr/people/lambert/work2/data/cmip6/{self.model}/{self.run}/thkcello/*.nc') + + # Select spatial and temporal domain + ds = ds.sel(time=slice(f'{self.y0}-01-01',f'{self.y1}-01-01'),j=slice(self.j0,self.j1),i=slice(self.i0,self.i1),lev=slice(self.k0,self.k1)) + + # Get bunch of variables, time-averaged + self.lev = np.average(ds.lev_bnds.values,axis=0,weights=dpm(self.y1-self.y0)) + self.depth = -np.average(self.lev,axis=1) + self.lon = ds.longitude.values + self.lat = ds.latitude.values + self.thkcello = np.average(ds.thkcello,axis=0,weights=dpm(self.y1-self.y0)) + ds.close() + + # Read horizontal area of cells + ds = xr.open_dataset(f'/usr/people/lambert/work2/data/cmip6/{self.model}/areacello_Ofx_UKESM1-0-LL_piControl_r1i1p1f2_gn.nc') + ds = ds.sel(j=slice(self.j0,self.j1),i=slice(self.i0,self.i1)) + + self.areacello = ds.areacello.values + ds.close() + + # Multiply thickness and area to get volume + self.V = self.thkcello*self.areacello + + # Set invalid cells to zero + self.V = np.where(np.isnan(self.V),0,self.V) + + elif self.model == 'WOA23': + # Read data file and select temporal and spatial domain + ds = xr.open_dataset('/usr/people/lambert/work/projects/data/woa23/woa23_decav91C0_t00_04.nc',decode_cf=False) + ds = ds.isel(time=0) + ds = ds.sel(lon=slice(self.i0,self.i1),lat=slice(self.j0,self.j1),depth=slice(self.k0,self.k1)) + + # Read temperature and set invalid values to nan + self.T = np.where(ds.t_an<1e30,ds.t_an,np.nan) + + self.lon = ds.lon.values + self.lat = ds.lat.values + self.depth = -ds.depth + + # Construct volume array + self.V = np.zeros((len(ds.depth),len(ds.lat),len(ds.lon))) + R = 6.371e6 + + # Define horizontal area + A = np.ones((len(ds.lat),len(ds.lon))) + for l in range(len(ds.lat)): + dy = 2*np.pi*R*(ds.lat_bnds[l,1]-ds.lat_bnds[l,0])/360 + dx = 2*np.pi*R*(ds.lon_bnds[:,1]-ds.lon_bnds[:,0])/360 * np.cos(2*np.pi*ds.lat[l]/360) + A[l,:] = dx*dy + + # Get thickness and multiply with area to get volume + for d in range(len(ds.depth)): + D = np.where(np.isnan(self.T[d,:,:]),0,ds.depth_bnds[d,1]-ds.depth_bnds[d,0]) + self.V[d,:,:] = D*A + ds.close() + + def get_TS(self): + # Read temperature and salinity from required model + # Fields are time-averaged over the required period + + if self.model in ['CESM2']: + ds = xr.open_mfdataset(f'/usr/people/lambert/work2/data/cmip6/{self.model}/{self.run}/thetao/*.nc',combine='by_coords') + ds = ds.sel(time=slice(f'{self.y0}-01-01',f'{self.y1}-01-01'),nlat=slice(self.j0,self.j1),nlon=slice(self.i0,self.i1),lev=slice(self.k0*100.,self.k1*100.)) + self.T = np.average(ds.thetao,axis=0,weights=dpm(self.y1-self.y0)) + ds.close() + + ds = xr.open_mfdataset(f'/usr/people/lambert/work2/data/cmip6/{self.model}/{self.run}/so/*.nc',combine='by_coords') + ds = ds.sel(time=slice(f'{self.y0}-01-01',f'{self.y1}-01-01'),nlat=slice(self.j0,self.j1),nlon=slice(self.i0,self.i1),lev=slice(self.k0*100.,self.k1*100.)) + self.S = np.average(ds.so,axis=0,weights=dpm(self.y1-self.y0)) + ds.close() + + elif self.model in ['EC-Earth3']: + if self.y1-self.y0==1: + ds = xr.open_dataset(f'/usr/people/lambert/work2/data/cmip6/{self.model}/{self.run}/thetao/thetao_Omon_{self.model}_{self.run}_r1i1p1f1_gn_{self.y0}01-{self.y0}12.nc') + else: + ds = xr.open_mfdataset(f'/usr/people/lambert/work2/data/cmip6/{self.model}/{self.run}/thetao/*.nc',combine='by_coords') + ds = ds.sel(time=slice(f'{self.y0}-01-01',f'{self.y1}-01-01'),j=slice(self.j0,self.j1),i=slice(self.i0,self.i1),lev=slice(self.k0,self.k1)) + self.T = np.average(ds.thetao,axis=0,weights=dpm(self.y1-self.y0)) + ds.close() + + if self.y1-self.y0==1: + ds = xr.open_dataset(f'/usr/people/lambert/work2/data/cmip6/{self.model}/{self.run}/so/so_Omon_{self.model}_{self.run}_r1i1p1f1_gn_{self.y0}01-{self.y0}12.nc') + else: + ds = xr.open_mfdataset(f'/usr/people/lambert/work2/data/cmip6/{self.model}/{self.run}/so/*.nc',combine='by_coords') + ds = ds.sel(time=slice(f'{self.y0}-01-01',f'{self.y1}-01-01'),j=slice(self.j0,self.j1),i=slice(self.i0,self.i1),lev=slice(self.k0,self.k1)) + self.S = np.average(ds.so,axis=0,weights=dpm(self.y1-self.y0)) + ds.close() + + elif self.model in ['UKESM1-0-LL']: + ds = xr.open_mfdataset(f'/usr/people/lambert/work2/data/cmip6/{self.model}/{self.run}/thetao/*.nc',combine='by_coords') + ds = ds.sel(time=slice(f'{self.y0}-01-01',f'{self.y1}-01-01'),j=slice(self.j0,self.j1),i=slice(self.i0,self.i1),lev=slice(self.k0,self.k1)) + self.T = np.average(ds.thetao,axis=0,weights=dpm(self.y1-self.y0)) + ds.close() + + ds = xr.open_mfdataset(f'/usr/people/lambert/work2/data/cmip6/{self.model}/{self.run}/so/*.nc',combine='by_coords') + ds = ds.sel(time=slice(f'{self.y0}-01-01',f'{self.y1}-01-01'),j=slice(self.j0,self.j1),i=slice(self.i0,self.i1),lev=slice(self.k0,self.k1)) + self.S = np.average(ds.so,axis=0,weights=dpm(self.y1-self.y0)) + ds.close() + + elif self.model == 'WOA23': + # Already got T, so only need to read S + ds = xr.open_dataset('/usr/people/lambert/work/projects/data/woa23/woa23_decav91C0_s00_04.nc',decode_cf=False) + ds = ds.isel(time=0) + ds = ds.sel(lon=slice(self.i0,self.i1),lat=slice(self.j0,self.j1),depth=slice(self.k0,self.k1)) + self.S = np.where(ds.s_an<1e30,ds.s_an,np.nan) + ds.close() + + def get_bins(self): + # Create 2D histogram of volume Vb in terms of binned salinity Sb and temperature Tb + self.Vb,self.Sb,self.Tb = np.histogram2d(self.S.flatten()[self.V.flatten()>0],self.T.flatten()[self.V.flatten()>0],bins=100,weights=self.V.flatten()[self.V.flatten()>0])