diff --git a/coast/__init__.py b/coast/__init__.py index 0e9ad3a9..f7050e2a 100644 --- a/coast/__init__.py +++ b/coast/__init__.py @@ -6,7 +6,8 @@ from .diagnostics.eof import compute_eofs, compute_hilbert_eofs from .diagnostics.internal_tide import InternalTide from .diagnostics.climatology import Climatology -from .diagnostics.annual_hydrographic_climatology import Annual_Climatology + +# from .diagnostics.annual_hydrographic_climatology import Annual_Climatology from ._utils import logging_util, general_utils, plot_util, crps_util from .data.index import Indexed from .data.profile import Profile @@ -27,7 +28,8 @@ from .data.copernicus import Copernicus, Product from ._utils.experiments_file_handling import experiments from ._utils.experiments_file_handling import nemo_filename_maker -from .diagnostics.hydrographic_profiles import Hydrographic_Profiles + +# from .diagnostics.hydrographic_profiles import Hydrographic_Profiles # Set default for logging level when coast is imported import logging diff --git a/coast/diagnostics/annual_hydrographic_climatology.py b/coast/diagnostics/annual_hydrographic_climatology.py deleted file mode 100644 index 36f7d69e..00000000 --- a/coast/diagnostics/annual_hydrographic_climatology.py +++ /dev/null @@ -1,76 +0,0 @@ -from ..data.gridded import Gridded -from ..diagnostics.internal_tide import InternalTide -import numpy as np -import xarray as xr - - -class Annual_Climatology(Gridded): - """ - Calculates a mean annual cycle from multi-annual monthly data - Because it calculates dervied properties (e.g PEA), data must be loaded. - Currently hardwired to calculate SST, SSS and PEA, placing these in the Gridded Objected - """ - - def __init__(self, gridded_t, gridded_t_out, Zmax=200.0): - """ - Args: - gridded_t: Input gridded object. - gridded_t: - Zmax: Max z. - """ - - # calculate a depth mask - Zd_mask, _, _ = gridded_t.calculate_vertical_mask(Zmax) - - ny = gridded_t.dataset.dims["y_dim"] - nx = gridded_t.dataset.dims["x_dim"] - - nt = gridded_t.dataset.dims["t_dim"] - - SSTy = np.zeros((12, ny, nx)) - SSSy = np.zeros((12, ny, nx)) - PEAy = np.zeros((12, ny, nx)) - # NBTy=np.zeros((12,ny,nx)) #will add near bed temperature later - - PEAy = np.zeros((12, ny, nx)) - - nyear = int(nt / 12) # hard wired for monthly data starting in Jan - for iy in range(nyear): - print("Calc PEA", iy) - it = np.arange((iy) * 12, (iy) * 12 + 12).astype(int) - for im in range(12): - itt = [it[im]] - print(itt) - gridded_t2 = gridded_t.subset_as_copy(t_dim=itt) - print("copied", im) - PEA = InternalTide(gridded_t2, gridded_t2) - PEA.calc_pea(gridded_t2, Zd_mask) - PEAy[im, :, :] = PEAy[im, :, :] + PEA.dataset["PEA"].values - PEAy = PEAy / nyear - - # need to find efficient method for bottom temperature - # NBT=np.zeros((nt,ny,nx)) - # for it in range(nt): - # NBT[it,:,:]=np.reshape(tmp[it,:,:,:].values.ravel()[Ikmax],(ny,nx)) - SST = gridded_t.dataset.variables["temperature"][:, 0, :, :] - SSS = gridded_t.dataset.variables["salinity"][:, 0, :, :] - - for im in range(12): - print("Month", im) - it = np.arange(im, nt, 12).astype(int) - SSTy[im, :, :] = np.mean(SST[it, :, :], axis=0) - SSSy[im, :, :] = np.mean(SSS[it, :, :], axis=0) - # NBTy[im,:,:]=np.mean(NBT[it,:,:],axis=0) - # save hard work in netcdf file - coords = { - "Months": (("mon_dim"), np.arange(12).astype(int)), - "latitude": (("y_dim", "x_dim"), gridded_t.dataset.latitude.values), - "longitude": (("y_dim", "x_dim"), gridded_t.dataset.longitude.values), - } - dims = ["mon_dim", "y_dim", "x_dim"] - attributes_SST = {"units": "o^C", "standard name": "Conservative Sea Surface Temperature"} - attributes_SSS = {"units": "", "standard name": "Absolute Sea Surface Salinity"} - attributes_PEA = {"units": "Jm^-3", "standard name": "Potential Energy Anomaly to " + str(Zmax) + "m"} - gridded_t_out.dataset["SSTy"] = xr.DataArray(np.squeeze(SSTy), coords=coords, dims=dims, attrs=attributes_SST) - gridded_t_out.dataset["SSSy"] = xr.DataArray(np.squeeze(SSSy), coords=coords, dims=dims, attrs=attributes_SSS) - gridded_t_out.dataset["PEAy"] = xr.DataArray(np.squeeze(PEAy), coords=coords, dims=dims, attrs=attributes_PEA) diff --git a/coast/diagnostics/hydrographic_profiles.py b/coast/diagnostics/hydrographic_profiles.py deleted file mode 100644 index 9f456be7..00000000 --- a/coast/diagnostics/hydrographic_profiles.py +++ /dev/null @@ -1,427 +0,0 @@ -import os -import numpy as np -import xarray as xr -import gsw -from typing import List - -from ..data.gridded import Gridded -from ..data.profile import Profile -from ..data.index import Indexed -from dask.diagnostics import ProgressBar - -# -Re = 6367456 * np.pi / 180 - - -class Hydrographic_Profiles(Indexed): - - ############################################################################### - def __init__(self, filename="none", datasetnames="none", config="", regionbounds=[]): - """Reads and manipulates lists of hydrographic profiles. - - Reads and manipulates lists of hydrographic profiles if called with datasetnames and regionbounds, - extract profiles in these bounds, and if a filenames is provided, saves them there. - """ - if datasetnames != "none" and len(regionbounds) == 4: - self.extractprofiles(datasetnames, regionbounds, config) - if filename != "none": - self.saveprofiles(filename) - - def extractprofiles(self, datasetnames, regionbounds, config): - """ - Args: - datasetnames: list of file names. - regionbounds: [lon min, lon max, lat min lat max] - config : a configuration file (optional) - """ - x_min = regionbounds[0] - x_max = regionbounds[1] - y_min = regionbounds[2] - y_max = regionbounds[3] - self.profile = Profile(config=config) - self.profile.read_en4(datasetnames, multiple=True) - self.profile = self.profile.subset_indices_lonlat_box(lonbounds=[x_min, x_max], latbounds=[y_min, y_max]) - self.profile = self.profile.process_en4() - - ######################################################################################## - def saveprofiles(self, filename): - """Saves profile and gridded objects to netcdf.""" - filename_profile = filename[:-3] + "_profile.nc" - filename_gridded = filename[:-3] + "_gridded.nc" - - print("saving Profile data") - with ProgressBar(): - self.profile.dataset.to_netcdf(filename_profile) - print("saving gridded data") - with ProgressBar(): - self.gridded.dataset.to_netcdf(filename_gridded) - - def loadprofiles(self, filename): - filename_profile = filename[:-3] + "_profile.nc" - filename_gridded = filename[:-3] + "_gridded.nc" - self.profile = Profile() - dataset = xr.load_dataset(filename_profile) - self.profile.insert_dataset(dataset) - dataset = xr.load_dataset(filename_gridded) - self.gridded.dataset = dataset - - ############################################################################## - def match_to_grid(self, gridded: Gridded, limits: List = [0, 0, 0, 0], rmax: int = 7000) -> None: - """Match profiles locations to grid, finding 4 nearest neighbours for each profile. - - Args: - gridded (Gridded): Gridded object. - limits (List): [jmin,jmax,imin,imax] - Subset to this region. - rmax (int): 7000 m maxmimum search distance. - """ - self.gridded = gridded - if sum(limits) != 0: - gridded.subset(ydim=range(limits[0], limits[1] + 0), xdim=range(limits[2], limits[3] + 1)) - # keep the grid or subset on the hydrographic profiles object - gridded.dataset["limits"] = limits - self.gridded = gridded - lon_prf = self.profile.dataset.longitude.values - lat_prf = self.profile.dataset.latitude.values - - # Find 4 nearest neighbours on grid - j_prf, i_prf, rmin_prf = gridded.find_j_i_list(lat=lat_prf, lon=lon_prf, n_nn=4) - - self.profile.dataset["i_min"] = limits[0] # reference back to origianl grid - self.profile.dataset["j_min"] = limits[2] - - i_min = self.profile.dataset.i_min.values - j_min = self.profile.dataset.j_min.values - - # Sort 4 NN by distance on grid - ii = np.nonzero(np.isnan(lon_prf)) - i_prf[ii, :] = 0 - j_prf[ii, :] = 0 - ip = np.where(np.logical_or(i_prf[:, 0] != 0, j_prf[:, 0] != 0))[0] - lon_prf4 = np.repeat(lon_prf[ip, np.newaxis], 4, axis=1).ravel() - lat_prf4 = np.repeat(lat_prf[ip, np.newaxis], 4, axis=1).ravel() - r = np.ones(i_prf.shape) * np.nan - lon_grd = gridded.dataset.longitude.values - lat_grd = gridded.dataset.latitude.values - - rr = Hydrographic_Profiles.distance_on_grid( - lat_grd, lon_grd, j_prf[ip, :].ravel(), i_prf[ip, :].ravel(), lat_prf4, lon_prf4 - ) - r[ip, :] = np.reshape(rr, (ip.size, 4)) - # sort by distance - ii = np.argsort(r, axis=1) - rmin_prf = np.take_along_axis(r, ii, axis=1) - i_prf = np.take_along_axis(i_prf, ii, axis=1) - j_prf = np.take_along_axis(j_prf, ii, axis=1) - - ii = np.nonzero(np.logical_or(np.min(r, axis=1) > rmax, np.isnan(lon_prf))) - i_prf = i_prf + i_min - j_prf = j_prf + j_min - i_prf[ii, :] = 0 # should the be nan? - j_prf[ii, :] = 0 - - self.profile.dataset["i_prf"] = xr.DataArray(i_prf, dims=["id_dim", "4"]) - self.profile.dataset["j_prf"] = xr.DataArray(j_prf, dims=["id_dim", "4"]) - self.profile.dataset["rmin_prf"] = xr.DataArray(rmin_prf, dims=["id_dim", "4"]) - - ############################################################################### - def stratificationmetrics(self, Zmax: int = 200, DZMAX: int = 30) -> None: - """Calculates various stratification metrics for observed profiles. - - Currently: PEA, PEAT, SST, SSS, NBT. - - Args: - Zmax = 200 m maximum depth of integration. - DZMAX = 30 m depth of surface layer. - """ - i_prf = self.profile.dataset.i_prf - self.profile.dataset.i_min - j_prf = self.profile.dataset.j_prf - self.profile.dataset.j_min - D = self.gridded.dataset.bathymetry # uses bathymetry from gridded object - i_prf = np.ma.masked_less(i_prf, 0) - j_prf = np.ma.masked_less(j_prf, 0) - - nprof = self.profile.dataset.dims["id_dim"] - nz = self.profile.dataset.dims["z_dim"] - sst = np.ones((nprof)) * np.nan - sss = np.ones((nprof)) * np.nan - nbt = np.ones((nprof)) * np.nan - kbot = np.ones((nprof), dtype=int) * np.nan - PEA = np.ones((nprof)) * np.nan - PEAT = np.ones((nprof)) * np.nan - quart = [0, 0.25, 0.5, 0.75, 1] - # fix memory issues for very large data sets, if this still needed with xarray? - if nprof < 1000000: - npr = nprof - else: - npr = int(nprof / 10) - - for ichnk in Hydrographic_Profiles.chunks(range(0, nprof), npr): - Ichnk = list(ichnk) - print(min(Ichnk), max(Ichnk)) - tmp = self.profile.dataset.potential_temperature[Ichnk, :].values - sal = self.profile.dataset.practical_salinity[Ichnk, :].values - ZZ = -self.profile.dataset.depth[Ichnk, :].values - lat = self.profile.dataset.latitude[Ichnk].values - lon = self.profile.dataset.longitude[Ichnk].values - rmin = self.profile.dataset.rmin_prf[Ichnk, :].values - nprof = len(Ichnk) - Zd_mask = np.zeros((nprof, nz)) - - ################################################################################ - # define interface layers and DZ associated with Z - Zw = np.empty_like(ZZ) - DZ = np.empty_like(ZZ) - Zw[:, 0] = 0.0 - I = np.arange(0, nz - 1) - Zw[:, I + 1] = 0.5 * (ZZ[:, I] + ZZ[:, I + 1]) - DZ[:, I] = Zw[:, I] - Zw[:, I + 1] - DZ[~np.isfinite(DZ)] = 0.0 - ZZ[~np.isfinite(ZZ)] = 0.0 - DP = np.ones((nprof)) * np.nan - # depth from model - print("Depth from model") - for ip in range(nprof): - - DP[ip] = 0.0 - rr = 0.0 - for iS in range(0, 4): - if D[j_prf[ip, iS], i_prf[ip, iS]] != 0: - DP[ip] = DP[ip] + D[j_prf[ip, iS], i_prf[ip, iS]] / rmin[ip, iS] - rr = rr + 1 / rmin[ip, iS] - if rr != 0.0: - DP[ip] = DP[ip] / rr - print("define good profiles") - good_profile = np.zeros((nprof)) - sstc = np.ones((nprof)) * np.nan - sssc = np.ones((nprof)) * np.nan - nbtc = np.ones((nprof)) * np.nan - kbot = np.ones((nprof), dtype=int) * np.nan - T = np.zeros(nz) * np.nan - S = np.zeros(nz) * np.nan - Z = np.zeros(nz) * np.nan - ZW = np.zeros(nz) * np.nan - DP[DP == 0] = np.nan - - for ip in range(nprof): - - Dp = DP[ip] - T[:] = tmp[ip, :] - S[:] = sal[ip, :] - # Z always -ve downwards - Z[:] = -np.abs(ZZ[ip, :]) - ZW[:] = -np.abs(Zw[ip, :]) - I = np.nonzero(np.isfinite(T))[0] - - if np.size(I) > 0 and np.isfinite(Dp): - kbot[ip] = np.max(I) - - # SST - if -Z[np.min(I)] < np.min([DZMAX, 0.25 * Dp]): - sstc[ip] = T[np.min(I)] - # SSS - if -Z[np.min(I)] < np.min([DZMAX, 0.25 * Dp]): - sssc[ip] = S[np.min(I)] - # Near bototm or ~Zmax temp. - if Dp < Zmax: - if Dp + Z[int(kbot[ip])] < np.min([DZMAX, 0.25 * Dp]): - nbtc[ip] = T[np.max(I)] - elif kbot[ip] == nz - 1: - nbtc[ip] = T[int(kbot[ip])] - elif Z[int(kbot[ip])] < -Zmax and np.size(np.nonzero(Z[I] > -Zmax)) != 0: - k = np.max(np.nonzero(Z[I] > -Zmax)[0]) - k = int(I[k]) - r = (-Zmax - Z[k]) / (Z[k + 1] - Z[k]) - nbt[ip] = T[k] * r + T[k + 1] * (1.0 - r) - - # Depth mask - Zd_mask[ip, 0 : int(kbot[ip])] = 1 - Imax = np.max(np.nonzero(ZW > -Zmax)[0]) # note ZW index - if Imax < kbot[ip]: - Zd_mask[ip, Imax:nz] = 0 - Zd_mask[ip, Imax] = (ZW[Imax] - (-Zmax)) / (ZW[Imax] - ZW[Imax + 1]) - if Zd_mask[ip, Imax - 1] < 0 or Zd_mask[ip, Imax - 1] > 1: - print("error", ip, Zd_mask[ip, Imax - 1]), Imax, kbot[ip] - # find good profiles - - DD = np.min([Dp, Zmax]) - good_profile[ip] = 1 - for iq in range(len(quart) - 1): - I = np.nonzero( - np.all(np.concatenate(([Z <= -DD * quart[iq]], [Z >= -DD * quart[iq + 1]]), axis=0), axis=0) - ) - - if np.size(I) == 0: - good_profile[ip] = 0 - elif ~(np.any((np.isfinite(S[I]))) and np.any((np.isfinite(S[I])))): - good_profile[ip] = 0 - ### - - T = Hydrographic_Profiles.fillholes(T) - S = Hydrographic_Profiles.fillholes(S) - tmp[ip, :] = T - sal[ip, :] = S - - ############################################################################### - print("Calculate metrics") - metrics = Hydrographic_Profiles.profile_metrics(tmp, sal, ZZ, DZ, Zd_mask, lon, lat) - - PEAc = metrics["PEA"] - PEATc = metrics["PEAT"] - PEAc[good_profile == 0] = np.nan - PEATc[good_profile == 0] = np.nan - sst[Ichnk] = sstc - sss[Ichnk] = sssc - nbt[Ichnk] = nbtc - PEA[Ichnk] = PEAc - PEAT[Ichnk] = PEATc - # Next chunk - - DT = sst - nbt - self.profile.dataset["PEA"] = xr.DataArray(PEA, dims=["id_dim"]) - self.profile.dataset["PEAT"] = xr.DataArray(PEAT, dims=["id_dim"]) - self.profile.dataset["SST"] = xr.DataArray(sst, dims=["id_dim"]) - self.profile.dataset["SSS"] = xr.DataArray(sss, dims=["id_dim"]) - self.profile.dataset["NBT"] = xr.DataArray(nbt, dims=["id_dim"]) - self.profile.dataset["DT"] = xr.DataArray(DT, dims=["id_dim"]) - - def grid_hydro_mnth(self): - i_prf = self.profile.dataset.i_prf.values[:, 0] - j_prf = self.profile.dataset.j_prf.values[:, 0] - varnames = ["SST", "SSS", "PEA", "PEAT", "DT", "NBT"] - for varname in varnames: - print("Gridding", varname) - mnth = self.profile.dataset.time.values.astype("datetime64[M]").astype(int) % 12 + 1 - var, nvar = Hydrographic_Profiles.grid_vars_mnth(self, varname, i_prf, j_prf, mnth) - self.gridded.dataset[varname] = xr.DataArray(var, dims=["12", "y_dim", "x_dim"]) - self.gridded.dataset["n" + varname] = xr.DataArray(nvar, dims=["12", "y_dim", "x_dim"]) - - ############################################################################### - @staticmethod - def makefilenames(path, dataset, yr_start, yr_stop): - if dataset == "EN4": - datasetnames = [] - january = 1 - december = 13 # range is non-inclusive so we need 12 + 1 - for yr in range(yr_start, yr_stop + 1): - for im in range(january, december): - name = os.path.join(path, f"EN.4.2.1.f.profiles.l09.{yr}{im:02}.nc") - datasetnames.append(name) - return datasetnames - print("Data set not coded") - - # Functions - ############################################################################### - # Functions for match to grid - @staticmethod - def subsetgrid(var_dom, limits): - i_min = limits[0] - i_max = limits[1] - j_min = limits[2] - j_max = limits[3] - if i_max > i_min: - return var_dom[i_min : i_max + 1, j_min : j_max + 1] - # special case for wrap-around - gvar1 = var_dom[i_min:, j_min : j_max + 1] - gvar2 = var_dom[:i_max, j_min : j_max + 1] - var_dom = np.concatenate((gvar1, gvar2), axis=0) - return var_dom - - ############################################################################### - ########################################### - def distance_on_grid(Y, X, jpts, ipts, Ypts, Xpts): - DX = (Xpts - X[jpts, ipts]) * Re * np.cos(Ypts * np.pi / 180.0) - DY = (Ypts - Y[jpts, ipts]) * Re - r = np.sqrt(DX**2 + DY**2) - return r - - ############################################################################### - # Functions for stratification metrics - @staticmethod - def fillholes(Y): - YY = np.ones(np.shape(Y)) - YY[:] = Y - I = np.nonzero(np.isfinite(YY)) - N = len(YY) - - if np.size(I) > 0: - if not np.isfinite(YY[0]): - YY[0 : np.min(I) + 1] = YY[np.min(I)] - - if ~np.isfinite(YY[N - 1]): - YY[np.max(I) : N] = YY[np.max(I)] - I = np.array(np.nonzero(~np.isfinite(YY))) - YY[I] = 0.5 * (YY[I - 1] + YY[I + 1]) - YYp = YY[0] - ip = 0 - for i in range(N): - if np.isfinite(YY[i]): - YYp = YY[i] - ip = i - else: - j = i - while ~np.isfinite(YY[j]): - j = j + 1 - Jp = np.arange(ip + 1, j - 1 + 1) - - pT = np.arange(1.0, (j - ip - 1.0) + 1.0) / (j - ip) - YY[Jp] = YYp + (YY[j] - YYp) * pT - return YY - - ########################################### - def chunks(lst, n): - """Yield successive n-sized chunks from lst.""" - for i in range(0, len(lst), n): - yield lst[i : i + n] - - ########################################### - @staticmethod - def profile_metrics(tmp, sal, Z, DZ, Zd_mask, lon, lat): - - metrics = {} - g = 9.81 - DD = np.sum(DZ * Zd_mask, axis=1) - nz = Z.shape[1] - lat = np.repeat(lat[:, np.newaxis], nz, axis=1) - lon = np.repeat(lon[:, np.newaxis], nz, axis=1) - pressure_absolute = gsw.p_from_z(Z, lat) - salinity_absolute = gsw.SA_from_SP(sal, pressure_absolute, lon, lat) - temp_conservative = gsw.CT_from_pt(salinity_absolute, tmp) - rho = np.ma.masked_invalid(gsw.rho(salinity_absolute, temp_conservative, 0.0)) - - Tbar = np.sum(temp_conservative * DZ * Zd_mask, axis=1) / DD - Sbar = np.sum(salinity_absolute * DZ * Zd_mask, axis=1) / DD - - rhobar = np.ma.masked_invalid(gsw.rho(Sbar, Tbar, 0.0)) - rhobar_2d = np.repeat(rhobar[:, np.newaxis], nz, axis=1) - Sbar_2d = np.repeat(Sbar[:, np.newaxis], nz, axis=1) - rhoT = np.ma.masked_invalid(gsw.rho(Sbar_2d, temp_conservative, 0.0)) # density with constant salinity - - PEA = -np.sum(Z * (rho - rhobar_2d) * DZ * Zd_mask, axis=1) * g / DD - PEAT = -np.sum(Z * (rhoT - rhobar_2d) * DZ * Zd_mask, axis=1) * g / DD - - metrics["PEA"] = PEA - metrics["PEAT"] = PEAT - - return metrics - - ########################################### - - def grid_vars_mnth(self, var, i_var, j_var, mnth_var): - VAR = self.profile.dataset[var].values - nx = self.gridded.dataset.dims["x_dim"] - ny = self.gridded.dataset.dims["y_dim"] - - Ig = np.nonzero(np.isfinite(VAR))[0] - - var = VAR[Ig] - VAR_g = np.zeros((12, ny, nx)) - nVAR_g = np.zeros((12, ny, nx)) - for ip in range(0, np.size(Ig)): - i = i_var[Ig[ip]] - j = j_var[Ig[ip]] - im = int(mnth_var[Ig[ip]]) - 1 - - VAR_g[im, j, i] = VAR_g[im, j, i] + var[ip] - nVAR_g[im, j, i] = nVAR_g[im, j, i] + 1 - - VAR_g = VAR_g / nVAR_g - return VAR_g, nVAR_g