diff --git a/README.rst b/README.rst index 39e5c72c..afde56f9 100644 --- a/README.rst +++ b/README.rst @@ -3,6 +3,8 @@ Community Analysis Pipeline (CAP) Welcome to the Mars Climate Modeling Center (MCMC) **Community Analysis Pipeline (CAP)**. +For instructions and documentation please see our `online documentation `_. + About ----- **CAP** is a set of Python3 libraries and command-line executables that streamline downloading, processing, and plotting output from the NASA Ames Mars Global Climate Models: diff --git a/amescap/FV3_utils.py b/amescap/FV3_utils.py index b53396a8..8acf91fd 100644 --- a/amescap/FV3_utils.py +++ b/amescap/FV3_utils.py @@ -1275,7 +1275,13 @@ def mass_stream(v_avg, lat, level, type="pstd", psfc=700, H=8000., .g. ⌡ f(z)dz = (Zn-Zn-1){f(Zn) + f(Zn-1)}/2 n-1 """ - + reverting = False + if level[0] < level[-1]: + reverting = True + print("Reversing pstd array for mass stream function calculation...") + level = level[::-1] + v_avg = v_avg[::-1, ...] + g = 3.72 # m/s2 a = 3400*1000 # m nlev = len(level) @@ -1335,6 +1341,11 @@ def mass_stream(v_avg, lat, level, type="pstd", psfc=700, H=8000., MSF[mask] = np.nan if isMasked: MSF = np.ma.array(MSF, mask = mask) + + if reverting: + print("Reversing pstd dimension of MSF array for compatibility...") + MSF = MSF[::-1, ...] + return MSF.reshape(shape_out) @@ -1916,6 +1927,12 @@ def daily_to_diurn(varIN, time_in): iperday = int(np.round(1/dt_in)) vshape_in = varIN.shape + + # Add safety check for integer sols + if not(np.mod(vshape_in[0],iperday) == 0): + print("Error: File cannot be split evenly into sols") + return None + vreshape = np.append([-1, iperday], vshape_in[1:]).astype(int) varOUT = varIN.reshape(vreshape) diff --git a/amescap/Ncdf_wrapper.py b/amescap/Ncdf_wrapper.py index 3f4c7e88..e78a52df 100644 --- a/amescap/Ncdf_wrapper.py +++ b/amescap/Ncdf_wrapper.py @@ -119,9 +119,9 @@ def add_constant(self, variable_name, value, longname_txt="", # Private Definitions def _def_variable(self, variable_name, dim_array, longname_txt="", - units_txt=""): + units_txt="",datatype="f4"): self.var_dict[variable_name] = self.f_Ncdf.createVariable(variable_name, - "f4", + datatype, dim_array) self.var_dict[variable_name].units = units_txt self.var_dict[variable_name].long_name = longname_txt @@ -129,9 +129,9 @@ def _def_variable(self, variable_name, dim_array, longname_txt="", def _def_axis1D(self, variable_name, dim_array, longname_txt="", - units_txt="", cart_txt=""): + units_txt="", cart_txt="",datatype="f8"): self.var_dict[variable_name] = self.f_Ncdf.createVariable(variable_name, - "f4", + datatype, dim_array) self.var_dict[variable_name].units = units_txt self.var_dict[variable_name].long_name = longname_txt @@ -167,17 +167,16 @@ def _is_cart_axis(self, Ncvar): def log_variable(self, variable_name, DATAin, dim_array, longname_txt="", - units_txt=""): + units_txt="",datatype="f4"): """ EX:: Log.log_variable("sfcT", sfcT, ("time", "Nx"), "soil temperature", "K") """ - if variable_name not in self.var_dict.keys(): self._def_variable(variable_name, dim_array, longname_txt, - units_txt) + units_txt,datatype) self.var_dict[variable_name].long_name = longname_txt self.var_dict[variable_name].dim_name = str(dim_array) self.var_dict[variable_name].units = units_txt @@ -263,7 +262,7 @@ def copy_Ncvar(self, Ncvar, swap_array=None): longname_txt = getattr(Ncvar, "long_name", Ncvar._name) units_txt = getattr(Ncvar, "units", "") self._def_variable(Ncvar._name, Ncvar.dimensions, longname_txt, - units_txt) + units_txt,Ncvar.dtype) if np.any(swap_array): self.log_variable(Ncvar._name, swap_array[:], Ncvar.dimensions, longname_txt, units_txt) diff --git a/amescap/Script_utils.py b/amescap/Script_utils.py index 166b0ebf..5bb0e483 100644 --- a/amescap/Script_utils.py +++ b/amescap/Script_utils.py @@ -928,7 +928,6 @@ def filter_vars(fNcdf, include_list=None, giveExclude=False): else: print(f"{Yellow}***Warning*** from filter_vars(): " f"{ivar} not found in file.\n") - exit() baseline_var = [] for ivar in var_list: diff --git a/amescap/Spectral_utils.py b/amescap/Spectral_utils.py index ad270cea..66c4b5c6 100644 --- a/amescap/Spectral_utils.py +++ b/amescap/Spectral_utils.py @@ -31,27 +31,40 @@ # DEFINITIONS # ====================================================================== -# Try to import pyshtools with proper error handling -try: - import pyshtools - PYSHTOOLS_AVAILABLE = True -except ImportError: - PYSHTOOLS_AVAILABLE = False - print( - f"{Yellow}__________________\n" - f"Zonal decomposition relies on the pyshtools library, " - f"referenced at:\n\n" - f"Mark A. Wieczorek and Matthias Meschede (2018). " - f"SHTools - Tools for working with spherical harmonics," - f"Geochemistry, Geophysics, Geosystems, 2574-2592, " - f"doi:10.1029/2018GC007529\n\nPlease consult pyshtools " - f"documentation at:\n" - f" {Cyan}https://pypi.org/project/pyshtools\n" - f"{Yellow}And installation instructions for CAP with pyshtools:\n" - f" {Cyan}https://amescap.readthedocs.io/en/latest/installation." - f"html#_spectral_analysis{Yellow}\n" - f"__________________{Nclr}\n\n" - ) +def import_pyshtools(): + """ + Attempt to import pyshtools and set the global variable + PYSHTOOLS_AVAILABLE accordingly. + Check if pyshtools is available and return its availability status + + :return: True if pyshtools is available, False otherwise + :rtype: bool + """ + + global PYSHTOOLS_AVAILABLE + # Try to import pyshtools with proper error handling + try: + import pyshtools + global pyshtools + PYSHTOOLS_AVAILABLE = True + except ImportError: + PYSHTOOLS_AVAILABLE = False + print( + f"{Yellow}__________________\n" + f"Tidal decomposition relies on the pyshtools library, " + f"referenced at:\n" + f"{Cyan}Mark A. Wieczorek and Matthias Meschede (2018). " + f"SHTools - Tools for working with spherical harmonics, " + f"Geochemistry, Geophysics, Geosystems, 2574-2592, " + f"doi:10.1029/2018GC007529\n\n" + f"Please consult pyshtools documentation at:\n" + f"{Cyan}https://pypi.org/project/pyshtools\n\n" + f"{Yellow}And installation instructions for CAP with pyshtools:\n" + f"{Cyan}https://amescap.readthedocs.io/en/latest/installation." + f"html#_spectral_analysis{Yellow}\n" + f"__________________{Nclr}\n\n" + ) + return PYSHTOOLS_AVAILABLE def diurn_extract(VAR, N, tod, lon): @@ -74,7 +87,8 @@ def diurn_extract(VAR, N, tod, lon): (e.g., size ``[Nh, time, lat, lon]``) :rtype: ND arrays """ - + import_pyshtools() + dimsIN = VAR.shape nsteps = len(tod) period = 24 @@ -140,6 +154,7 @@ def diurn_extract(VAR, N, tod, lon): # Return the phase and amplitude return amp.reshape(dimsOUT), phas.reshape(dimsOUT) + def reconstruct_diurn(amp, phas, tod, lon, sumList=[]): """ Reconstructs a field wave based on its diurnal harmonics @@ -162,7 +177,8 @@ def reconstruct_diurn(amp, phas, tod, lon, sumList=[]): aggregated (i.e., size = ``[tod, time, lat, lon]``) :rtype: _type_ """ - + import_pyshtools() + dimsIN = amp.shape N = dimsIN[0] dimsSUM = np.append([len(tod)], dimsIN[1:]) @@ -207,21 +223,23 @@ def reconstruct_diurn(amp, phas, tod, lon, sumList=[]): # Return harmonics separately return varOUT -def space_time(lon, timex, varIN, kmx, tmx): + +def extract_diurnal_harmonics(kmx, tmx, varIN, tod, lon): """ Obtain west and east propagating waves. This is a Python implementation of John Wilson's ``space_time`` routine. - Alex Kling 2019. + Richard Urata 2025. :param lon: longitude [°] (0-360) :type lon: 1D array - :param timex: time [sol] (e.g., 1.5 days sampled every hour is - ``[0/24, 1/24, 2/24,.. 1,.. 1.5]``) - :type timex: 1D array - :param varIN: variable for the Fourier analysis. First axis must be - ``lon`` and last axis must be ``time`` (e.g., - ``varIN[lon, time]``, ``varIN[lon, lat, time]``, or - ``varIN[lon, lev, lat, time]``) + :param tod: time_of_day [hour] (e.g., a 24 hour sol sampled every hour could be + ``[0.5, 1.5, 2.5, ..., 23.5]``) + :type tod: 1D array + :param varIN: variable for the Fourier analysis, expected to be scaled by zonal mean. + First axis must be + ``time`` and last axis must be ``lon`` (e.g., + ``varIN[time, tod, lon]``, ``varIN[time, tod, lat, lon]``, or + ``varIN[time, tod, lat, lon, lev]``) :type varIN: array :param kmx: the number of longitudinal wavenumbers to extract (max = ``nlon/2``) @@ -236,7 +254,7 @@ def space_time(lon, timex, varIN, kmx, tmx): .. note:: 1. ``ampe``, ``ampw``, ``phasee``, and ``phasew`` have dimensions ``[kmx, tmx]`` or ``[kmx, tmx, lat]`` or - ``[kmx, tmx, lev, lat]`` etc.\n + ``[kmx, tmx, lat, lev]`` etc.\n 2. The x and y axes may be constructed as follows, which will display the eastern and western modes:: @@ -247,125 +265,71 @@ def space_time(lon, timex, varIN, kmx, tmx): phase = np.concatenate((phasew[:, ::-1], phasee), axis = 1) """ - # Get input variable dimensions - dims = varIN.shape - lon_id = dims[0] - time_id = dims[-1] - - # Additional dimensions stacked in the middle - dim_sup_id = dims[1:-1] - - # jd = total number of dimensions in the middle (``varIN > 3D``) - jd = int(np.prod(dim_sup_id)) - - # Flatten the middle dimensions, if any - varIN = np.reshape(varIN, (lon_id, jd, time_id)) - - # Initialize 4 empty arrays - ampw, ampe, phasew, phasee = ( - [np.zeros((kmx, tmx, jd)) for _x in range(0, 4)] - ) - - #TODO not implemented yet: - # zamp, zphas = [np.zeros((jd, tmx)) for _x in range(0, 2)] - - tpi = 2*np.pi - # Normalize longitude array - argx = lon * 2*np.pi / 360 - rnorm = 2. / len(argx) - # If timex = [0/24, 1/24, 2/24,.. 1] arg cycles for m [0, 2Pi] - arg = timex * 2*np.pi - # Nyquist cut off - rnormt = 2. / len(arg) - - for kk in range(0, kmx): - progress(kk, kmx) - cosx = np.cos(kk * argx) * rnorm - sinx = np.sin(kk * argx) * rnorm - - # Inner product calculates the Fourier coefficients of the - # cosine and sine contributions of the spatial variation - acoef = np.dot(varIN.T, cosx) - bcoef = np.dot(varIN.T, sinx) - - for nn in range(0, tmx): - # Get the cosine and sine series expansions of the temporal - # variations of the acoef and bcoef spatial terms - cosray = rnormt * np.cos(nn * arg) - sinray = rnormt * np.sin(nn * arg) - - cosA = np.dot(acoef.T, cosray) - sinA = np.dot(acoef.T, sinray) - cosB = np.dot(bcoef.T, cosray) - sinB = np.dot(bcoef.T, sinray) - - wr = 0.5*(cosA - sinB) - wi = 0.5*(-sinA - cosB) - er = 0.5*(cosA + sinB) - ei = 0.5*(sinA - cosB) - - aw = np.sqrt(wr**2 + wi**2) - ae = np.sqrt(er**2 + ei**2) - pe = np.arctan2(ei, er) * 180/np.pi - pw = np.arctan2(wi, wr) * 180/np.pi - - pe = np.mod(-np.arctan2(ei, er) + tpi, tpi) * 180/np.pi - pw = np.mod(-np.arctan2(wi, wr) + tpi, tpi) * 180/np.pi - - ampw[kk, nn, :] = aw.T - ampe[kk, nn, :] = ae.T - phasew[kk, nn, :] = pw.T - phasee[kk, nn, :] = pe.T - - ampw = np.reshape(ampw, (kmx, tmx) + dim_sup_id) - ampe = np.reshape(ampe, (kmx, tmx) + dim_sup_id) - phasew = np.reshape(phasew, (kmx, tmx) + dim_sup_id) - phasee = np.reshape(phasee, (kmx, tmx) + dim_sup_id) - - # TODO implement zonal mean: zamp, zphas (standing wave k = 0, - # zonally averaged) stamp, stphs (stationary component ktime = 0) - - # # varIN = reshape(varIN, dims) - # # if nargout < 5: - # # # only ampe, ampw, phasee, phasew are requested - # # return - - # # Now calculate the axisymmetric tides zamp,zphas - - # zvarIN = np.mean(varIN, axis=0) - # zvarIN = np.reshape(zvarIN, (jd, time_id)) - - # arg = timex * 2*np.pi - # arg = np.reshape(arg, (len(arg), 1)) - # rnorm = 2/time_id - - # for nn in range(0, tmx): - # cosray = rnorm * np.cos(nn*arg) - # sinray = rnorm * np.sin(nn*arg) - # cosser = np.dot(zvarIN, cosray) - # sinser = np.dot(zvarIN, sinray) - - # zamp[:, nn] = np.sqrt(cosser[:]**2 + sinser[:]**2).T - # zphas[:, nn] = np.mod(-np.arctan2(sinser, cosser) - # + tpi, tpi).T * 180/np.pi - - # zamp = zamp.T - # zphas = zphas.T - - # if len(dims) > 2: - # zamp = np.reshape(zamp, (tmx,) + dim_sup_id) - # zamp = np.reshape(zphas, (tmx,) + dim_sup_id) - - # # if nargout < 7: - # # return - - # sxx = np.mean(varIN, ndims(varIN)) - # [stamp, stphs] = amp_phase(sxx, lon, kmx) - - # if len(dims) > 2: - # stamp = reshape(stamp, [kmx dims(2:end-1)]) - # stphs = reshape(stphs, [kmx dims(2:end-1)]) - return ampe, ampw, phasee, phasew + # Ensure time is in days + if np.max(tod) > 1: + tod = tod / 24 + + # Reshape varIN to (...,tod) + if varIN.ndim == 3: + varIN = np.transpose(varIN, (0, 2, 1)) + ntime, nlon, ntod = varIN.shape + output_shape = (ntime, kmx, tmx) + if varIN.ndim == 4: + varIN = np.transpose(varIN, (0, 3, 2, 1)) + ntime, nlon, nlat, ntod = varIN.shape + output_shape = (ntime, kmx, tmx, nlat) + if varIN.ndim == 5: + varIN = np.transpose(varIN, (0, 4, 3, 2, 1)) + ntime, nlon, nlat, nlev, ntod = varIN.shape + output_shape = (ntime, kmx, tmx, nlat, nlev) + + # Convert longitude to radians + if np.max(lon) > 100: + lon_rad = np.deg2rad(lon) + else: + lon_rad = lon + + # Prepare arrays for results + ampe = np.zeros(output_shape) + ampw = np.zeros(output_shape) + phasee = np.zeros(output_shape) + phasew = np.zeros(output_shape) + + # Constants + tpi = 2 * np.pi + rnorm = 2 / nlon + rnormt = 2 / ntod + + # Main calculation loop + for t in range(ntime): + for k in range(kmx): # Changed to start from 0 + cosx = np.cos((k + 1) * lon_rad) * rnorm # Add 1 to k for the calculation + sinx = np.sin((k + 1) * lon_rad) * rnorm + + acoef = np.tensordot(varIN[t,...], cosx, axes=([0], [0])) + bcoef = np.tensordot(varIN[t,...], sinx, axes=([0], [0])) + + for n in range(tmx): # Changed to start from 0 + arg = (n + 1) * tpi * tod # Add 1 to n for the calculation + cosray = rnormt * np.cos(arg) + sinray = rnormt * np.sin(arg) + + cosA = np.dot(acoef, cosray) + sinA = np.dot(acoef, sinray) + cosB = np.dot(bcoef, cosray) + sinB = np.dot(bcoef, sinray) + + wr = 0.5 * (cosA - sinB) + wi = 0.5 * (-sinA - cosB) + er = 0.5 * (cosA + sinB) + ei = 0.5 * (sinA - cosB) + + ampw[t, k, n, ...] = np.sqrt(wr**2 + wi**2) + ampe[t, k, n, ...] = np.sqrt(er**2 + ei**2) + phasew[t, k, n, ...] = np.mod(-np.arctan2(wi, wr) + tpi, tpi) * 180 / np.pi + phasee[t, k, n, ...] = np.mod(-np.arctan2(ei, er) + tpi, tpi) * 180 / np.pi + + return np.squeeze(ampe), np.squeeze(ampw), np.squeeze(phasee), np.squeeze(phasew) def zeroPhi_filter(VAR, btype, low_highcut, fs, axis=0, order=4, @@ -397,6 +361,7 @@ def zeroPhi_filter(VAR, btype, low_highcut, fs, axis=0, order=4, .. note:: ``Wn=[low, high]`` are expressed as a function of the Nyquist frequency """ + import_pyshtools() # Create the filter low_highcut = np.array(low_highcut) @@ -434,6 +399,7 @@ def zonal_decomposition(VAR): .. note:: Output size is (``[...,lat/2, lat/2]``) as lat is the smallest dimension. This matches the Nyquist frequency. """ + import_pyshtools() if not PYSHTOOLS_AVAILABLE: raise ImportError( @@ -493,6 +459,7 @@ def zonal_construct(COEFFS_flat, VAR_shape, btype=None, low_highcut=None): print("(kmin,kmax) = ({kmin}, {kmax}) -> dx min = {L_min} km, dx max = {L_max} km") """ + import_pyshtools() if not PYSHTOOLS_AVAILABLE: raise ImportError( @@ -525,12 +492,3 @@ def zonal_construct(COEFFS_flat, VAR_shape, btype=None, low_highcut=None): return VAR.reshape(VAR_shape) -def init_shtools(): - """ - Check if pyshtools is available and return its availability status - - :return: True if pyshtools is available, False otherwise - :rtype: bool - """ - - return PYSHTOOLS_AVAILABLE diff --git a/bin/MarsFiles.py b/bin/MarsFiles.py index 93d38713..3aee1fe8 100755 --- a/bin/MarsFiles.py +++ b/bin/MarsFiles.py @@ -27,6 +27,7 @@ * ``[-tide, --tide_decomp]`` Extract diurnal tide and its harmonics * ``[-recon, --reconstruct]`` Reconstruct the first N harmonics * ``[-norm, --normalize]`` Provide ``-tide`` result in % amplitude + * ``[-prop, --prop_tides]`` Extract propagating tide harmonics * ``[-regrid, --regrid_XY_to_match]`` Regrid a target file to match a source file * ``[-zavg, --zonal_average]`` Zonally average all variables in a file * ``[-incl, --include]`` Only include specific variables in a calculation @@ -62,6 +63,7 @@ import shutil # For OS-friendly file operations import functools # For function decorators import traceback # For printing stack traces +import shutil # For copy/pasting fixed file after -split # Load amesCAP modules from amescap.Ncdf_wrapper import (Ncdf, Fort) @@ -283,7 +285,10 @@ def wrapper(*args, **kwargs): ) ) -parser.add_argument('-c', '--concatenate', action='store_true', +parser.add_argument('-c', '--concatenate', action=ExtAction, + ext_content='_concatenated', + parser=parser, + nargs=0, help=( f"Combine sequential files of the same type into one file.\n" f"Works on 'daily', 'diurn', and 'average' files.\n" @@ -291,13 +296,6 @@ def wrapper(*args, **kwargs): f"> ls\n" f"00334.atmos_average.nc 00668.atmos_average.nc\n" f"> MarsFiles *.atmos_average.nc -c\n" - f"{Blue}Overwrites 00334.atmos_average.nc with concatenated " - f"files:{Green}\n" - f"> ls\n" - f" 00334.atmos_average.nc\n" - f"{Yellow}To preserve original files, use [-ext --extension]:" - f"{Green}\n" - f"> MarsFiles *.atmos_average.nc -c -ext _concatenated\n" f"{Blue}Produces 00334.atmos_average_concatenated.nc and " f"preserves all other files:{Green}\n" f"> ls\n" @@ -446,7 +444,7 @@ def wrapper(*args, **kwargs): f"in Sols.\n" f"{Yellow}Generates a new file ending in ``_hps.nc``\n" f"{Green}Example:\n" - f"> MarsFiles 01336.atmos_daily.nc -hps 10 -add_trend\n" + f"> MarsFiles 01336.atmos_daily.nc -hps 10\n" f"{Nclr}\n\n" ) ) @@ -466,7 +464,7 @@ def wrapper(*args, **kwargs): f"cutoff frequency in Sols.\n" f"{Yellow}Generates a new file ending in ``_lps.nc``\n" f"{Green}Example:\n" - f"> MarsFiles 01336.atmos_daily.nc -lps 20 -add_trend\n" + f"> MarsFiles 01336.atmos_daily.nc -lps 20\n" f"{Nclr}\n\n" ) ) @@ -486,7 +484,7 @@ def wrapper(*args, **kwargs): f"cutoff frequency in Sols.\nData detrended before filtering.\n" f"{Yellow}Generates a new file ending in ``_bps.nc``\n" f"{Green}Example:\n" - f"> MarsFiles 01336.atmos_daily.nc -bps 10 20 -add_trend\n" + f"> MarsFiles 01336.atmos_daily.nc -bps 10 20\n" f"{Nclr}\n\n" ) ) @@ -536,6 +534,34 @@ def wrapper(*args, **kwargs): ) ) +parser.add_argument('-prop', '--prop_tides', action=ExtAction, + ext_content='_prop_tides', + parser=parser, + nargs=2, type=int, + help=( + f"{Yellow}This function is separate distinct from [-tide " + f"--tide_decomp] and therefore does not return total amplitude \n" + f"and phase nor does it work with [-norm --normalize] or [-recon " + f"--reconstruct].\n" + f"For 'diurn' files only.\n" + f"{Nclr}\n" + f"Use fourier decomposition to break down a variable into `kmx` " + f"longitudinal (spatial) harmonics and `tmx` \n" + f"diurnal (time) harmonics. This returns the normalized phases and " + f"amplitudes (not percent) of the \n" + f"propagating tides for a variable.\n" + f"{Yellow}Generates a new file ending in ``_prop_tides.nc``{Nclr}\n" + f"`kmx = 1` for wavenumber 1, `kmx = 2` for wavenumber 2, etc.\n" + f"`tmx = 1` for diurnal tide, `tmx = 2` for semi-diurnal tide, etc.\n" + f"{Green}Example:\n" + f"> MarsFiles 01336.atmos_diurn.nc -prop kmx tmx -incl ps temp\n" + f"> MarsFiles 01336.atmos_diurn.nc -prop 2 2 -incl ps temp\n" + f"{Blue}(extracts the eastward and westward tide components of ps and" + f"temp up to semi-diurnal wavenumber 2)" + f"{Nclr}\n\n" + ) +) + parser.add_argument('-zavg', '--zonal_average', action=ExtAction, ext_content='_zavg', parser=parser, @@ -692,7 +718,7 @@ def wrapper(*args, **kwargs): args.low_pass_temporal, args.band_pass_temporal, args.high_pass_spatial, args.low_pass_spatial, args.band_pass_spatial, args.tide_decomp, args.normalize, - args.regrid_XY_to_match, args.zonal_average] + args.regrid_XY_to_match, args.zonal_average, args.prop_tides] if (all(v is None or v is False for v in all_args) and args.include is not None): @@ -729,8 +755,10 @@ def wrapper(*args, **kwargs): f"{args.tide_decomp_ext}" f"{args.reconstruct_ext}" f"{args.normalize_ext}" + f"{args.prop_tides_ext}" f"{args.regrid_XY_to_match_ext}" f"{args.zonal_average_ext}" + f"{args.concatenate_ext}" ) if args.extension: @@ -810,15 +838,10 @@ def concatenate_files(file_list, full_file_list): ls_end = file_list[-1][18:21] merged_file = f"LegacyGCM_Ls{ls_ini}_Ls{ls_end}.nc" else: - merged_file = full_file_list[0] + output_file_name = full_file_list[0] + merged_file = (f"{output_file_name[:-3]}{out_ext}.nc") - # Delete the files that were concatenated. # Apply the new name created above - for file in full_file_list: - try: - os.remove(file) - except OSError as e: - print(f"Warning: Could not remove {file}: {e}") try: shutil.move(tmp_file, merged_file) @@ -968,32 +991,43 @@ def split_files(file_list, split_dim): ) exit() - if split_dim in ('time', 'areo'): - time_dim = (np.squeeze(fNcdf.variables['time'][:]))[indices] - print(f"time_dim = {time_dim}") - fpath = os.path.dirname(input_file_name) fname = os.path.basename(input_file_name) + + if split_dim in ('time', 'areo'): + time_dim = (np.squeeze(fNcdf.variables['time'][:]))[indices] + print(f"time_dim = {time_dim}\n") + + try: + org_fixed_file = (os.path.normpath(os.path.join(fpath, f"{original_date}.fixed.nc"))) + new_fixed_file = (os.path.normpath(os.path.join(fpath, f"{int(time_dim[0]):05d}.fixed.nc"))) + shutil.copyfile(org_fixed_file, new_fixed_file) + print(f"File {original_date}.fixed.nc copied to {int(time_dim[0]):05d}.fixed.nc.\n") + except FileNotFoundError: + print(f"{Red}No compatible fixed file for {fname} (e.g., {original_date}.fixed.nc) was found in {fpath}{Nclr}\n") + except Exception as e: + print(f"{Red}An error occurred: {e}\n") + if split_dim == 'time': if len(np.atleast_1d(bounds)) < 2: base_name = (f"{int(time_dim):05d}{fname[5:-3]}_nearest_sol" - f"{int(bounds_in[0]):03d}.nc") + f"{int(bounds_in[0]):03d}") output_file_name = (os.path.normpath(os.path.join(fpath, f"{base_name}.nc"))) else: base_name = (f"{int(time_dim[0]):05d}{fname[5:-3]}_sol" - f"{int(bounds_in[0]):05d}_{int(bounds_in[1]):05d}.nc") + f"{int(bounds_in[0]):05d}_{int(bounds_in[1]):05d}") output_file_name = (os.path.normpath(os.path.join(fpath, f"{base_name}.nc"))) elif split_dim =='areo': if len(np.atleast_1d(bounds)) < 2: base_name = (f"{int(time_dim):05d}{fname[5:-3]}_nearest_Ls" - f"{int(bounds_in[0]):03d}.nc") + f"{int(bounds_in[0]):03d}") output_file_name = (os.path.normpath(os.path.join(fpath, f"{base_name}.nc"))) else: base_name = (f"{int(time_dim[0]):05d}{fname[5:-3]}_" - f"Ls{int(bounds_in[0]):03d}_{int(bounds_in[1]):03d}.nc") + f"Ls{int(bounds_in[0]):03d}_{int(bounds_in[1]):03d}") output_file_name = (os.path.normpath(os.path.join(fpath, f"{base_name}.nc"))) split_dim = 'time' @@ -1005,14 +1039,14 @@ def split_files(file_list, split_dim): ] if len(np.atleast_1d(bounds)) < 2: base_name = (f"{original_date}{fname[5:-3]}_nearest_{split_dim}_" - f"{new_bounds[0]}.nc") + f"{new_bounds[0]}") output_file_name = (os.path.normpath(os.path.join(fpath, f"{base_name}.nc"))) else: print(f"{Yellow}bounds = {bounds[0]} {bounds[1]}") print(f"{Yellow}new_bounds = {new_bounds[0]} {new_bounds[1]}") base_name = (f"{original_date}{fname[5:-3]}_{split_dim}_" - f"{new_bounds[0]}_{new_bounds[1]}.nc") + f"{new_bounds[0]}_{new_bounds[1]}") output_file_name = (os.path.normpath(os.path.join(fpath, f"{base_name}.nc"))) elif split_dim == interp_type: @@ -1035,25 +1069,25 @@ def split_files(file_list, split_dim): print(f"{Yellow}bounds = {bounds[0]}") print(f"{Yellow}new_bounds = {new_bounds[0]}") base_name = (f"{original_date}{fname[5:-3]}_nearest_" - f"{new_bounds[0]}.nc") + f"{new_bounds[0]}") output_file_name = (os.path.normpath(os.path.join(fpath, f"{base_name}.nc"))) else: print(f"{Yellow}bounds = {bounds[0]} {bounds[1]}") print(f"{Yellow}new_bounds = {new_bounds[0]} {new_bounds[1]}") base_name = (f"{original_date}{fname[5:-3]}_{new_bounds[0]}_" - f"{new_bounds[1]}.nc") + f"{new_bounds[1]}") output_file_name = (os.path.normpath(os.path.join(fpath, f"{base_name}.nc"))) else: if len(np.atleast_1d(bounds)) < 2: base_name = (f"{original_date}{fname[5:-3]}_nearest_{split_dim}_" - f"{int(bounds[0]):03d}.nc") + f"{int(bounds[0]):03d}") output_file_name = (os.path.normpath(os.path.join(fpath, f"{base_name}.nc"))) else: base_name = (f"{original_date}{fname[5:-3]}_{split_dim}_" - f"{int(bounds[0]):03d}_{int(bounds[1]):03d}.nc") + f"{int(bounds[0]):03d}_{int(bounds[1]):03d}") output_file_name = (os.path.normpath(os.path.join(fpath, f"{base_name}.nc"))) @@ -1765,9 +1799,9 @@ def main(): args.band_pass_spatial): from amescap.Spectral_utils import (zonal_decomposition, zonal_construct, - init_shtools) + import_pyshtools) # Load the module - init_shtools() + import_pyshtools() if args.high_pass_spatial: btype = "high" nk = np.asarray(args.high_pass_spatial).astype(int) @@ -2068,6 +2102,144 @@ def main(): ) fnew.close() + # ------------------------------------------------------------------ + # Propagating Tides + # Richard Urata & R. J. Wilson + # ------------------------------------------------------------------ + elif args.prop_tides: + from amescap.Spectral_utils import extract_diurnal_harmonics + kmx = args.prop_tides[0] + tmx = args.prop_tides[1] + if len(args.prop_tides) != 2: + print(f"{Red}***Error*** prop_tides accepts only two values") + exit() + + for file in file_list: + # Add path unless full path is provided + if not ("/" in file): + input_file_name = os.path.normpath(os.path.join(data_dir, file)) + else: + input_file_name = file + + base_name = os.path.splitext(input_file_name)[0] + output_file_name = os.path.normpath(f"{base_name}{out_ext}.nc") + + fdiurn = Dataset(input_file_name, "r", format="NETCDF4_CLASSIC") + + var_list = filter_vars(fdiurn, args.include) + + # Find time_of_day variable name + tod_name = find_tod_in_diurn(fdiurn) + + target_tod = fdiurn.variables[tod_name][:] + lon = fdiurn.variables["lon"][:] + areo = fdiurn.variables["areo"][:] + numt = areo.shape[0] + + # Define a netcdf object from the netcdf wrapper module + fnew = Ncdf(output_file_name) + # Copy all dims but time_of_day from the old file to the + # new file + + fnew.copy_all_dims_from_Ncfile( + fdiurn, exclude_dim = [tod_name] + ) + # Create new dimension holding the harmonics. We reuse + # the time_of_day name to facilitate. Compatible with + # other routines, but keep in mind this is the harmonic + # number + fnew.add_dim_with_content( + dimension_name = f"diurnal_harmonics", + DATAin = np.arange(1, tmx+1), + longname_txt = "diurnal harmonics", + units_txt = "Diurnal harmonic number", + cart_txt = "tmx" + ) + fnew.add_dim_with_content( + dimension_name = f"lon_harmonics", + DATAin = np.arange(1, kmx+1), + longname_txt = "longitudinal harmonics", + units_txt = "Longitudinal harmonic number", + cart_txt = "kmx" + ) + + # Loop over all variables in the file + for ivar in var_list: + varNcf = fdiurn.variables[ivar] + varIN = varNcf[:] + longname_txt, units_txt = get_longname_unit(fdiurn, ivar) + var_unit = getattr(varNcf, "units", "") + + if (tod_name in varNcf.dimensions and + ivar not in [tod_name, "areo"] and + len(varNcf.shape) > 2): + print(f"{Cyan}Processing: {ivar}{Nclr}") + + # Normalize the data by diurnal mean + norm = np.mean(varIN, axis = 1)[:, np.newaxis, ...] + varIN = varIN/norm + ampe, ampw, phasee, phasew = extract_diurnal_harmonics(kmx, tmx, varIN, target_tod, lon) + + new_dim = list(varNcf.dimensions) + index = new_dim.index(tod_name) + new_dim[index:index+1] = [f"lon_harmonics", f"diurnal_harmonics"] + new_dim.remove("lon") + fnew.log_variable( + f"{ivar}_ampE", + ampe, + new_dim, + f"eastward tidal amplitude for {longname_txt}", + var_unit + ) + fnew.log_variable( + f"{ivar}_phaseE", + phasee, + new_dim, + f"eastward tidal phase for {longname_txt}", + "hr" + ) + fnew.log_variable( + f"{ivar}_ampW", + ampw, + new_dim, + f"westward tidal amplitude for {longname_txt}", + var_unit + ) + fnew.log_variable( + f"{ivar}_phaseW", + phasew, + new_dim, + f"westward tidal phase for {longname_txt}", + "hr" + ) + + elif ivar in ["pfull", "lat", "lon", "phalf", "pk", + "bk", "pstd", "zstd", "zagl", "time"]: + print(f"{Cyan}Copying axis: {ivar}...{Nclr}") + fnew.copy_Ncaxis_with_content(fdiurn.variables[ivar]) + elif ivar in ["areo"]: + print(f"{Cyan}Processing: {ivar}...{Nclr}") + # Create areo variable reflecting the + # new shape + areo_new = np.zeros((areo.shape[0], kmx, tmx, 1)) + # Copy areo + for xx in range(kmx): + for yy in range(tmx): + areo_new[:, xx, yy, :] = areo[:, 0, :] + # Update the dimensions + new_dim = list(varNcf.dimensions) + index = new_dim.index(tod_name) + new_dim[index:index+1] = [f"lon_harmonics", f"diurnal_harmonics"] + # fnew.log_variable(ivar, bareo_new, new_dim, + # longname_txt, units_txt) + fnew.log_variable( + ivar, + areo_new, + new_dim, + longname_txt, + var_unit + ) + fnew.close() # ------------------------------------------------------------------ # Regridding Routine diff --git a/bin/MarsInterp.py b/bin/MarsInterp.py index f4c3dee7..8bcdd27b 100755 --- a/bin/MarsInterp.py +++ b/bin/MarsInterp.py @@ -141,8 +141,8 @@ def wrapper(*args, **kwargs): f"(zstd), or altitude above ground level (zagl).\nWorks on " f"'daily', 'average', and 'diurn' files.\n" f"{Green}Example:\n" - f"> MarsInterp 01336.atmos_average.nc\n" f"> MarsInterp 01336.atmos_average.nc -t pstd\n" + f"> MarsInterp 01336.atmos_average.nc -t pstd -v pstd_default\n" f"{Nclr}\n\n" ) ) @@ -300,7 +300,7 @@ def main(): # Load all of the netcdf files file_list = file_list = [f.name for f in args.input_file] interp_type = args.interp_type # e.g. pstd - custom_level = args.vertical_grid # e.g. p44 + custom_level = args.vertical_grid # e.g. pstd_default grid_out = args.print_grid # Create a namespace with numpy available @@ -322,7 +322,7 @@ def main(): if custom_level: lev_in = eval(f"np.array({custom_level})", namespace) else: - lev_in = np.array(namespace['pstd_default']) + lev_in = eval("np.array(pstd_default)", namespace) # =========================== zstd =========================== elif interp_type == "zstd": @@ -499,8 +499,6 @@ def main(): long_name_txt = getattr(fNcdf.variables[ivar], "long_name", "") units_txt = getattr(fNcdf.variables[ivar], "units", "") - # long_name_txt=fNcdf.variables[ivar].long_name - # units_txt=fNcdf.variables[ivar].units) if not do_diurn: if "tile" in ifile: @@ -538,13 +536,18 @@ def main(): fnew.copy_Ncaxis_with_content(fNcdf.variables[ivar]) else: fnew.copy_Ncvar(fNcdf.variables[ivar]) - + + with Dataset(newname, 'r') as nc_file: + # Print the global attributes of the NetCDF file + print(f"\nGlobal File Attributes for {newname}:") + for attr_name in nc_file.ncattrs(): + print(f" {attr_name}: {getattr(nc_file, attr_name)}") + print("\r ", end="") fNcdf.close() fnew.close() print(f"Completed in {(time.time() - start_time):3f} sec") - - + # ====================================================================== # END OF PROGRAM # ====================================================================== diff --git a/bin/MarsNest.py b/bin/MarsNest.py new file mode 100755 index 00000000..9aa738ef --- /dev/null +++ b/bin/MarsNest.py @@ -0,0 +1,433 @@ +#!/usr/bin/env python +""" +The MarsNest executable does not accept any input. + +The executable requires the presence of grid_spec files in the current directory. +Test grid_spec files can be generated by running the AmesGCM with 0 days as the input. + + +Third-party Requirements: + + * ``numpy`` + * ``argparse`` + * ``netCDF4`` + * ``matplotlib`` + * ``os`` +""" +#============Import modules============== + +# Make print statements appear in color +from amescap.Script_utils import ( + Yellow, Nclr, Green, Red +) + +import numpy as np #numerics module +import matplotlib.pyplot as plt #graphical libraries +from netCDF4 import Dataset +import argparse # parse arguments +import os + +from matplotlib.ticker import MultipleLocator #adjust tic spacing on graphs +from matplotlib.backends.backend_pdf import PdfPages +import matplotlib.path as mpltPath + +#Might be needed on PFE if issues with thr Qt back-end +import matplotlib +matplotlib.use('Agg') + +def debug_wrapper(func): + """ + A decorator that wraps a function with error handling + based on the --debug flag. + + If the --debug flag is set, it prints the full traceback + of any exception that occurs. Otherwise, it prints a + simplified error message. + + :param func: The function to wrap. + :type func: function + :return: The wrapped function. + :rtype: function + :raises Exception: If an error occurs during the function call. + :raises TypeError: If the function is not callable. + :raises ValueError: If the function is not found. + :raises NameError: If the function is not defined. + :raises AttributeError: If the function does not have the + specified attribute. + :raises ImportError: If the function cannot be imported. + :raises RuntimeError: If the function cannot be run. + :raises KeyError: If the function does not have the + specified key. + :raises IndexError: If the function does not have the + specified index. + :raises IOError: If the function cannot be opened. + :raises OSError: If the function cannot be accessed. + :raises EOFError: If the function cannot be read. + :raises MemoryError: If the function cannot be allocated. + :raises OverflowError: If the function cannot be overflowed. + :raises ZeroDivisionError: If the function cannot be divided by zero. + :raises StopIteration: If the function cannot be stopped. + :raises KeyboardInterrupt: If the function cannot be interrupted. + :raises SystemExit: If the function cannot be exited. + :raises AssertionError: If the function cannot be asserted. + """ + + @functools.wraps(func) + def wrapper(*args, **kwargs): + global debug + try: + return func(*args, **kwargs) + except Exception as e: + if debug: + # In debug mode, show the full traceback + print(f"{Red}ERROR in {func.__name__}: {str(e)}{Nclr}") + traceback.print_exc() + else: + # In normal mode, show a clean error message + print(f"{Red}ERROR in {func.__name__}: {str(e)}\nUse " + f"--debug for more information.{Nclr}") + return 1 # Error exit code + return wrapper + +# ====================================================================== +# ARGUMENT PARSER +# ====================================================================== + +parser = argparse.ArgumentParser( + prog=('MarsNest'), + description=( + f"{Yellow}Creates a global map of Mars from grid_spec files\n" + f"including any nests. This tool can be used to orient nest\n" + f"positions. The executable requires the presence of grid_spec\n" + f"files in the current directory. Test grid_spec files can be\n" + f"generated by running the AmesGCM with 0 days as the input.\n" + f"{Nclr}\n\n" + ), + formatter_class=argparse.RawTextHelpFormatter +) + +# Arguments: Only debug + +parser.add_argument('--debug', action='store_true', + help=( + f"Use with any other argument to pass all Python errors and\n" + f"status messages to the screen when running CAP.\n" + f"{Green}Example:\n" + f"> MarsNest --debug" + f"{Nclr}\n\n" + ) + ) + +args = parser.parse_args() +debug = args.debug + + +# ====================================================================== +# DEFINITIONS +# ====================================================================== + + +def extract_grid(lon,lat,varIN,lon_target,lat_target,window_lon,window_lat): + ''' + This function extract a subsample of a larger grid variable. It deals with wrapping around the longitudes on the edges of the grid. + Args: + lon: 1D array of longitude -180/+180 + lat: 1D array of latitudes -90/+90 + varIN: 2D array var2IN(lat,lon) + lon_target: float, center longitude of subgrid (e.g =-136.) + lat_target: float, center latitude of subgrid (e.g =40.) + window_lon: width of window in degree (e.g =10.) + window_lat: heigth of window in degree (e.g. =5.) + Returns: + lon_window: 1D array of truncated longitude (may not be increasing) + lat_window: 1D array of truncated latitudes (may not be increasing) + varIN_out: 2D array of subsample of varIN + ''' + nlat=len(lat) + nlon=len(lon) + Dangle=lon[1]-lon[0] #angle increment + wlon=int(0.5*window_lon/Dangle)# degree window, half side + wlat=int(0.5*window_lat/Dangle)# degree window, hal f side + + #print("(wlat,wlon)= (%i,%i)"%(2*wlat+1,2*wlon+1)) + + ilat_c=int((lat_target-lat[0])/Dangle) #lat index of center of requested tile + ilon_c=int((lon_target-lon[0])/Dangle) #lon index of center of requested tile + + #case where the request is close to the south pole (no issue at the north pole) + # For example --> if lat() has 15 elements, lat[-5:15] will raise an error but lat[0:37] will just return lat[0:15] + if lat_target-window_lat/2lon[-1]: + window_right=int(((lon_target+window_lon/2.)-lon[-1])/Dangle) #number of points outside left of grid + topo_l=varIN[ilat_c-wlat:ilat_c+wlat+1,nlon-(2*wlon-window_right):nlon] #left + topo_r=varIN[ilat_c-wlat:ilat_c+wlat+1,0:window_right+1] #right + varIN_out=np.hstack((topo_l,topo_r)) #combine left and right to wrap around + lon_window=np.append(lon[nlon-(2*wlon-window_right):nlon],lon[0:window_right+1]) + + else: + #---------rest of points------------------ + varIN_out=varIN[ilat_c-wlat:ilat_c+wlat+1,ilon_c-wlon:ilon_c+wlon+1] + lon_window=lon[ilon_c-wlon:ilon_c+wlon+1] + + return lon_window,lat_window,varIN_out + +def extract_path_basename(filename): + ''' + Return the path and basename of a file. If only the filename is provided, assume it is the current directory + Args: + filename: e.g. 'XXXXX.fixed.nc', './history/XXXXX.fixed.nc' or '/home/user/history/XXXXX.fixed.nc' + Returns: + filepath : '/home/user/history/XXXXX.fixed.nc' in all the cases above + basename: XXXXX.fixed.nc in all the cases above + + ***NOTE*** + This routine does not check for file existence and only operates on the provided input string. + ''' + #Get the filename without the path + if '/' in filename or '\\' in filename : + filepath,basename=os.path.split(filename) + else: + filepath=os.getcwd() + basename= filename + + # Is the home ('~') symbol is included, expend the user path + if '~' in filepath:filepath= os.path.expanduser(filepath) + return filepath,basename + +def list_full_paths(input_path): + ''' + List all files in directory. + ''' + return [os.path.join(input_path, file) for file in os.listdir(input_path)] + + + +def find_tile_number(filename): + ''' + Return the line number given a file name. + (Spit grid_spec.nest03.tile8.nc as "grid_spec.nest03" ".tile" "8.nc" and remove ".nc" + ''' + return filename.split('.tile')[1][:-3] + + + +def get_poly_path(gridspec_filename): + ''' + Return the North West and South East corner of a gridscspec file. + Args: + gridspec_filename: full path to a gridscpec file + Return: + poly = a numpy array size (N,(lat,lon) )with the corners + ''' + f=Dataset(gridspec_filename, 'r', format='NETCDF4') + + #Extract 4 corners. Tiles may be rotated to North West is not necessarily at (0,0) + C1=np.asarray([f.variables['grid_lat'][0,0],f.variables['grid_lon'][0,0]]) + C2=np.asarray([f.variables['grid_lat'][0,-1],f.variables['grid_lon'][0,-1]]) + C3=np.asarray([f.variables['grid_lat'][-1,-1],f.variables['grid_lon'][-1,-1]]) + C4=np.asarray([f.variables['grid_lat'][-1,0],f.variables['grid_lon'][-1,0]]) + f.close() + lats=np.array([C1[0],C2[0],C3[0],C4[0]]) + lons=np.array([C1[1],C2[1],C3[1],C4[1]]) + #X,Y couples + return np.array([[C1[1],C1[0]],[C2[1],C2[0]],[C3[1],C3[0]],[C4[1],C4[0]]]) + +def get_corners(poly_path): + ''' + Return the North West and South East corner of a gridscspec file. + Args: + gridspec_filename: full path to a gridscpec file + Return: + NWc,SEc =(lat,lon) of the two corners + ''' + + lons=poly_path[:,0] + lats=poly_path[:,1] + + + #Calculate angular distance between the NW and SE corners + dist_NW=np.abs(lats-lats.max())**2+np.abs(lons-lons.min()) + dist_SE=np.abs(lats-lats.min())**2+np.abs(lons-lons.max()) + + #The corners that shows the min distances are the NW and SE corners + ii_NW=np.argmin(dist_NW) + ii_SE=np.argmin(dist_SE) + + NWc=np.array([lats[ii_NW],lons[ii_NW]]) + SEc=np.array([lats[ii_SE],lons[ii_SE]]) + return NWc,SEc + +def is_child(poly_path_child,poly_path_parent): + ''' + Test if a grid is nested in another grid. + Args: + poly_path_child: a polygon object wit the 4 corners of the child tile + poly_path_parent: a polygon object wit the 4 corners of the parent tile + Returns: + bool: True if all the corners of the child grids are contained within the path defined by the paretn grid + + ''' + path = mpltPath.Path(poly_path_parent) + inside = path.contains_points(poly_path_child) + return inside.all() + + +nskip=2 #skip a couple points in the topography file to speed-up plotting +thick_line=0.5 #line thickness for individual cells + +#TODO development map_dir='/Users/akling/Data/Mars_topo/mars_topo_mola_16.nc' +map_dir='/nobackup/rurata/FMS_MARS_data/mars_topo.nc' + +f=Dataset(map_dir,'r') +lon=f.variables['lon'][::nskip] +lat=f.variables['lat'][::nskip] +topo=f.variables['topo'][::nskip,::nskip] +f.close() + + +#===List of tiles to includes==== + +# ====================================================== +# ARGUMENT PARSER +# ====================================================== +parser = argparse.ArgumentParser(description="""\033[93m Nested grid utility, V%s\033[00m """, + formatter_class=argparse.RawTextHelpFormatter) + +parser.add_argument('input_path', nargs='?', default=os.getcwd(), + help='Optional path to /history directory \n' + '> Usage: MarsNest.py path_to_history/') + + +def main(): + input_path=parser.parse_args().input_path + if input_path=='.':input_path=os.getcwd() #Special case where current dir ('.') is provided + list_dir = list_full_paths(input_path) #Files in directory + avail_grid_specs = [k for k in list_dir if 'grid_spec' in k] + + #List of nest only, note that the directory, e.g. C24_nest/history may contain 'nest' so the test is done on the basename, not the fullpath + avail_grid_nest=[] + for name_tmp in avail_grid_specs: + _,basename = extract_path_basename(name_tmp) #Extract path and basename at two objects + if 'nest' in basename:avail_grid_nest.append(name_tmp) + + outfile=input_path+'/nest_layout.pdf' + + plt.close('all') + with PdfPages(outfile) as pdf: + #Mother grid + plt.figure(figsize=(12,10)) + ax=plt.subplot(111) + plt.contourf(lon,lat,topo,32,cmap='jet') + for name in avail_grid_specs: + + f=Dataset(name, 'r', format='NETCDF4') + grid_lon=f.variables['grid_lon'][:] + grid_lat=f.variables['grid_lat'][:] + f.close() + tileN=find_tile_number(name) + #Get NWc corner + poly=get_poly_path(name) + NWc,_= get_corners(poly) + + x, y = grid_lon,grid_lat + #Filter projection singularities + x[x>0.9*1e+30]=np.NaN + y[y>0.9*1e+30]=np.NaN + _,basename = extract_path_basename(name) #Extract path and basename at two objects + #Plot the grid only for the mother + if 'nest' not in basename: + for i in range(0,grid_lon.shape[0]): + plt.plot(x[:,i], y[:,i], '-k',lw= thick_line) #vertical + plt.plot(x[i,:], y[i,:], '-k',lw= thick_line) #horizontal + #Add tile number + latNWc,lonNWc=NWc + plt.text(lonNWc+4,latNWc-4,tileN,fontsize='18',color='r') + + plt.plot(x[0,:], y[0,:], 'r',lw=3 )#S + plt.plot(x[:,-1], y[:,-1], 'r',lw=3 )#E + plt.plot(x[-1,:], y[-1,:], 'r',lw=3 )#N + plt.plot(x[:,0], y[:,0], 'r',lw=3 )#W + + ax.xaxis.set_major_locator(MultipleLocator(30)) + ax.xaxis.set_minor_locator(MultipleLocator(10)) + ax.yaxis.set_major_locator(MultipleLocator(30)) + ax.yaxis.set_minor_locator(MultipleLocator(10)) + ax.set_xlim([0,360]) + ax.set_ylim([-90,90]) + pdf.savefig() # saves the current figure into a pdf page + plt.close() + + + #====Add one page for each nest===== + for name in avail_grid_nest: + _,basename = extract_path_basename(name) #Extract path and basename at two objects + plt.figure(figsize=(12,10)) + ax=plt.subplot(111) + f=Dataset(name, 'r', format='NETCDF4') + grid_lon=f.variables['grid_lon'][:] + grid_lat=f.variables['grid_lat'][:] + f.close() + + #Get grid center lat/lon and angular width of the tile + N_tile=grid_lon.shape[0] + lon_cent=grid_lon[N_tile//2,N_tile//2] + lat_cent=grid_lat[N_tile//2,N_tile//2] + Dlon=(grid_lon.max()-grid_lon.min()) + + #Get high res topo + lon_tmp,lat_tmp,topo_tmp=extract_grid(lon=lon,lat=lat,varIN=topo,lon_target=lon_cent,lat_target=lat_cent,window_lon=Dlon,window_lat=Dlon) + + + poly_parent=get_poly_path(name) + NWc,_= get_corners(poly_parent) + x, y = grid_lon,grid_lat + #Filter projection singularities + x[x>0.9*1e+30]=np.NaN + y[y>0.9*1e+30]=np.NaN + + tileN=find_tile_number(name) + + plt.contourf(lon_tmp,lat_tmp,topo_tmp,32,cmap='jet') + plt.plot(x[0,:], y[0,:], 'r',lw=3 )#S + plt.plot(x[:,-1], y[:,-1], 'r',lw=3 )#E + plt.plot(x[-1,:], y[-1,:], 'r',lw=3 )#N + plt.plot(x[:,0], y[:,0], 'r',lw=3 )#W + + #Add tile number + latNWc,lonNWc=NWc + plt.text(lonNWc+0.5,latNWc-0.5,tileN,fontsize='18',color='r') + #Plot the nested nest + for name2 in avail_grid_nest: + poly_child=get_poly_path(name2) + if name2!=name and is_child(poly_child,poly_parent): + f=Dataset(name2, 'r', format='NETCDF4') + grid_lon=f.variables['grid_lon'][:] + grid_lat=f.variables['grid_lat'][:] + f.close() + x, y = grid_lon,grid_lat + #Filter projection singularities + x[x>0.9*1e+30]=np.NaN + y[y>0.9*1e+30]=np.NaN + plt.plot(x[0,:], y[0,:], 'r',lw=3 )#S + plt.plot(x[:,-1], y[:,-1], 'r',lw=3 )#E + plt.plot(x[-1,:], y[-1,:], 'r',lw=3 )#N + plt.plot(x[:,0], y[:,0], 'r',lw=3 )#W + + pdf.savefig() # saves the current figure into a pdf page + plt.close() + + + print(outfile +' was created') + +if __name__ == '__main__': + main() diff --git a/bin/MarsVars.py b/bin/MarsVars.py index fbcf552d..6261c194 100755 --- a/bin/MarsVars.py +++ b/bin/MarsVars.py @@ -254,6 +254,18 @@ def wrapper(*args, **kwargs): ['Vg_sed', 'w'], ['pfull', 'pstd', 'zstd', 'zagl'] ], + 'dustref_per_pa': [ + "Dust opacity per Pascal [dustref/delp]", + 'op/Pa', + ['dustref', 'delp'], + ['pfull'] + ], + 'dustref_per_km': [ + "Dust opacity per kilometer [dustref/delz]", + 'op/km', + ['dustref', 'delz'], + ['pfull'] + ], 'wdir': [ "Wind direction", 'degree', @@ -299,7 +311,7 @@ def wrapper(*args, **kwargs): 'msf': [ "Mass stream function", '1.e8 kg/s', - ['vcomp'], + ['vcomp', 'temp'], ['pstd', 'zstd', 'zagl'] ], 'mx': [ @@ -1173,6 +1185,59 @@ def compute_w_net(Vg, wvar): return w_net + +# ===================================================================== +def compute_dustref_per_pa(dustref, delp): + """ + Computes visible dust opacity per Pascal from dustref and delp. + dustref is visible dust opacity per level (model layer). + + opacity per Pa = opacity per layer / layer thickness in Pa:: + + dustref_per_pa = dustref/delp + + [Courtney Batterson, 2025] + + :param dustref: Visible dust opacity [op/model layer] + :type dustref: array [time, lev, lat, lon] + :param delp: Layer thickness [Pa] + :type delp: array [time, lev, lat, lon] + :return: `dustref_per_pa` Visible dust opacity [op/Pa] + :rtype: array [time, lev, lat, lon] + :raises ValueError: If the input dimensions are not compatible + :raises TypeError: If the input types are not compatible + :raises Exception: If any other error occurs + """ + + dustref_per_pa = dustref/delp + return dustref_per_pa + +# ===================================================================== +def compute_dustref_per_z(dustref, delz): + """ + Computes visible dust opacity per kilometer from dustref and delz. + dustref is visible dust opacity per level (model layer). + + opacity per km = opacity per layer / layer thickness in m * 1000:: + + dustref_per_z = dustref/delz*1000 + + [Courtney Batterson, 2025] + + :param dustref: Visible dust opacity [op/model layer] + :type dustref: array [time, lev, lat, lon] + :param delz: Layer thickness [m] + :type delz: array [time, lev, lat, lon] + :return: `dustref_per_z` Visible dust opacity [op/km] + :rtype: array [time, lev, lat, lon] + :raises ValueError: If the input dimensions are not compatible + :raises TypeError: If the input types are not compatible + :raises Exception: If any other error occurs + """ + + dustref_per_z = dustref/delz*1000 + return dustref_per_z + # ===================================================================== def compute_theta(p_3D, ps, T, f_type): """ @@ -1906,9 +1971,9 @@ def process_add_variables(file_name, add_list, master_list, debug=False): f"file:{Nclr}") for var, actual_var in already_in_file: if var == actual_var: - print(f"{Yellow} - {var}{Nclr}") + print(f"{Yellow} - {var}{Nclr}") else: - print(f"{Yellow} - {var} (as '{actual_var}'){Nclr}") + print(f"{Yellow} - {var} (as '{actual_var}'){Nclr}") return @@ -2101,6 +2166,16 @@ def add_with_dependencies(var): wvar = f.variables["w"][:] OUT = compute_w_net(Vg, wvar) + if var == "dustref_per_pa": + dustref = f.variables["dustref"][:] + delp = f.variables["delp"][:] + OUT = compute_dustref_per_pa(dustref, delp) + + if var == "dustref_per_z": + dustref = f.variables["dustref"][:] + delz = f.variables["delz"][:] + OUT = compute_dustref_per_z(dustref, delz) + if var == "pfull3D": OUT = p_3D @@ -2709,6 +2784,14 @@ def __init__(self, name): if interp_type == "pfull": ak, bk = ak_bk_loader(f) + else: + print(f"{Red}column integration error: the file " + f"{input_file} does not have the requisite variables " + f"(ak and bk), please use a non-column integrated file " + f"to use this function (e.g., atmos_average.nc and not " + f"of atmos_average_pstd.nc).{Nclr}") + f.close() + continue # Use check_variable_exists instead of direct key lookup if not check_variable_exists(icol, f.variables.keys()): @@ -2800,7 +2883,8 @@ def __init__(self, name): lname_text = getattr(var_Ncdf, "long_name", "") unit_text = getattr(var_Ncdf, "units", "") cart_text = getattr(var_Ncdf, "cartesian_axis", "") - + datatype=var_Ncdf.dtype + if args.rename: name_text = args.rename if args.longname: @@ -2809,16 +2893,16 @@ def __init__(self, name): unit_text = args.unit if args.multiply: vals *= args.multiply - + if cart_text == "": Log.log_variable( - name_text, vals, dim_out, lname_text, unit_text + name_text, vals, dim_out, lname_text, unit_text, datatype ) else: Log.log_axis1D( - name_text, vals, dim_out, lname_text, unit_text, cart_text + name_text, vals, dim_out, lname_text, unit_text, cart_text, datatype ) - + # Close files to release handles f_IN.close() Log.close() diff --git a/docs/NAS_INSTRUCTIONS.txt b/docs/NAS_INSTRUCTIONS.txt deleted file mode 100644 index c57212f9..00000000 --- a/docs/NAS_INSTRUCTIONS.txt +++ /dev/null @@ -1,31 +0,0 @@ -#!/bin/bash - -#================================================== -# AmesCAP Install on NAS -# Updated: June 2023 -# Last Editor: Kling -#================================================== - -# 1. Remove the AmesCAP virtual environment folder entirely: -rm -r AmesCAP # or whatever the name is for you - -# 2. Load the necessary modules: -module purge -module load python3/3.11.5 - -# 3. (Re)create the virtual environment -python3.11 -m venv --system-site-packages AmesCAP - -# 4. Activate the virtual environment using the appropriate syntax for your shell: -#______________________________________________________________________________________________ -# CSH/TSCH | BASH -#________________________________________________|_____________________________________________ -source AmesCAP/bin/activate.csh | source AmesCAP/bin/activate -#---------------------------------------------------------------------------------------------- - - -# 5. Install CAP in the virtual environment: -pip install git+https://github.com/NASA-Planetary-Science/AmesCAP.git - -# 6. Deactivate CAP with: -deactivate diff --git a/docs/source/description.rst b/docs/source/description.rst index f2c4f5da..09759a52 100644 --- a/docs/source/description.rst +++ b/docs/source/description.rst @@ -143,6 +143,7 @@ Data Transformation Functions ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ * Perform **tidal analysis** on diurnal composite files: ``[-tide --tide_decomp]`` +* Perform **propagating tideal analysis** on diurnal composite files: ``[-prop --prop_tides]`` * Apply **temporal filters** to time-varying fields: * Low pass: ``[-lpt --low_pass_temporal]`` @@ -299,8 +300,6 @@ Open ``~/.amesgcm_profile`` with any text editor to see customizable grid defini 1.0e-02, 5.0e-03, 3.0e-03, 5.0e-04, 3.0e-04, 1.0e-04, 5.0e-05, 3.0e-05, 1.0e-05] - phalf_mb=[50] - Use your custom grid with the ``[-v --vertical_grid]`` argument: .. code-block:: bash diff --git a/docs/source/installation.rst b/docs/source/installation.rst index 508d2c8a..2bd98770 100644 --- a/docs/source/installation.rst +++ b/docs/source/installation.rst @@ -469,7 +469,7 @@ The pip installation method is less recommended as it requires manual installati # Activate your virtual environment # Install CAP with spectral analysis support - pip install "amescap[spectral] @ git+https://github.com/NASA-Planetary-Science/AmesCAP.git@pyshtools" + pip install "amescap[spectral] @ git+https://github.com/NASA-Planetary-Science/AmesCAP.git" # Don't forget to copy the profile file to your home directory cp amescap/mars_templates/amescap_profile ~/.amescap_profile diff --git a/pyproject.toml b/pyproject.toml index 5ecf683c..ffb25a8c 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -42,6 +42,7 @@ MarsVars = "bin.MarsVars:main" MarsFiles = "bin.MarsFiles:main" MarsFormat = "bin.MarsFormat:main" MarsCalendar = "bin.MarsCalendar:main" +MarsNest = "bin.MarsNest:main" [tool.setuptools] packages = ["amescap", "bin"] diff --git a/tests/test_marsfiles.py b/tests/test_marsfiles.py index b962dc74..1eebbecd 100644 --- a/tests/test_marsfiles.py +++ b/tests/test_marsfiles.py @@ -659,7 +659,26 @@ def test_tide_decomposition_with_reconstruct(self): nc.close() print("Include argument succeeded") - + + def test_prop_tides(self): + """Test propagating tidal analysis on diurn file""" + + result = self.run_mars_files(['01336.atmos_diurn.nc', '-prop', '2', '2', '-incl', 'ps', 'temp']) + + # Check for successful execution + self.assertEqual(result.returncode, 0, "Propagating tide analysis command failed") + + # Check that output file was created + output_file = self.check_file_exists('01336.atmos_diurn_prop_tides.nc') + + # Verify that the output file has amplitude and phase variables + self.verify_netcdf_has_variable(output_file, 'ps_ampE') + self.verify_netcdf_has_variable(output_file, 'ps_phaseE') + self.verify_netcdf_has_variable(output_file, 'temp_ampW') + self.verify_netcdf_has_variable(output_file, 'temp_phaseW') + + print("Propagating tide analysis succeeded") + def test_regrid(self): """Test regridding operation""" result = self.run_mars_files([