diff --git a/doc/source/notebooks/GRACE-Spatial-Error.ipynb b/doc/source/notebooks/GRACE-Spatial-Error.ipynb index 3177750..a40e8f7 100644 --- a/doc/source/notebooks/GRACE-Spatial-Error.ipynb +++ b/doc/source/notebooks/GRACE-Spatial-Error.ipynb @@ -306,9 +306,9 @@ "mm = np.arange(0, MMAX + 1)\n", "PLM2 = PLM[:, mm, :] ** 2\n", "# Calculating cos(m*phi)^2 and sin(m*phi)^2\n", - "phi = np.radians(grid.lon[np.newaxis, :])\n", - "ccos = np.cos(np.dot(mm[:, np.newaxis], phi)) ** 2\n", - "ssin = np.sin(np.dot(mm[:, np.newaxis], phi)) ** 2\n", + "phi = np.radians(grid.lon)\n", + "mp = np.einsum('m...,p...->mp...', mm, phi)\n", + "m_phi2 = (1.0 - 1j) * (np.exp(-2j * mp) + np.exp(2j * mp) + 2j) / 4.0\n", "\n", "# read load love numbers file\n", "# PREM outputs from Han and Wahr (1995)\n", @@ -346,7 +346,7 @@ "\n", "# dfactor is the degree dependent coefficients\n", "# for converting to spherical harmonic output units\n", - "factors = gravtk.units(lmax=LMAX).harmonic(hl, kl, ll).mmwe\n", + "factors = gravtk.units(lmax=LMAX).harmonic(hl, kl, ll)\n", "# mmwe, millimeters water equivalent\n", "dfactor = factors.get('mmwe')\n", "# units strings for output plots\n", @@ -388,17 +388,14 @@ "delta_Ylms = delta_Ylms.convolve(dfactor * wt)\n", "# smooth harmonics and convert to output units\n", "YLM2 = delta_Ylms.power(2.0).scale(1.0 / nsmth)\n", + "\n", "# Calculate fourier coefficients\n", - "d_cos = np.zeros((MMAX + 1, nlat)) # [m,th]\n", - "d_sin = np.zeros((MMAX + 1, nlat)) # [m,th]\n", - "# Calculating delta spatial values\n", - "for k in range(0, nlat):\n", - " # summation over all spherical harmonic degrees\n", - " d_cos[:, k] = np.sum(PLM2[:, :, k] * YLM2.clm, axis=0)\n", - " d_sin[:, k] = np.sum(PLM2[:, :, k] * YLM2.slm, axis=0)\n", + "# summation over all spherical harmonic degrees\n", + "pconv2 = np.einsum('lmh...,lm...->mh...', PLM2, YLM2.ilm)\n", "\n", - "# Multiplying by c/s(phi#m) to get spatial maps (lon,lat)\n", - "grid.data = np.sqrt(np.dot(ccos.T, d_cos) + np.dot(ssin.T, d_sin)).T\n", + "# Multiplying by c/s(phi#m) to get spatial error map\n", + "# take the square root and drop imaginary component\n", + "grid.data = np.sqrt(np.einsum('mp...,mh...->hp...', m_phi2, pconv2)).real\n", "grid.mask = np.zeros_like(grid.data, dtype=bool)" ] }, @@ -514,11 +511,9 @@ } ], "metadata": { - "interpreter": { - "hash": "31f2aee4e71d21fbe5cf8b01ff0e069b9275f58929596ceb00d14d90e3e16cd6" - }, "kernelspec": { - "display_name": "Python 3.8.10 64-bit", + "display_name": "py13", + "language": "python", "name": "python3" }, "language_info": { @@ -531,7 +526,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.10.6" + "version": "3.13.0" } }, "nbformat": 4, diff --git a/gravity_toolkit/grace_date.py b/gravity_toolkit/grace_date.py index 9388e30..e8b4e18 100644 --- a/gravity_toolkit/grace_date.py +++ b/gravity_toolkit/grace_date.py @@ -47,6 +47,7 @@ UPDATE HISTORY: Updated 08/2026: added verbosity argument to print debugging information + use numpy datetime arrays to save start and end dates Updated 05/2023: use pathlib to define and operate on paths Updated 03/2023: use f-strings for formatting output date lines added regex formatting for CNES GRGS harmonics @@ -110,7 +111,9 @@ import pathlib import argparse import numpy as np +import collections import gravity_toolkit as gravtk +from datetime import datetime # PURPOSE: parses GRACE/GRACE-FO data files and assigns month numbers @@ -177,10 +180,8 @@ def grace_date(base_dir, PROC='', DREL='', DSET='', OUTPUT=True, MODE=0o775): n_files = len(input_files) # define date variables - start_yr = np.zeros((n_files)) # year start date - end_yr = np.zeros((n_files)) # year end date - start_day = np.zeros((n_files)) # day number start date - end_day = np.zeros((n_files)) # day number end date + start_date = np.zeros((n_files), dtype='datetime64[s]') # start date + end_date = np.zeros((n_files), dtype='datetime64[s]') # end date mid_day = np.zeros((n_files)) # mid-month day tot_days = np.zeros((n_files)) # number of days since Jan 2002 tdec = np.zeros((n_files)) # date in decimal form @@ -188,74 +189,47 @@ def grace_date(base_dir, PROC='', DREL='', DSET='', OUTPUT=True, MODE=0o775): # for each data file for t, infile in enumerate(input_files): - if PROC in ( - 'GRAZ', - 'Swarm', - ): + if PROC in ('GRAZ', 'Swarm'): # get date lists for the start and end of fields - start_date, end_date = gravtk.time.parse_gfc_file( + start_date[t], end_date[t] = gravtk.time.parse_gfc_file( infile, PROC, DSET ) - # start and end year - start_yr[t] = np.float64(start_date[0]) - end_yr[t] = np.float64(end_date[0]) - # number of days in each month for the calendar year - dpm = gravtk.time.calendar_days(start_yr[t]) - # start and end day of the year - start_day[t] = ( - np.sum(dpm[: start_date[1] - 1]) - + start_date[2] - + start_date[3] / 24.0 - + start_date[4] / 1440.0 - + start_date[5] / 86400.0 - ) - end_day[t] = ( - np.sum(dpm[: end_date[1] - 1]) - + end_date[2] - + end_date[3] / 24.0 - + end_date[4] / 1440.0 - + end_date[5] / 86400.0 - ) else: # get date lists for the start and end of fields - start_date, end_date = gravtk.time.parse_grace_file(infile) - # start and end year - start_yr[t] = np.float64(start_date[0]) - end_yr[t] = np.float64(end_date[0]) - # start and end day of the year - start_day[t] = np.float64(start_date[1]) - end_day[t] = np.float64(end_date[1]) + start_date[t], end_date[t] = gravtk.time.parse_grace_file(infile) + # start year and day of the year + start_struct = start_date[t].astype(datetime).timetuple() + start_yr = start_struct.tm_year + start_day = start_struct.tm_yday + # end year and day of the year + end_struct = end_date[t].astype(datetime).timetuple() + end_yr = end_struct.tm_year + end_day = end_struct.tm_yday # number of days in the starting year for leap and standard years - dpy = gravtk.time.calendar_days(start_yr[t]).sum() + dpy = gravtk.time.calendar_days(start_yr).sum() # end date taking into account measurements taken on different years - end_cyclic = (end_yr[t] - start_yr[t]) * dpy + end_day[t] + end_cyclic = (end_yr - start_yr) * dpy + end_day # calculate mid-month value - mid_day[t] = np.mean([start_day[t], end_cyclic]) + mid_day[t] = np.mean([start_day, end_cyclic]) # calculate Modified Julian Day from start_yr and mid_day MJD = gravtk.time.convert_calendar_dates( - start_yr[t], 1.0, mid_day[t], epoch=(1858, 11, 17, 0, 0, 0) + start_yr, 1.0, mid_day[t], epoch=(1858, 11, 17, 0, 0, 0) ) # convert from Modified Julian Days to calendar dates cal_date = gravtk.time.convert_julian(MJD + 2400000.5) # Calculating the mid-month date in decimal form - tdec[t] = start_yr[t] + mid_day[t] / dpy + tdec[t] = start_yr + mid_day[t] / dpy # Calculation of total days since start of campaign - count = 0 - n_yrs = np.int64(start_yr[t] - 2002) - # for each of the GRACE years up to the file year - for iyr in range(n_yrs): - # year - year = 2002 + iyr - # add all days from prior years to count - # number of days in year i (if leap year or standard year) - count += gravtk.time.calendar_days(year).sum() - - # calculating the total number of days since 2002 - tot_days[t] = np.mean([count + start_day[t], count + end_cyclic]) + # add all days from prior years to count + # number of days in year i (if leap year or standard year) + count = np.sum( + [gravtk.time.calendar_days(y) for y in range(2002, start_yr)] + ) + tot_days[t] = count + mid_day[t] # Calculates the month number (or 10-day number for CNES RL01,RL02) if (PROC == 'CNES') and (DREL in ('RL01', 'RL02')): @@ -285,18 +259,33 @@ def grace_date(base_dir, PROC='', DREL='', DSET='', OUTPUT=True, MODE=0o775): print('{0} {1:>10} {2:>11} {3:>10} {4:>13}'.format(*args), file=fid) # create python dictionary mapping input file names with GRACE months - grace_files = {} + grace_files = collections.OrderedDict() + # add attributes + grace_files.attrs = { + 'start_date': start_date, + 'end_date': end_date, + 'time_decimal': tdec, + } # for each data file for t, infile in enumerate(input_files): # add file to python dictionary mapped to GRACE/GRACE-FO month grace_files[mon[t]] = grace_dir.joinpath(infile) # print to GRACE dates ascii file (NOTE: tot_days will be rounded) if OUTPUT: + # start year and day of the year + start_struct = start_date[t].astype(datetime).timetuple() + start_yr = start_struct.tm_year + start_day = start_struct.tm_yday + # end year and day of the year + end_struct = end_date[t].astype(datetime).timetuple() + end_yr = end_struct.tm_year + end_day = end_struct.tm_yday + # print to GRACE dates ascii file print( ( f'{tdec[t]:13.8f} {mon[t]:03d} ' - f'{start_yr[t]:8.0f} {start_day[t]:03.0f} ' - f'{end_yr[t]:8.0f} {end_day[t]:03.0f} ' + f'{start_yr:8.0f} {start_day:03.0f} ' + f'{end_yr:8.0f} {end_day:03.0f} ' f'{tot_days[t]:8.0f}' ), file=fid, diff --git a/gravity_toolkit/grace_input_months.py b/gravity_toolkit/grace_input_months.py index 4df1db5..4b13990 100644 --- a/gravity_toolkit/grace_input_months.py +++ b/gravity_toolkit/grace_input_months.py @@ -110,6 +110,7 @@ UPDATE HISTORY: Updated 08/2026: use upstream file logger for verbose output + save start and end date of files as numpy datetime64[s] arrays Updated 10/2023: standardize ocean model for UCI degree 1 coefficients Updated 05/2023: use pathlib to define and operate on paths Updated 04/2023: use release-03 GFZ GravIS SLR and geocenter files @@ -363,7 +364,7 @@ def grace_input_months( n_cons = len(months) # Initializing input data matrices - grace_Ylms = {} + grace_Ylms = collections.OrderedDict() grace_Ylms['clm'] = np.zeros((LMAX + 1, MMAX + 1, n_cons)) grace_Ylms['slm'] = np.zeros((LMAX + 1, MMAX + 1, n_cons)) grace_Ylms['eclm'] = np.zeros((LMAX + 1, MMAX + 1, n_cons)) @@ -373,6 +374,11 @@ def grace_input_months( # output dimensions grace_Ylms['l'] = np.arange(LMAX + 1) grace_Ylms['m'] = np.arange(MMAX + 1) + # datetime arrays for start and end dates + grace_Ylms.attrs = { + 'start_date': np.zeros((n_cons), dtype='datetime64[s]'), + 'end_date': np.zeros((n_cons), dtype='datetime64[s]'), + } # attributes for processing run attributes = collections.OrderedDict() @@ -383,11 +389,14 @@ def grace_input_months( grace_files = grace_date( base_dir, PROC=PROC, DREL=DREL, DSET=DSET, OUTPUT=False ) + # get the list of GRACE/GRACE-FO months for the input files + file_months = list(grace_files.keys()) # importing data from GRACE/GRACE-FO files for i, grace_month in enumerate(months): # read spherical harmonic data products infile = grace_files[grace_month] + j = file_months.index(grace_month) # log input file if debugging logger.debug(f'Reading file {i:d}: {str(infile)}') # read GRACE/GRACE-FO/Swarm file @@ -408,6 +417,9 @@ def grace_input_months( # copy date variables grace_Ylms['time'][i] = np.copy(Ylms['time']) grace_Ylms['month'][i] = np.int64(grace_month) + # copy start and end dates from grace_files attributes + grace_Ylms.attrs['start_date'][i] = grace_files.attrs['start_date'][j] + grace_Ylms.attrs['end_date'][i] = grace_files.attrs['end_date'][j] # copy input file basename attributes['lineage'][i] = pathlib.Path(infile).stem diff --git a/gravity_toolkit/mapping/plot_AIS_GrIS_maps.py b/gravity_toolkit/mapping/plot_AIS_GrIS_maps.py index eaaaa51..5fc511d 100644 --- a/gravity_toolkit/mapping/plot_AIS_GrIS_maps.py +++ b/gravity_toolkit/mapping/plot_AIS_GrIS_maps.py @@ -39,6 +39,7 @@ https://pypi.python.org/pypi/GDAL/ UPDATE HISTORY: + Updated 08/2026: use upstream file logger for verbose output Updated 10/2023: increase size of scale bar for both maps Updated 05/2023: use pathlib to define and operate on paths added option to set the input variable names or column order @@ -236,23 +237,28 @@ # PURPOSE: keep track of threads def info(args): - logging.info(pathlib.Path(sys.argv[0]).name) - logging.info(args) - logging.info(f'module name: {__name__}') + # get logger + logger = logging.getLogger(__name__) + logger.info(pathlib.Path(sys.argv[0]).name) + logger.info(args) + logger.info(f'module name: {__name__}') if hasattr(os, 'getppid'): - logging.info(f'parent process: {os.getppid():d}') - logging.info(f'process id: {os.getpid():d}') + logger.info(f'parent process: {os.getppid():d}') + logger.info(f'process id: {os.getpid():d}') # PURPOSE: plot Rignot 2012 drainage basin polylines def plot_rignot_basins(ax, base_dir, HEM, projection): + # get logger + logger = logging.getLogger(__name__) region_directory = base_dir.joinpath(*region_dir) + # regional drainage basin files + region_files = region_filename[HEM] # for each region for reg in region_title: # read the regional polylines - region_file = region_directory.joinpath( - region_filename[HEM].format(reg) - ) + region_file = region_directory.joinpath(region_files.format(reg)) + logger.debug(str(region_file)) region_ll = np.loadtxt(region_file, dtype=region_dtype) # converting region lat/lon into plot coordinates points = projection.transform_points( @@ -263,9 +269,11 @@ def plot_rignot_basins(ax, base_dir, HEM, projection): # PURPOSE: plot Greenland and Antarctic drainage basins from IMBIE2 (Mouginot) def plot_IMBIE2_basins(ax, base_dir, HEM): + # get logger + logger = logging.getLogger(__name__) # read drainage basin polylines from shapefile (using splat operator) basin_shapefile = base_dir.joinpath(*IMBIE_basin_file[HEM]) - logging.debug(str(basin_shapefile)) + logger.debug(str(basin_shapefile)) shape_input = shapefile.Reader(str(basin_shapefile)) shape_entities = shape_input.shapes() shape_attributes = shape_input.records() @@ -310,10 +318,12 @@ def plot_IMBIE2_basins(ax, base_dir, HEM): # PURPOSE: plot Antarctic drainage sub-basins from IMBIE-2 (Mouginot) def plot_IMBIE2_subbasins(ax, base_dir): + # get logger + logger = logging.getLogger(__name__) # read drainage basin polylines from shapefile (using splat operator) IMBIE_subbasin_file = ['Basins_20Oct2016_v1.7', 'Basins_v1.7.shp'] basin_shapefile = base_dir.joinpath('masks', *IMBIE_subbasin_file) - logging.debug(str(basin_shapefile)) + logger.debug(str(basin_shapefile)) shape_input = shapefile.Reader(str(basin_shapefile)) shape_entities = shape_input.shapes() shape_attributes = shape_input.records() @@ -336,7 +346,10 @@ def plot_IMBIE2_subbasins(ax, base_dir): # PURPOSE: plot Greenland and Antarctic grounded ice delineation def plot_grounded_ice(ax, base_dir, HEM, START=1): + # get logger + logger = logging.getLogger(__name__) grounded_ice_shape_file = base_dir.joinpath(*coast_file[HEM]) + logger.debug(str(grounded_ice_shape_file)) shape_input = shapefile.Reader(str(grounded_ice_shape_file)) shape_entities = shape_input.shapes() shape_attributes = shape_input.records() @@ -359,6 +372,8 @@ def plot_grounded_ice(ax, base_dir, HEM, START=1): # PURPOSE plot coastlines and islands (GSHHS with G250 Greenland) def plot_coastline(ax, base_dir): + # get logger + logger = logging.getLogger(__name__) # read the coastline shape file coastline_dir = base_dir.joinpath('masks', 'G250') coastline_shape_files = [] @@ -366,7 +381,7 @@ def plot_coastline(ax, base_dir): coastline_shape_files.append('greenland_coastline_islands.shp') for fi, S in zip(coastline_shape_files, [1000, 200]): coast_shapefile = coastline_dir.joinpath(fi) - logging.debug(str(coast_shapefile)) + logger.debug(str(coast_shapefile)) shape_input = shapefile.Reader(str(coast_shapefile)) shape_entities = shape_input.shapes() # for each entity within the shapefile @@ -378,8 +393,12 @@ def plot_coastline(ax, base_dir): # PURPOSE: plot MODIS mosaic of Antarctica and Greenland as background image def plot_image_mosaic(ax, base_dir, HEM, MASKED=True): + # get logger + logger = logging.getLogger(__name__) # read MODIS mosaic of Antarctica and Greenland - ds = osgeo.gdal.Open(base_dir.joinpath(*image_file[HEM])) + image_geotiff_file = base_dir.joinpath(*image_file[HEM]) + logger.debug(str(image_geotiff_file)) + ds = osgeo.gdal.Open(str(image_geotiff_file)) # get dimensions xsize = ds.RasterXSize ysize = ds.RasterYSize @@ -433,12 +452,31 @@ def add_plot_scale(ax, X, Y, dx, dy, masked, fc1='w', fc2='k'): Y + 3.2 * dy, ] ax.fill( - [x1, x2, x2, x1, x1], [y1, y1, y2, y2, y1], fc1, alpha=0.5, zorder=4 + [x1, x2, x2, x1, x1], + [y1, y1, y2, y2, y1], + fc1, + alpha=0.5, + zorder=4, ) for i, c in enumerate([fc1, fc2, fc1, fc2]): - x1, x2, y1, y2 = [X + 0.25 * i * dx, X + 0.25 * (i + 1) * dx, Y, Y + dy] - ax.fill([x1, x2, x2, x1, x1], [y1, y1, y2, y2, y1], c, zorder=5) - ax.plot([X, X + dx, X + dx, X, X], [Y, Y, Y + dy, Y + dy, Y], fc2, zorder=6) + x1, x2, y1, y2 = [ + X + 0.25 * i * dx, + X + 0.25 * (i + 1) * dx, + Y, + Y + dy, + ] + ax.fill( + [x1, x2, x2, x1, x1], + [y1, y1, y2, y2, y1], + c, + zorder=5, + ) + ax.plot( + [X, X + dx, X + dx, X, X], + [Y, Y, Y + dy, Y + dy, Y], + fc2, + zorder=6, + ) for i in range(3): ax.plot( [X + 0.5 * i * dx, X + 0.5 * i * dx], @@ -504,6 +542,8 @@ def plot_grid( FIGURE_DPI=None, MODE=0o775, ): + # get logger + logger = logging.getLogger(__name__) # extend list if a single format was entered for all files if len(DATAFORM) < len(FILENAMES): DATAFORM = DATAFORM * len(FILENAMES) @@ -604,7 +644,7 @@ def plot_grid( for i, ax in ax1.items(): # set hemisphere flag HEM = hem_flag[i] - logging.info(f'Hemisphere: {HEM}') + logger.info(f'Hemisphere: {HEM}') # x and y limits for region xlimits, ylimits = (region_xlimit[HEM], region_ylimit[HEM]) # plot image background (MODIS Mosaic of Antarctica or Greenland) @@ -865,7 +905,7 @@ def plot_grid( # create output directory if non-existent FIGURE_FILE.parent.mkdir(mode=MODE, parents=True, exist_ok=True) # save to file - logging.info(str(FIGURE_FILE)) + logger.info(str(FIGURE_FILE)) plt.savefig( FIGURE_FILE, metadata={'Title': pathlib.Path(sys.argv[0]).name}, @@ -1124,7 +1164,9 @@ def main(): # create logger loglevels = [logging.CRITICAL, logging.INFO, logging.DEBUG] - logging.basicConfig(level=loglevels[args.verbose]) + logger = gravtk.utilities.build_logger( + __name__, level=loglevels[args.verbose] + ) # try to run the analysis with listed parameters try: @@ -1167,8 +1209,8 @@ def main(): # if there has been an error exception # print the type, value, and stack trace of the # current exception being handled - logging.critical(f'process id {os.getpid():d} failed') - logging.error(traceback.format_exc()) + logger.critical(f'process id {os.getpid():d} failed') + logger.error(traceback.format_exc()) # run main program diff --git a/gravity_toolkit/mapping/plot_AIS_grid_3maps.py b/gravity_toolkit/mapping/plot_AIS_grid_3maps.py index 701e468..d940372 100644 --- a/gravity_toolkit/mapping/plot_AIS_grid_3maps.py +++ b/gravity_toolkit/mapping/plot_AIS_grid_3maps.py @@ -1,7 +1,7 @@ #!/usr/bin/env python """ plot_AIS_grid_3maps.py -Written by Tyler Sutterley (05/2023) +Written by Tyler Sutterley (08/2026) Creates 3 GMT-like plots of the Antarctic Ice Sheet on a polar stereographic south (3031) projection @@ -31,6 +31,7 @@ https://pypi.python.org/pypi/GDAL/ UPDATE HISTORY: + Updated 08/2026: use upstream file logger for verbose output Updated 05/2023: use pathlib to define and operate on paths added option to set the input variable names or column order Updated 03/2023: switch from parameter files to argparse arguments @@ -172,21 +173,26 @@ # PURPOSE: keep track of threads def info(args): - logging.info(pathlib.Path(sys.argv[0]).name) - logging.info(args) - logging.info(f'module name: {__name__}') + # get logger + logger = logging.getLogger(__name__) + logger.info(pathlib.Path(sys.argv[0]).name) + logger.info(args) + logger.info(f'module name: {__name__}') if hasattr(os, 'getppid'): - logging.info(f'parent process: {os.getppid():d}') - logging.info(f'process id: {os.getpid():d}') + logger.info(f'parent process: {os.getppid():d}') + logger.info(f'process id: {os.getpid():d}') # PURPOSE: plot Rignot 2012 drainage basin polylines def plot_rignot_basins(ax, base_dir): + # get logger + logger = logging.getLogger(__name__) region_directory = base_dir.joinpath(*region_dir) # for each region for reg in region_title: # read the regional polylines region_file = region_directory.joinpath(region_filename.format(reg)) + logger.debug(str(region_file)) region_ll = np.loadtxt(region_file, dtype=region_dtype) # converting region lat/lon into plot coordinates points = projection.transform_points( @@ -197,9 +203,11 @@ def plot_rignot_basins(ax, base_dir): # PURPOSE: plot Antarctic drainage basins from IMBIE2 (Mouginot) def plot_IMBIE2_basins(ax, base_dir): + # get logger + logger = logging.getLogger(__name__) # read drainage basin polylines from shapefile (using splat operator) basin_shapefile = base_dir.joinpath(*IMBIE_basin_file) - logging.debug(str(basin_shapefile)) + logger.debug(str(basin_shapefile)) shape_input = shapefile.Reader(str(basin_shapefile)) shape_entities = shape_input.shapes() shape_attributes = shape_input.records() @@ -221,10 +229,12 @@ def plot_IMBIE2_basins(ax, base_dir): # PURPOSE: plot Antarctic drainage sub-basins from IMBIE-2 (Mouginot) def plot_IMBIE2_subbasins(ax, base_dir): + # get logger + logger = logging.getLogger(__name__) # read drainage basin polylines from shapefile (using splat operator) IMBIE_subbasin_file = ['Basins_20Oct2016_v1.7', 'Basins_v1.7.shp'] basin_shapefile = base_dir.joinpath('masks', *IMBIE_subbasin_file) - logging.debug(str(basin_shapefile)) + logger.debug(str(basin_shapefile)) shape_input = shapefile.Reader(str(basin_shapefile)) shape_entities = shape_input.shapes() shape_attributes = shape_input.records() @@ -244,8 +254,10 @@ def plot_IMBIE2_subbasins(ax, base_dir): # PURPOSE: plot Antarctic grounded ice delineation def plot_grounded_ice(ax, base_dir, START=1): + # get logger + logger = logging.getLogger(__name__) grounded_ice_shapefile = base_dir.joinpath(*coast_file) - logging.debug(str(grounded_ice_shapefile)) + logger.debug(str(grounded_ice_shapefile)) shape_input = shapefile.Reader(str(grounded_ice_shapefile)) shape_entities = shape_input.shapes() shape_attributes = shape_input.records() @@ -258,9 +270,11 @@ def plot_grounded_ice(ax, base_dir, START=1): # PURPOSE: plot MODIS mosaic of Antarctica as background image def plot_image_mosaic(ax, base_dir, MASKED=True): + # get logger + logger = logging.getLogger(__name__) # read MODIS mosaic of Antarctica image_geotiff_file = base_dir.joinpath(*image_file) - logging.debug(str(image_geotiff_file)) + logger.debug(str(image_geotiff_file)) ds = osgeo.gdal.Open(str(image_geotiff_file)) # get dimensions xsize = ds.RasterXSize @@ -312,11 +326,31 @@ def add_plot_scale(ax, X, Y, dx, dy, masked, fc1='w', fc2='k'): Y - 2.5 * dy, Y + 3.2 * dy, ] - ax.fill([x1, x2, x2, x1, x1], [y1, y1, y2, y2, y1], fc1, zorder=4) + ax.fill( + [x1, x2, x2, x1, x1], + [y1, y1, y2, y2, y1], + fc1, + zorder=4, + ) for i, c in enumerate([fc1, fc2, fc1, fc2]): - x1, x2, y1, y2 = [X + 0.25 * i * dx, X + 0.25 * (i + 1) * dx, Y, Y + dy] - ax.fill([x1, x2, x2, x1, x1], [y1, y1, y2, y2, y1], c, zorder=5) - ax.plot([X, X + dx, X + dx, X, X], [Y, Y, Y + dy, Y + dy, Y], fc2, zorder=6) + x1, x2, y1, y2 = [ + X + 0.25 * i * dx, + X + 0.25 * (i + 1) * dx, + Y, + Y + dy, + ] + ax.fill( + [x1, x2, x2, x1, x1], + [y1, y1, y2, y2, y1], + c, + zorder=5, + ) + ax.plot( + [X, X + dx, X + dx, X, X], + [Y, Y, Y + dy, Y + dy, Y], + fc2, + zorder=6, + ) for i in range(3): ax.plot( [X + 0.5 * i * dx, X + 0.5 * i * dx], @@ -382,6 +416,8 @@ def plot_grid( FIGURE_DPI=None, MODE=0o775, ): + # get logger + logger = logging.getLogger(__name__) # extend list if a single format was entered for all files if len(DATAFORM) < len(FILENAMES): DATAFORM = DATAFORM * len(FILENAMES) @@ -726,7 +762,7 @@ def plot_grid( # create output directory if non-existent FIGURE_FILE.parent.mkdir(mode=MODE, parents=True, exist_ok=True) # save to file - logging.info(str(FIGURE_FILE)) + logger.info(str(FIGURE_FILE)) plt.savefig( FIGURE_FILE, metadata={'Title': pathlib.Path(sys.argv[0]).name}, @@ -982,7 +1018,9 @@ def main(): # create logger loglevels = [logging.CRITICAL, logging.INFO, logging.DEBUG] - logging.basicConfig(level=loglevels[args.verbose]) + logger = gravtk.utilities.build_logger( + __name__, level=loglevels[args.verbose] + ) # try to run the analysis with listed parameters try: @@ -1025,8 +1063,8 @@ def main(): # if there has been an error exception # print the type, value, and stack trace of the # current exception being handled - logging.critical(f'process id {os.getpid():d} failed') - logging.error(traceback.format_exc()) + logger.critical(f'process id {os.getpid():d} failed') + logger.error(traceback.format_exc()) # run main program diff --git a/gravity_toolkit/mapping/plot_AIS_grid_4maps.py b/gravity_toolkit/mapping/plot_AIS_grid_4maps.py index ea9cbff..33fcda0 100644 --- a/gravity_toolkit/mapping/plot_AIS_grid_4maps.py +++ b/gravity_toolkit/mapping/plot_AIS_grid_4maps.py @@ -1,7 +1,7 @@ #!/usr/bin/env python """ plot_AIS_grid_4maps.py -Written by Tyler Sutterley (05/2023) +Written by Tyler Sutterley (08/2026) Creates 3 GMT-like plots of the Antarctic Ice Sheet on a polar stereographic south (3031) projection @@ -31,6 +31,7 @@ https://pypi.python.org/pypi/GDAL/ UPDATE HISTORY: + Updated 08/2026: use upstream file logger for verbose output Updated 05/2023: use pathlib to define and operate on paths added option to set the input variable names or column order Updated 03/2023: switch from parameter files to argparse arguments @@ -173,21 +174,26 @@ # PURPOSE: keep track of threads def info(args): - logging.info(pathlib.Path(sys.argv[0]).name) - logging.info(args) - logging.info(f'module name: {__name__}') + # get logger + logger = logging.getLogger(__name__) + logger.info(pathlib.Path(sys.argv[0]).name) + logger.info(args) + logger.info(f'module name: {__name__}') if hasattr(os, 'getppid'): - logging.info(f'parent process: {os.getppid():d}') - logging.info(f'process id: {os.getpid():d}') + logger.info(f'parent process: {os.getppid():d}') + logger.info(f'process id: {os.getpid():d}') # PURPOSE: plot Rignot 2012 drainage basin polylines def plot_rignot_basins(ax, base_dir): + # get logger + logger = logging.getLogger(__name__) region_directory = base_dir.joinpath(*region_dir) # for each region for reg in region_title: # read the regional polylines region_file = region_directory.joinpath(region_filename.format(reg)) + logger.debug(str(region_file)) region_ll = np.loadtxt(region_file, dtype=region_dtype) # converting region lat/lon into plot coordinates points = projection.transform_points( @@ -198,9 +204,11 @@ def plot_rignot_basins(ax, base_dir): # PURPOSE: plot Antarctic drainage basins from IMBIE2 (Mouginot) def plot_IMBIE2_basins(ax, base_dir): + # get logger + logger = logging.getLogger(__name__) # read drainage basin polylines from shapefile (using splat operator) basin_shapefile = base_dir.joinpath(*IMBIE_basin_file) - logging.debug(str(basin_shapefile)) + logger.debug(str(basin_shapefile)) shape_input = shapefile.Reader(str(basin_shapefile)) shape_entities = shape_input.shapes() shape_attributes = shape_input.records() @@ -222,10 +230,12 @@ def plot_IMBIE2_basins(ax, base_dir): # PURPOSE: plot Antarctic drainage sub-basins from IMBIE-2 (Mouginot) def plot_IMBIE2_subbasins(ax, base_dir): + # get logger + logger = logging.getLogger(__name__) # read drainage basin polylines from shapefile (using splat operator) IMBIE_subbasin_file = ['Basins_20Oct2016_v1.7', 'Basins_v1.7.shp'] basin_shapefile = base_dir.joinpath('masks', *IMBIE_subbasin_file) - logging.debug(str(basin_shapefile)) + logger.debug(str(basin_shapefile)) shape_input = shapefile.Reader(str(basin_shapefile)) shape_entities = shape_input.shapes() shape_attributes = shape_input.records() @@ -245,8 +255,10 @@ def plot_IMBIE2_subbasins(ax, base_dir): # PURPOSE: plot Antarctic grounded ice delineation def plot_grounded_ice(ax, base_dir, START=1): + # get logger + logger = logging.getLogger(__name__) grounded_ice_shapefile = base_dir.joinpath(*coast_file) - logging.debug(str(grounded_ice_shapefile)) + logger.debug(str(grounded_ice_shapefile)) shape_input = shapefile.Reader(str(grounded_ice_shapefile)) shape_entities = shape_input.shapes() shape_attributes = shape_input.records() @@ -259,9 +271,11 @@ def plot_grounded_ice(ax, base_dir, START=1): # PURPOSE: plot MODIS mosaic of Antarctica as background image def plot_image_mosaic(ax, base_dir, MASKED=True): + # get logger + logger = logging.getLogger(__name__) # read MODIS mosaic of Antarctica image_geotiff_file = base_dir.joinpath(*image_file) - logging.debug(str(image_geotiff_file)) + logger.debug(str(image_geotiff_file)) ds = osgeo.gdal.Open(str(image_geotiff_file)) # get dimensions xsize = ds.RasterXSize @@ -313,11 +327,31 @@ def add_plot_scale(ax, X, Y, dx, dy, masked, fc1='w', fc2='k'): Y - 2.5 * dy, Y + 3.2 * dy, ] - ax.fill([x1, x2, x2, x1, x1], [y1, y1, y2, y2, y1], fc1, zorder=4) + ax.fill( + [x1, x2, x2, x1, x1], + [y1, y1, y2, y2, y1], + fc1, + zorder=4, + ) for i, c in enumerate([fc1, fc2, fc1, fc2]): - x1, x2, y1, y2 = [X + 0.25 * i * dx, X + 0.25 * (i + 1) * dx, Y, Y + dy] - ax.fill([x1, x2, x2, x1, x1], [y1, y1, y2, y2, y1], c, zorder=5) - ax.plot([X, X + dx, X + dx, X, X], [Y, Y, Y + dy, Y + dy, Y], fc2, zorder=6) + x1, x2, y1, y2 = [ + X + 0.25 * i * dx, + X + 0.25 * (i + 1) * dx, + Y, + Y + dy, + ] + ax.fill( + [x1, x2, x2, x1, x1], + [y1, y1, y2, y2, y1], + c, + zorder=5, + ) + ax.plot( + [X, X + dx, X + dx, X, X], + [Y, Y, Y + dy, Y + dy, Y], + fc2, + zorder=6, + ) for i in range(3): ax.plot( [X + 0.5 * i * dx, X + 0.5 * i * dx], @@ -383,6 +417,8 @@ def plot_grid( FIGURE_DPI=None, MODE=0o775, ): + # get logger + logger = logging.getLogger(__name__) # extend list if a single format was entered for all files if len(DATAFORM) < len(FILENAMES): DATAFORM = DATAFORM * len(FILENAMES) @@ -737,7 +773,7 @@ def plot_grid( # create output directory if non-existent FIGURE_FILE.parent.mkdir(mode=MODE, parents=True, exist_ok=True) # save to file - logging.info(str(FIGURE_FILE)) + logger.info(str(FIGURE_FILE)) plt.savefig( FIGURE_FILE, metadata={'Title': pathlib.Path(sys.argv[0]).name}, @@ -993,7 +1029,9 @@ def main(): # create logger loglevels = [logging.CRITICAL, logging.INFO, logging.DEBUG] - logging.basicConfig(level=loglevels[args.verbose]) + logger = gravtk.utilities.build_logger( + __name__, level=loglevels[args.verbose] + ) # try to run the analysis with listed parameters try: @@ -1036,8 +1074,8 @@ def main(): # if there has been an error exception # print the type, value, and stack trace of the # current exception being handled - logging.critical(f'process id {os.getpid():d} failed') - logging.error(traceback.format_exc()) + logger.critical(f'process id {os.getpid():d} failed') + logger.error(traceback.format_exc()) # run main program diff --git a/gravity_toolkit/mapping/plot_AIS_grid_maps.py b/gravity_toolkit/mapping/plot_AIS_grid_maps.py index 3d8c6c2..314a8c7 100644 --- a/gravity_toolkit/mapping/plot_AIS_grid_maps.py +++ b/gravity_toolkit/mapping/plot_AIS_grid_maps.py @@ -1,7 +1,7 @@ #!/usr/bin/env python """ plot_AIS_grid_maps.py -Written by Tyler Sutterley (05/2023) +Written by Tyler Sutterley (08/2026) Creates GMT-like plots for the Antarctic Ice Sheet on a polar stereographic south (3031) projection @@ -31,6 +31,7 @@ https://pypi.python.org/pypi/GDAL/ UPDATE HISTORY: + Updated 08/2026: use upstream file logger for verbose output Updated 05/2023: use pathlib to define and operate on paths added option to set the input variable names or column order Updated 03/2023: switch from parameter files to argparse arguments @@ -186,21 +187,26 @@ # PURPOSE: keep track of threads def info(args): - logging.info(pathlib.Path(sys.argv[0]).name) - logging.info(args) - logging.info(f'module name: {__name__}') + # get logger + logger = logging.getLogger(__name__) + logger.info(pathlib.Path(sys.argv[0]).name) + logger.info(args) + logger.info(f'module name: {__name__}') if hasattr(os, 'getppid'): - logging.info(f'parent process: {os.getppid():d}') - logging.info(f'process id: {os.getpid():d}') + logger.info(f'parent process: {os.getppid():d}') + logger.info(f'process id: {os.getpid():d}') # PURPOSE: plot Rignot 2012 drainage basin polylines def plot_rignot_basins(ax, base_dir): + # get logger + logger = logging.getLogger(__name__) region_directory = base_dir.joinpath(*region_dir) # for each region for reg in region_title: # read the regional polylines region_file = region_directory.joinpath(region_filename.format(reg)) + logger.debug(str(region_file)) region_ll = np.loadtxt(region_file, dtype=region_dtype) # converting region lat/lon into plot coordinates points = projection.transform_points( @@ -211,10 +217,11 @@ def plot_rignot_basins(ax, base_dir): # PURPOSE: plot Antarctic drainage basins from IMBIE2 (Mouginot) def plot_IMBIE2_basins(ax, base_dir): + # get logger + logger = logging.getLogger(__name__) # read drainage basin polylines from shapefile (using splat operator) basin_shapefile = base_dir.joinpath(*IMBIE_basin_file) - logging.debug(str(basin_shapefile)) - logging.debug(str(basin_shapefile)) + logger.debug(str(basin_shapefile)) shape_input = shapefile.Reader(str(basin_shapefile)) shape_entities = shape_input.shapes() shape_attributes = shape_input.records() @@ -236,11 +243,12 @@ def plot_IMBIE2_basins(ax, base_dir): # PURPOSE: plot Antarctic drainage sub-basins from IMBIE-2 (Mouginot) def plot_IMBIE2_subbasins(ax, base_dir): + # get logger + logger = logging.getLogger(__name__) # read drainage basin polylines from shapefile (using splat operator) IMBIE_subbasin_file = ['Basins_20Oct2016_v1.7', 'Basins_v1.7.shp'] basin_shapefile = base_dir.joinpath('masks', *IMBIE_subbasin_file) - logging.debug(str(basin_shapefile)) - logging.debug(str(basin_shapefile)) + logger.debug(str(basin_shapefile)) shape_input = shapefile.Reader(str(basin_shapefile)) shape_entities = shape_input.shapes() shape_attributes = shape_input.records() @@ -260,8 +268,10 @@ def plot_IMBIE2_subbasins(ax, base_dir): # PURPOSE: plot Antarctic grounded ice delineation def plot_grounded_ice(ax, base_dir, START=1): + # get logger + logger = logging.getLogger(__name__) grounded_ice_shapefile = base_dir.joinpath(*coast_file) - logging.debug(str(grounded_ice_shapefile)) + logger.debug(str(grounded_ice_shapefile)) shape_input = shapefile.Reader(str(grounded_ice_shapefile)) shape_entities = shape_input.shapes() shape_attributes = shape_input.records() @@ -274,9 +284,11 @@ def plot_grounded_ice(ax, base_dir, START=1): # PURPOSE: plot MODIS mosaic of Antarctica as background image def plot_image_mosaic(ax, base_dir, MASKED=True): + # get logger + logger = logging.getLogger(__name__) # read MODIS mosaic of Antarctica image_geotiff_file = base_dir.joinpath(*image_file) - logging.debug(str(image_geotiff_file)) + logger.debug(str(image_geotiff_file)) ds = osgeo.gdal.Open(str(image_geotiff_file)) # get dimensions xsize = ds.RasterXSize @@ -328,11 +340,31 @@ def add_plot_scale(ax, X, Y, dx, dy, masked, fc1='w', fc2='k'): Y - 2.5 * dy, Y + 3.2 * dy, ] - ax.fill([x1, x2, x2, x1, x1], [y1, y1, y2, y2, y1], fc1, zorder=4) + ax.fill( + [x1, x2, x2, x1, x1], + [y1, y1, y2, y2, y1], + fc1, + zorder=4, + ) for i, c in enumerate([fc1, fc2, fc1, fc2]): - x1, x2, y1, y2 = [X + 0.25 * i * dx, X + 0.25 * (i + 1) * dx, Y, Y + dy] - ax.fill([x1, x2, x2, x1, x1], [y1, y1, y2, y2, y1], c, zorder=5) - ax.plot([X, X + dx, X + dx, X, X], [Y, Y, Y + dy, Y + dy, Y], fc2, zorder=6) + x1, x2, y1, y2 = [ + X + 0.25 * i * dx, + X + 0.25 * (i + 1) * dx, + Y, + Y + dy, + ] + ax.fill( + [x1, x2, x2, x1, x1], + [y1, y1, y2, y2, y1], + c, + zorder=5, + ) + ax.plot( + [X, X + dx, X + dx, X, X], + [Y, Y, Y + dy, Y + dy, Y], + fc2, + zorder=6, + ) for i in range(3): ax.plot( [X + 0.5 * i * dx, X + 0.5 * i * dx], @@ -397,6 +429,8 @@ def plot_grid( FIGURE_DPI=None, MODE=0o775, ): + # get logger + logger = logging.getLogger(__name__) # read CPT or use color map if CPT_FILE is not None: # cpt file @@ -740,7 +774,7 @@ def plot_grid( # create output directory if non-existent FIGURE_FILE.parent.mkdir(mode=MODE, parents=True, exist_ok=True) # save to file - logging.info(str(FIGURE_FILE)) + logger.info(str(FIGURE_FILE)) plt.savefig( FIGURE_FILE, metadata={'Title': pathlib.Path(sys.argv[0]).name}, @@ -993,7 +1027,9 @@ def main(): # create logger loglevels = [logging.CRITICAL, logging.INFO, logging.DEBUG] - logging.basicConfig(level=loglevels[args.verbose]) + logger = gravtk.utilities.build_logger( + __name__, level=loglevels[args.verbose] + ) # try to run the analysis with listed parameters try: @@ -1036,8 +1072,8 @@ def main(): # if there has been an error exception # print the type, value, and stack trace of the # current exception being handled - logging.critical(f'process id {os.getpid():d} failed') - logging.error(traceback.format_exc()) + logger.critical(f'process id {os.getpid():d} failed') + logger.error(traceback.format_exc()) # run main program diff --git a/gravity_toolkit/mapping/plot_AIS_grid_movie.py b/gravity_toolkit/mapping/plot_AIS_grid_movie.py index b11409b..a93e556 100644 --- a/gravity_toolkit/mapping/plot_AIS_grid_movie.py +++ b/gravity_toolkit/mapping/plot_AIS_grid_movie.py @@ -1,7 +1,7 @@ #!/usr/bin/env python """ plot_AIS_grid_movie.py -Written by Tyler Sutterley (05/2023) +Written by Tyler Sutterley (08/2026) Creates GMT-like animations for the Antarctic Ice Sheet on a polar stereographic south (3031) projection @@ -31,6 +31,7 @@ https://pypi.python.org/pypi/GDAL/ UPDATE HISTORY: + Updated 08/2026: use upstream file logger for verbose output Updated 05/2023: use pathlib to define and operate on paths Updated 03/2023: switch from parameter files to argparse arguments updated inputs to spatial from_file function @@ -192,21 +193,26 @@ # PURPOSE: keep track of threads def info(args): - logging.info(pathlib.Path(sys.argv[0]).name) - logging.info(args) - logging.info(f'module name: {__name__}') + # get logger + logger = logging.getLogger(__name__) + logger.info(pathlib.Path(sys.argv[0]).name) + logger.info(args) + logger.info(f'module name: {__name__}') if hasattr(os, 'getppid'): - logging.info(f'parent process: {os.getppid():d}') - logging.info(f'process id: {os.getpid():d}') + logger.info(f'parent process: {os.getppid():d}') + logger.info(f'process id: {os.getpid():d}') # PURPOSE: plot Rignot 2012 drainage basin polylines def plot_rignot_basins(ax, base_dir): + # get logger + logger = logging.getLogger(__name__) region_directory = base_dir.joinpath(*region_dir) # for each region for reg in region_title: # read the regional polylines region_file = region_directory.joinpath(region_filename.format(reg)) + logger.debug(str(region_file)) region_ll = np.loadtxt(region_file, dtype=region_dtype) # converting region lat/lon into plot coordinates points = projection.transform_points( @@ -217,9 +223,11 @@ def plot_rignot_basins(ax, base_dir): # PURPOSE: plot Antarctic drainage basins from IMBIE2 (Mouginot) def plot_IMBIE2_basins(ax, base_dir): + # get logger + logger = logging.getLogger(__name__) # read drainage basin polylines from shapefile (using splat operator) basin_shapefile = base_dir.joinpath(*IMBIE_basin_file) - logging.debug(str(basin_shapefile)) + logger.debug(str(basin_shapefile)) shape_input = shapefile.Reader(str(basin_shapefile)) shape_entities = shape_input.shapes() shape_attributes = shape_input.records() @@ -241,10 +249,12 @@ def plot_IMBIE2_basins(ax, base_dir): # PURPOSE: plot Antarctic drainage sub-basins from IMBIE-2 (Mouginot) def plot_IMBIE2_subbasins(ax, base_dir): + # get logger + logger = logging.getLogger(__name__) # read drainage basin polylines from shapefile (using splat operator) IMBIE_subbasin_file = ['Basins_20Oct2016_v1.7', 'Basins_v1.7.shp'] basin_shapefile = base_dir.joinpath('masks', *IMBIE_subbasin_file) - logging.debug(str(basin_shapefile)) + logger.debug(str(basin_shapefile)) shape_input = shapefile.Reader(str(basin_shapefile)) shape_entities = shape_input.shapes() shape_attributes = shape_input.records() @@ -264,8 +274,10 @@ def plot_IMBIE2_subbasins(ax, base_dir): # PURPOSE: plot Antarctic grounded ice delineation def plot_grounded_ice(ax, base_dir, START=1): + # get logger + logger = logging.getLogger(__name__) grounded_ice_shapefile = base_dir.joinpath(*coast_file) - logging.debug(str(grounded_ice_shapefile)) + logger.debug(str(grounded_ice_shapefile)) shape_input = shapefile.Reader(str(grounded_ice_shapefile)) shape_entities = shape_input.shapes() shape_attributes = shape_input.records() @@ -278,9 +290,11 @@ def plot_grounded_ice(ax, base_dir, START=1): # PURPOSE: plot MODIS mosaic of Antarctica as background image def plot_image_mosaic(ax, base_dir, MASKED=True): + # get logger + logger = logging.getLogger(__name__) # read MODIS mosaic of Antarctica image_geotiff_file = base_dir.joinpath(*image_file) - logging.debug(str(image_geotiff_file)) + logger.debug(str(image_geotiff_file)) ds = osgeo.gdal.Open(str(image_geotiff_file)) # get dimensions xsize = ds.RasterXSize @@ -332,11 +346,31 @@ def add_plot_scale(ax, X, Y, dx, dy, masked, fc1='w', fc2='k'): Y - 2.5 * dy, Y + 3.2 * dy, ] - ax.fill([x1, x2, x2, x1, x1], [y1, y1, y2, y2, y1], fc1, zorder=4) + ax.fill( + [x1, x2, x2, x1, x1], + [y1, y1, y2, y2, y1], + fc1, + zorder=4, + ) for i, c in enumerate([fc1, fc2, fc1, fc2]): - x1, x2, y1, y2 = [X + 0.25 * i * dx, X + 0.25 * (i + 1) * dx, Y, Y + dy] - ax.fill([x1, x2, x2, x1, x1], [y1, y1, y2, y2, y1], c, zorder=5) - ax.plot([X, X + dx, X + dx, X, X], [Y, Y, Y + dy, Y + dy, Y], fc2, zorder=6) + x1, x2, y1, y2 = [ + X + 0.25 * i * dx, + X + 0.25 * (i + 1) * dx, + Y, + Y + dy, + ] + ax.fill( + [x1, x2, x2, x1, x1], + [y1, y1, y2, y2, y1], + c, + zorder=5, + ) + ax.plot( + [X, X + dx, X + dx, X, X], + [Y, Y, Y + dy, Y + dy, Y], + fc2, + zorder=6, + ) for i in range(3): ax.plot( [X + 0.5 * i * dx, X + 0.5 * i * dx], @@ -400,6 +434,8 @@ def animate_grid( FIGURE_DPI=None, MODE=0o775, ): + # get logger + logger = logging.getLogger(__name__) # read CPT or use color map if CPT_FILE is not None: # cpt file @@ -1041,7 +1077,9 @@ def main(): # create logger loglevels = [logging.CRITICAL, logging.INFO, logging.DEBUG] - logging.basicConfig(level=loglevels[args.verbose]) + logger = gravtk.utilities.build_logger( + __name__, level=loglevels[args.verbose] + ) # try to run the analysis with listed parameters try: @@ -1082,8 +1120,8 @@ def main(): # if there has been an error exception # print the type, value, and stack trace of the # current exception being handled - logging.critical(f'process id {os.getpid():d} failed') - logging.error(traceback.format_exc()) + logger.critical(f'process id {os.getpid():d} failed') + logger.error(traceback.format_exc()) # run main program diff --git a/gravity_toolkit/mapping/plot_AIS_regional_maps.py b/gravity_toolkit/mapping/plot_AIS_regional_maps.py index d44c340..51e4389 100644 --- a/gravity_toolkit/mapping/plot_AIS_regional_maps.py +++ b/gravity_toolkit/mapping/plot_AIS_regional_maps.py @@ -1,7 +1,7 @@ #!/usr/bin/env python """ plot_AIS_regional_maps.py -Written by Tyler Sutterley (05/2023) +Written by Tyler Sutterley (08/2026) Creates GMT-like plots for sub-regions of Antarctica on a polar stereographic south (3031) projection @@ -31,6 +31,7 @@ https://pypi.python.org/pypi/GDAL/ UPDATE HISTORY: + Updated 08/2026: use upstream file logger for verbose output Updated 05/2023: use pathlib to define and operate on paths added option to set the input variable names or column order Updated 03/2023: switch from parameter files to argparse arguments @@ -248,21 +249,26 @@ # PURPOSE: keep track of threads def info(args): - logging.info(pathlib.Path(sys.argv[0]).name) - logging.info(args) - logging.info(f'module name: {__name__}') + # get logger + logger = logging.getLogger(__name__) + logger.info(pathlib.Path(sys.argv[0]).name) + logger.info(args) + logger.info(f'module name: {__name__}') if hasattr(os, 'getppid'): - logging.info(f'parent process: {os.getppid():d}') - logging.info(f'process id: {os.getpid():d}') + logger.info(f'parent process: {os.getppid():d}') + logger.info(f'process id: {os.getpid():d}') # PURPOSE: plot Rignot 2012 drainage basin polylines def plot_rignot_basins(ax, base_dir): + # get logger + logger = logging.getLogger(__name__) region_directory = base_dir.joinpath(*region_dir) # for each region for reg in region_title: # read the regional polylines region_file = region_directory.joinpath(region_filename.format(reg)) + logger.debug(str(region_file)) region_ll = np.loadtxt(region_file, dtype=region_dtype) # converting region lat/lon into plot coordinates points = projection.transform_points( @@ -273,9 +279,11 @@ def plot_rignot_basins(ax, base_dir): # PURPOSE: plot Antarctic drainage basins from IMBIE2 (Mouginot) def plot_IMBIE2_basins(ax, base_dir): + # get logger + logger = logging.getLogger(__name__) # read drainage basin polylines from shapefile (using splat operator) basin_shapefile = base_dir.joinpath(*IMBIE_basin_file) - logging.debug(str(basin_shapefile)) + logger.debug(str(basin_shapefile)) shape_input = shapefile.Reader(str(basin_shapefile)) shape_entities = shape_input.shapes() shape_attributes = shape_input.records() @@ -297,10 +305,12 @@ def plot_IMBIE2_basins(ax, base_dir): # PURPOSE: plot Antarctic drainage sub-basins from IMBIE-2 (Mouginot) def plot_IMBIE2_subbasins(ax, base_dir): + # get logger + logger = logging.getLogger(__name__) # read drainage basin polylines from shapefile (using splat operator) IMBIE_subbasin_file = ['Basins_20Oct2016_v1.7', 'Basins_v1.7.shp'] basin_shapefile = base_dir.joinpath('masks', *IMBIE_subbasin_file) - logging.debug(str(basin_shapefile)) + logger.debug(str(basin_shapefile)) shape_input = shapefile.Reader(str(basin_shapefile)) shape_entities = shape_input.shapes() shape_attributes = shape_input.records() @@ -320,12 +330,14 @@ def plot_IMBIE2_subbasins(ax, base_dir): # PURPOSE: plot Amundsen Sea basins from Mouginot et al. (2014) def plot_amundsen_basins(ax, base_dir): + # get logger + logger = logging.getLogger(__name__) # read Amundsen Sea basin polylines from shapefile basin_shapefile = base_dir.joinpath( 'masks', 'Basins_Admunsen', 'Basins_admunsen_match_coastline_and_IS.shp' ) basin_title = ['pope_smith', 'haynes', 'thwaites', 'pine_island', 'kohler'] - logging.debug(str(basin_shapefile)) + logger.debug(str(basin_shapefile)) shape_input = shapefile.Reader(str(basin_shapefile)) shape_entities = shape_input.shapes() # for each shape entity @@ -337,8 +349,10 @@ def plot_amundsen_basins(ax, base_dir): # PURPOSE: plot Antarctic grounded ice delineation def plot_grounded_ice(ax, base_dir, START=1): + # get logger + logger = logging.getLogger(__name__) grounded_ice_shapefile = base_dir.joinpath(*coast_file) - logging.debug(str(grounded_ice_shapefile)) + logger.debug(str(grounded_ice_shapefile)) shape_input = shapefile.Reader(str(grounded_ice_shapefile)) shape_entities = shape_input.shapes() shape_attributes = shape_input.records() @@ -351,9 +365,11 @@ def plot_grounded_ice(ax, base_dir, START=1): # PURPOSE: plot MODIS mosaic of Antarctica as background image def plot_image_mosaic(ax, base_dir, xlimits, ylimits, MASKED=True): + # get logger + logger = logging.getLogger(__name__) # read MODIS mosaic of Antarctica image_geotiff_file = base_dir.joinpath(*image_file) - logging.debug(str(image_geotiff_file)) + logger.debug(str(image_geotiff_file)) ds = osgeo.gdal.Open(str(image_geotiff_file)) # get geotiff info info_geotiff = ds.GetGeoTransform() @@ -412,11 +428,31 @@ def add_plot_scale(ax, X, Y, dx, dy, masked, fc1='w', fc2='k'): Y - 1.8 * dy, Y + 2.4 * dy, ] - ax.fill([x1, x2, x2, x1, x1], [y1, y1, y2, y2, y1], fc1, zorder=4) + ax.fill( + [x1, x2, x2, x1, x1], + [y1, y1, y2, y2, y1], + fc1, + zorder=4, + ) for i, c in enumerate([fc1, fc2, fc1, fc2]): - x1, x2, y1, y2 = [X + 0.25 * i * dx, X + 0.25 * (i + 1) * dx, Y, Y + dy] - ax.fill([x1, x2, x2, x1, x1], [y1, y1, y2, y2, y1], c, zorder=5) - ax.plot([X, X + dx, X + dx, X, X], [Y, Y, Y + dy, Y + dy, Y], fc2, zorder=6) + x1, x2, y1, y2 = [ + X + 0.25 * i * dx, + X + 0.25 * (i + 1) * dx, + Y, + Y + dy, + ] + ax.fill( + [x1, x2, x2, x1, x1], + [y1, y1, y2, y2, y1], + c, + zorder=5, + ) + ax.plot( + [X, X + dx, X + dx, X, X], + [Y, Y, Y + dy, Y + dy, Y], + fc2, + zorder=6, + ) for i in range(3): ax.plot( [X + 0.5 * i * dx, X + 0.5 * i * dx], @@ -482,6 +518,8 @@ def plot_grid( FIGURE_DPI=None, MODE=0o775, ): + # get logger + logger = logging.getLogger(__name__) # read CPT or use color map if CPT_FILE is not None: # cpt file @@ -830,7 +868,7 @@ def plot_grid( # create output directory if non-existent FIGURE_FILE.parent.mkdir(mode=MODE, parents=True, exist_ok=True) # save to file - logging.info(str(FIGURE_FILE)) + logger.info(str(FIGURE_FILE)) plt.savefig( FIGURE_FILE, metadata={'Title': pathlib.Path(sys.argv[0]).name}, @@ -1093,7 +1131,9 @@ def main(): # create logger loglevels = [logging.CRITICAL, logging.INFO, logging.DEBUG] - logging.basicConfig(level=loglevels[args.verbose]) + logger = gravtk.utilities.build_logger( + __name__, level=loglevels[args.verbose] + ) # try to run the analysis with listed parameters try: @@ -1137,8 +1177,8 @@ def main(): # if there has been an error exception # print the type, value, and stack trace of the # current exception being handled - logging.critical(f'process id {os.getpid():d} failed') - logging.error(traceback.format_exc()) + logger.critical(f'process id {os.getpid():d} failed') + logger.error(traceback.format_exc()) # run main program diff --git a/gravity_toolkit/mapping/plot_AIS_regional_movie.py b/gravity_toolkit/mapping/plot_AIS_regional_movie.py index babbf1a..2847177 100644 --- a/gravity_toolkit/mapping/plot_AIS_regional_movie.py +++ b/gravity_toolkit/mapping/plot_AIS_regional_movie.py @@ -1,7 +1,7 @@ #!/usr/bin/env python """ plot_ASE_grid_movie.py -Written by Tyler Sutterley (05/2023) +Written by Tyler Sutterley (08/2026) Creates GMT-like animations for sub-regions of Antarctica on a polar stereographic south (3031) projection @@ -31,6 +31,7 @@ https://pypi.python.org/pypi/GDAL/ UPDATE HISTORY: + Updated 08/2026: use upstream file logger for verbose output Updated 05/2023: use pathlib to define and operate on paths Updated 03/2023: switch from parameter files to argparse arguments updated inputs to spatial from_file function @@ -253,21 +254,26 @@ # PURPOSE: keep track of threads def info(args): - logging.info(pathlib.Path(sys.argv[0]).name) - logging.info(args) - logging.info(f'module name: {__name__}') + # get logger + logger = logging.getLogger(__name__) + logger.info(pathlib.Path(sys.argv[0]).name) + logger.info(args) + logger.info(f'module name: {__name__}') if hasattr(os, 'getppid'): - logging.info(f'parent process: {os.getppid():d}') - logging.info(f'process id: {os.getpid():d}') + logger.info(f'parent process: {os.getppid():d}') + logger.info(f'process id: {os.getpid():d}') # PURPOSE: plot Rignot 2012 drainage basin polylines def plot_rignot_basins(ax, base_dir): + # get logger + logger = logging.getLogger(__name__) region_directory = base_dir.joinpath(*region_dir) # for each region for reg in region_title: # read the regional polylines region_file = region_directory.joinpath(region_filename.format(reg)) + logger.debug(str(region_file)) region_ll = np.loadtxt(region_file, dtype=region_dtype) # converting region lat/lon into plot coordinates points = projection.transform_points( @@ -278,9 +284,11 @@ def plot_rignot_basins(ax, base_dir): # PURPOSE: plot Antarctic drainage basins from IMBIE2 (Mouginot) def plot_IMBIE2_basins(ax, base_dir): + # get logger + logger = logging.getLogger(__name__) # read drainage basin polylines from shapefile (using splat operator) basin_shapefile = base_dir.joinpath(*IMBIE_basin_file) - logging.debug(str(basin_shapefile)) + logger.debug(str(basin_shapefile)) shape_input = shapefile.Reader(str(basin_shapefile)) shape_entities = shape_input.shapes() shape_attributes = shape_input.records() @@ -302,10 +310,12 @@ def plot_IMBIE2_basins(ax, base_dir): # PURPOSE: plot Antarctic drainage sub-basins from IMBIE-2 (Mouginot) def plot_IMBIE2_subbasins(ax, base_dir): + # get logger + logger = logging.getLogger(__name__) # read drainage basin polylines from shapefile (using splat operator) IMBIE_subbasin_file = ['Basins_20Oct2016_v1.7', 'Basins_v1.7.shp'] basin_shapefile = base_dir.joinpath('masks', *IMBIE_subbasin_file) - logging.debug(str(basin_shapefile)) + logger.debug(str(basin_shapefile)) shape_input = shapefile.Reader(str(basin_shapefile)) shape_entities = shape_input.shapes() shape_attributes = shape_input.records() @@ -325,12 +335,14 @@ def plot_IMBIE2_subbasins(ax, base_dir): # PURPOSE: plot Amundsen Sea basins from Mouginot et al. (2014) def plot_amundsen_basins(ax, base_dir): + # get logger + logger = logging.getLogger(__name__) # read Amundsen Sea basin polylines from shapefile basin_shapefile = base_dir.joinpath( 'masks', 'Basins_Admunsen', 'Basins_admunsen_match_coastline_and_IS.shp' ) basin_title = ['pope_smith', 'haynes', 'thwaites', 'pine_island', 'kohler'] - logging.debug(str(basin_shapefile)) + logger.debug(str(basin_shapefile)) shape_input = shapefile.Reader(str(basin_shapefile)) shape_entities = shape_input.shapes() # for each shape entity @@ -342,8 +354,10 @@ def plot_amundsen_basins(ax, base_dir): # PURPOSE: plot Antarctic grounded ice delineation def plot_grounded_ice(ax, base_dir, START=1): + # get logger + logger = logging.getLogger(__name__) grounded_ice_shapefile = base_dir.joinpath(*coast_file) - logging.debug(str(grounded_ice_shapefile)) + logger.debug(str(grounded_ice_shapefile)) shape_input = shapefile.Reader(str(grounded_ice_shapefile)) shape_entities = shape_input.shapes() shape_attributes = shape_input.records() @@ -356,9 +370,11 @@ def plot_grounded_ice(ax, base_dir, START=1): # PURPOSE: plot MODIS mosaic of Antarctica as background image def plot_image_mosaic(ax, base_dir, xlimits, ylimits, MASKED=True): + # get logger + logger = logging.getLogger(__name__) # read MODIS mosaic of Antarctica image_geotiff_file = base_dir.joinpath(*image_file) - logging.debug(str(image_geotiff_file)) + logger.debug(str(image_geotiff_file)) ds = osgeo.gdal.Open(str(image_geotiff_file)) # get geotiff info info_geotiff = ds.GetGeoTransform() @@ -417,11 +433,31 @@ def add_plot_scale(ax, X, Y, dx, dy, masked, fc1='w', fc2='k'): Y - 1.8 * dy, Y + 2.4 * dy, ] - ax.fill([x1, x2, x2, x1, x1], [y1, y1, y2, y2, y1], fc1, zorder=4) + ax.fill( + [x1, x2, x2, x1, x1], + [y1, y1, y2, y2, y1], + fc1, + zorder=4, + ) for i, c in enumerate([fc1, fc2, fc1, fc2]): - x1, x2, y1, y2 = [X + 0.25 * i * dx, X + 0.25 * (i + 1) * dx, Y, Y + dy] - ax.fill([x1, x2, x2, x1, x1], [y1, y1, y2, y2, y1], c, zorder=5) - ax.plot([X, X + dx, X + dx, X, X], [Y, Y, Y + dy, Y + dy, Y], fc2, zorder=6) + x1, x2, y1, y2 = [ + X + 0.25 * i * dx, + X + 0.25 * (i + 1) * dx, + Y, + Y + dy, + ] + ax.fill( + [x1, x2, x2, x1, x1], + [y1, y1, y2, y2, y1], + c, + zorder=5, + ) + ax.plot( + [X, X + dx, X + dx, X, X], + [Y, Y, Y + dy, Y + dy, Y], + fc2, + zorder=6, + ) for i in range(3): ax.plot( [X + 0.5 * i * dx, X + 0.5 * i * dx], @@ -486,6 +522,8 @@ def animate_grid( FIGURE_DPI=None, MODE=0o775, ): + # get logger + logger = logging.getLogger(__name__) # read CPT or use color map if CPT_FILE is not None: # cpt file @@ -1141,7 +1179,9 @@ def main(): # create logger loglevels = [logging.CRITICAL, logging.INFO, logging.DEBUG] - logging.basicConfig(level=loglevels[args.verbose]) + logger = gravtk.utilities.build_logger( + __name__, level=loglevels[args.verbose] + ) # try to run the analysis with listed parameters try: @@ -1183,8 +1223,8 @@ def main(): # if there has been an error exception # print the type, value, and stack trace of the # current exception being handled - logging.critical(f'process id {os.getpid():d} failed') - logging.error(traceback.format_exc()) + logger.critical(f'process id {os.getpid():d} failed') + logger.error(traceback.format_exc()) # run main program diff --git a/gravity_toolkit/mapping/plot_GrIS_grid_3maps.py b/gravity_toolkit/mapping/plot_GrIS_grid_3maps.py index 1c44d2f..172bbdc 100644 --- a/gravity_toolkit/mapping/plot_GrIS_grid_3maps.py +++ b/gravity_toolkit/mapping/plot_GrIS_grid_3maps.py @@ -1,7 +1,7 @@ #!/usr/bin/env python """ plot_GrIS_grid_3maps.py -Written by Tyler Sutterley (05/2023) +Written by Tyler Sutterley (08/2026) Creates 3 GMT-like plots for the Greenland ice sheet on a NSIDC polar stereographic north (3413) projection @@ -34,6 +34,7 @@ https://pypi.python.org/pypi/GDAL/ UPDATE HISTORY: + Updated 08/2026: use upstream file logger for verbose output Updated 05/2023: use pathlib to define and operate on paths added option to set the input variable names or column order Updated 03/2023: switch from parameter files to argparse arguments @@ -132,21 +133,26 @@ # PURPOSE: keep track of threads def info(args): - logging.info(pathlib.Path(sys.argv[0]).name) - logging.info(args) - logging.info(f'module name: {__name__}') + # get logger + logger = logging.getLogger(__name__) + logger.info(pathlib.Path(sys.argv[0]).name) + logger.info(args) + logger.info(f'module name: {__name__}') if hasattr(os, 'getppid'): - logging.info(f'parent process: {os.getppid():d}') - logging.info(f'process id: {os.getpid():d}') + logger.info(f'parent process: {os.getppid():d}') + logger.info(f'process id: {os.getpid():d}') # PURPOSE: plot Rignot 2012 drainage basin polylines def plot_rignot_basins(ax, base_dir): + # get logger + logger = logging.getLogger(__name__) region_directory = base_dir.joinpath(*region_dir) # for each region for reg in region_title: # read the regional polylines region_file = region_directory.joinpath(region_filename.format(reg)) + logger.debug(str(region_file)) region_ll = np.loadtxt(region_file, dtype=region_dtype) # converting region lat/lon into plot coordinates points = projection.transform_points( @@ -157,9 +163,11 @@ def plot_rignot_basins(ax, base_dir): # PURPOSE: plot Greenland drainage basins from IMBIE2 (Mouginot) def plot_IMBIE2_basins(ax, base_dir): + # get logger + logger = logging.getLogger(__name__) # read drainage basin polylines from shapefile (using splat operator) basin_shapefile = base_dir.joinpath(*IMBIE_basin_file) - logging.debug(str(basin_shapefile)) + logger.debug(str(basin_shapefile)) shape_input = shapefile.Reader(str(basin_shapefile)) shape_entities = shape_input.shapes() shape_attributes = shape_input.records() @@ -185,8 +193,10 @@ def plot_IMBIE2_basins(ax, base_dir): # PURPOSE: plot Greenland grounded ice delineation from GIMP def plot_grounded_ice(ax, base_dir, START=1, END=300, LINEWIDTH=0.6): + # get logger + logger = logging.getLogger(__name__) coast_shapefile = base_dir.joinpath(*coast_file) - logging.debug(str(coast_shapefile)) + logger.debug(str(coast_shapefile)) shape_input = shapefile.Reader(str(coast_shapefile)) shape_entities = shape_input.shapes() shape_attributes = shape_input.records() @@ -205,6 +215,9 @@ def plot_grounded_ice(ax, base_dir, START=1, END=300, LINEWIDTH=0.6): # PURPOSE: plot glaciated regions from Randolph Glacier Inventory def plot_glacier_inventory(ax, base_dir, START=0, END=30, LINEWIDTH=0.6): + # get logger + logger = logging.getLogger(__name__) + # RGI shapefiles for Arctic glaciers within plot region RGI_files = [] RGI_files.append('03_rgi60_ArcticCanadaNorth') RGI_files.append('04_rgi60_ArcticCanadaSouth') @@ -212,7 +225,7 @@ def plot_glacier_inventory(ax, base_dir, START=0, END=30, LINEWIDTH=0.6): RGI_files.append('07_rgi60_Svalbard') for f in RGI_files: RGI_shapefile = base_dir.joinpath('RGI', f, f'{f}_plot.shp') - logging.debug(str(RGI_shapefile)) + logger.debug(str(RGI_shapefile)) shape_input = shapefile.Reader(str(RGI_shapefile)) shape_entities = shape_input.shapes() shape_attributes = shape_input.records() @@ -231,6 +244,8 @@ def plot_glacier_inventory(ax, base_dir, START=0, END=30, LINEWIDTH=0.6): # PURPOSE plot coastlines and islands (GSHHS with G250 Greenland) def plot_coastline(ax, base_dir): + # get logger + logger = logging.getLogger(__name__) # read the coastline shape file coastline_dir = base_dir.joinpath('masks', 'G250') coastline_shape_files = [] @@ -238,7 +253,7 @@ def plot_coastline(ax, base_dir): coastline_shape_files.append('greenland_coastline_islands.shp') for fi, S in zip(coastline_shape_files, [1000, 200]): coast_shapefile = coastline_dir.joinpath(fi) - logging.debug(str(coast_shapefile)) + logger.debug(str(coast_shapefile)) shape_input = shapefile.Reader(str(coast_shapefile)) shape_entities = shape_input.shapes() # for each entity within the shapefile @@ -250,9 +265,11 @@ def plot_coastline(ax, base_dir): # plot the MODIS Mosaic of Greenland as a background image def plot_image_mosaic(ax, base_dir, MASKED=True): + # get logger + logger = logging.getLogger(__name__) # read MODIS mosaic of Greenland image_geotiff_file = base_dir.joinpath(*image_file) - logging.debug(str(image_geotiff_file)) + logger.debug(str(image_geotiff_file)) ds = osgeo.gdal.Open(str(image_geotiff_file)) # get dimensions xsize = ds.RasterXSize @@ -302,11 +319,31 @@ def add_plot_scale(ax, X, Y, dx, dy, masked, fc1='w', fc2='k'): Y - 2.5 * dy, Y + 3.2 * dy, ] - ax.fill([x1, x2, x2, x1, x1], [y1, y1, y2, y2, y1], fc1, zorder=4) + ax.fill( + [x1, x2, x2, x1, x1], + [y1, y1, y2, y2, y1], + fc1, + zorder=4, + ) for i, c in enumerate([fc1, fc2, fc1, fc2]): - x1, x2, y1, y2 = [X + 0.25 * i * dx, X + 0.25 * (i + 1) * dx, Y, Y + dy] - ax.fill([x1, x2, x2, x1, x1], [y1, y1, y2, y2, y1], c, zorder=5) - ax.plot([X, X + dx, X + dx, X, X], [Y, Y, Y + dy, Y + dy, Y], fc2, zorder=6) + x1, x2, y1, y2 = [ + X + 0.25 * i * dx, + X + 0.25 * (i + 1) * dx, + Y, + Y + dy, + ] + ax.fill( + [x1, x2, x2, x1, x1], + [y1, y1, y2, y2, y1], + c, + zorder=5, + ) + ax.plot( + [X, X + dx, X + dx, X, X], + [Y, Y, Y + dy, Y + dy, Y], + fc2, + zorder=6, + ) for i in range(3): ax.plot( [X + 0.5 * i * dx, X + 0.5 * i * dx], @@ -373,6 +410,8 @@ def plot_grid( FIGURE_DPI=None, MODE=0o775, ): + # get logger + logger = logging.getLogger(__name__) # extend list if a single format was entered for all files if len(DATAFORM) < len(FILENAMES): DATAFORM = DATAFORM * len(FILENAMES) @@ -721,7 +760,7 @@ def plot_grid( # create output directory if non-existent FIGURE_FILE.parent.mkdir(mode=MODE, parents=True, exist_ok=True) # save to file - logging.info(str(FIGURE_FILE)) + logger.info(str(FIGURE_FILE)) plt.savefig( FIGURE_FILE, metadata={'Title': pathlib.Path(sys.argv[0]).name}, @@ -983,7 +1022,9 @@ def main(): # create logger loglevels = [logging.CRITICAL, logging.INFO, logging.DEBUG] - logging.basicConfig(level=loglevels[args.verbose]) + logger = gravtk.utilities.build_logger( + __name__, level=loglevels[args.verbose] + ) # try to run the analysis with listed parameters try: @@ -1027,8 +1068,8 @@ def main(): # if there has been an error exception # print the type, value, and stack trace of the # current exception being handled - logging.critical(f'process id {os.getpid():d} failed') - logging.error(traceback.format_exc()) + logger.critical(f'process id {os.getpid():d} failed') + logger.error(traceback.format_exc()) # run main program diff --git a/gravity_toolkit/mapping/plot_GrIS_grid_5maps.py b/gravity_toolkit/mapping/plot_GrIS_grid_5maps.py index 4b2a7cf..b03400a 100644 --- a/gravity_toolkit/mapping/plot_GrIS_grid_5maps.py +++ b/gravity_toolkit/mapping/plot_GrIS_grid_5maps.py @@ -34,6 +34,7 @@ https://pypi.python.org/pypi/GDAL/ UPDATE HISTORY: + Updated 08/2026: use upstream file logger for verbose output Written 10/2023 """ @@ -119,21 +120,26 @@ # PURPOSE: keep track of threads def info(args): - logging.info(pathlib.Path(sys.argv[0]).name) - logging.info(args) - logging.info(f'module name: {__name__}') + # get logger + logger = logging.getLogger(__name__) + logger.info(pathlib.Path(sys.argv[0]).name) + logger.info(args) + logger.info(f'module name: {__name__}') if hasattr(os, 'getppid'): - logging.info(f'parent process: {os.getppid():d}') - logging.info(f'process id: {os.getpid():d}') + logger.info(f'parent process: {os.getppid():d}') + logger.info(f'process id: {os.getpid():d}') # PURPOSE: plot Rignot 2012 drainage basin polylines def plot_rignot_basins(ax, base_dir): + # get logger + logger = logging.getLogger(__name__) region_directory = base_dir.joinpath(*region_dir) # for each region for reg in region_title: # read the regional polylines region_file = region_directory.joinpath(region_filename.format(reg)) + logger.debug(str(region_file)) region_ll = np.loadtxt(region_file, dtype=region_dtype) # converting region lat/lon into plot coordinates points = projection.transform_points( @@ -144,9 +150,11 @@ def plot_rignot_basins(ax, base_dir): # PURPOSE: plot Greenland drainage basins from IMBIE2 (Mouginot) def plot_IMBIE2_basins(ax, base_dir): + # get logger + logger = logging.getLogger(__name__) # read drainage basin polylines from shapefile (using splat operator) basin_shapefile = base_dir.joinpath(*IMBIE_basin_file) - logging.debug(str(basin_shapefile)) + logger.debug(str(basin_shapefile)) shape_input = shapefile.Reader(str(basin_shapefile)) shape_entities = shape_input.shapes() shape_attributes = shape_input.records() @@ -172,8 +180,10 @@ def plot_IMBIE2_basins(ax, base_dir): # PURPOSE: plot Greenland grounded ice delineation from GIMP def plot_grounded_ice(ax, base_dir, START=1, END=300, LINEWIDTH=0.6): + # get logger + logger = logging.getLogger(__name__) coast_shapefile = base_dir.joinpath(*coast_file) - logging.debug(str(coast_shapefile)) + logger.debug(str(coast_shapefile)) shape_input = shapefile.Reader(str(coast_shapefile)) shape_entities = shape_input.shapes() shape_attributes = shape_input.records() @@ -192,6 +202,9 @@ def plot_grounded_ice(ax, base_dir, START=1, END=300, LINEWIDTH=0.6): # PURPOSE: plot glaciated regions from Randolph Glacier Inventory def plot_glacier_inventory(ax, base_dir, START=0, END=30, LINEWIDTH=0.6): + # get logger + logger = logging.getLogger(__name__) + # RGI shapefiles for Arctic glaciers within plot region RGI_files = [] RGI_files.append('03_rgi60_ArcticCanadaNorth') RGI_files.append('04_rgi60_ArcticCanadaSouth') @@ -199,7 +212,7 @@ def plot_glacier_inventory(ax, base_dir, START=0, END=30, LINEWIDTH=0.6): RGI_files.append('07_rgi60_Svalbard') for f in RGI_files: RGI_shapefile = base_dir.joinpath('RGI', f, f'{f}_plot.shp') - logging.debug(str(RGI_shapefile)) + logger.debug(str(RGI_shapefile)) shape_input = shapefile.Reader(str(RGI_shapefile)) shape_entities = shape_input.shapes() shape_attributes = shape_input.records() @@ -218,6 +231,8 @@ def plot_glacier_inventory(ax, base_dir, START=0, END=30, LINEWIDTH=0.6): # PURPOSE plot coastlines and islands (GSHHS with G250 Greenland) def plot_coastline(ax, base_dir): + # get logger + logger = logging.getLogger(__name__) # read the coastline shape file coastline_dir = base_dir.joinpath('masks', 'G250') coastline_shape_files = [] @@ -225,7 +240,7 @@ def plot_coastline(ax, base_dir): coastline_shape_files.append('greenland_coastline_islands.shp') for fi, S in zip(coastline_shape_files, [1000, 200]): coast_shapefile = coastline_dir.joinpath(fi) - logging.debug(str(coast_shapefile)) + logger.debug(str(coast_shapefile)) shape_input = shapefile.Reader(str(coast_shapefile)) shape_entities = shape_input.shapes() # for each entity within the shapefile @@ -237,9 +252,11 @@ def plot_coastline(ax, base_dir): # plot the MODIS Mosaic of Greenland as a background image def plot_image_mosaic(ax, base_dir, MASKED=True): + # get logger + logger = logging.getLogger(__name__) # read MODIS mosaic of Greenland image_geotiff_file = base_dir.joinpath(*image_file) - logging.debug(str(image_geotiff_file)) + logger.debug(str(image_geotiff_file)) ds = osgeo.gdal.Open(str(image_geotiff_file)) # get dimensions xsize = ds.RasterXSize @@ -289,11 +306,31 @@ def add_plot_scale(ax, X, Y, dx, dy, masked, fc1='w', fc2='k'): Y - 0.5 * dy, Y + 4.5 * dy, ] - ax.fill([x1, x2, x2, x1, x1], [y1, y1, y2, y2, y1], fc1, zorder=1) + ax.fill( + [x1, x2, x2, x1, x1], + [y1, y1, y2, y2, y1], + fc1, + zorder=1, + ) for i, c in enumerate([fc1, fc2, fc1, fc2, fc1]): - x1, x2, y1, y2 = [X + 0.2 * i * dx, X + 0.2 * (i + 1) * dx, Y, Y + dy] - ax.fill([x1, x2, x2, x1, x1], [y1, y1, y2, y2, y1], c, zorder=5) - ax.plot([X, X + dx, X + dx, X, X], [Y, Y, Y + dy, Y + dy, Y], fc2, zorder=4) + x1, x2, y1, y2 = [ + X + 0.2 * i * dx, + X + 0.2 * (i + 1) * dx, + Y, + Y + dy, + ] + ax.fill( + [x1, x2, x2, x1, x1], + [y1, y1, y2, y2, y1], + c, + zorder=5, + ) + ax.plot( + [X, X + dx, X + dx, X, X], + [Y, Y, Y + dy, Y + dy, Y], + fc2, + zorder=4, + ) ax.plot( [X, X, X + dx, X + dx], [Y + 1.5 * dy, Y, Y, Y + 1.5 * dy], @@ -348,6 +385,8 @@ def plot_grid( FIGURE_DPI=None, MODE=0o775, ): + # get logger + logger = logging.getLogger(__name__) # extend list if a single format was entered for all files if len(DATAFORM) < len(FILENAMES): DATAFORM = DATAFORM * len(FILENAMES) @@ -696,7 +735,7 @@ def plot_grid( # create output directory if non-existent FIGURE_FILE.parent.mkdir(mode=MODE, parents=True, exist_ok=True) # save to file - logging.info(str(FIGURE_FILE)) + logger.info(str(FIGURE_FILE)) plt.savefig( FIGURE_FILE, metadata={'Title': pathlib.Path(sys.argv[0]).name}, @@ -958,7 +997,9 @@ def main(): # create logger loglevels = [logging.CRITICAL, logging.INFO, logging.DEBUG] - logging.basicConfig(level=loglevels[args.verbose]) + logger = gravtk.utilities.build_logger( + __name__, level=loglevels[args.verbose] + ) # try to run the analysis with listed parameters try: @@ -1002,8 +1043,8 @@ def main(): # if there has been an error exception # print the type, value, and stack trace of the # current exception being handled - logging.critical(f'process id {os.getpid():d} failed') - logging.error(traceback.format_exc()) + logger.critical(f'process id {os.getpid():d} failed') + logger.error(traceback.format_exc()) # run main program diff --git a/gravity_toolkit/mapping/plot_GrIS_grid_maps.py b/gravity_toolkit/mapping/plot_GrIS_grid_maps.py index 729ff5a..ac238bc 100644 --- a/gravity_toolkit/mapping/plot_GrIS_grid_maps.py +++ b/gravity_toolkit/mapping/plot_GrIS_grid_maps.py @@ -1,7 +1,7 @@ #!/usr/bin/env python """ plot_GrIS_grid_maps.py -Written by Tyler Sutterley (05/2023) +Written by Tyler Sutterley (08/2026) Creates GMT-like plots for the Greenland ice sheet on a NSIDC polar stereographic north (3413) projection @@ -34,6 +34,7 @@ https://pypi.python.org/pypi/GDAL/ UPDATE HISTORY: + Updated 08/2026: use upstream file logger for verbose output Updated 05/2023: use pathlib to define and operate on paths added option to set the input variable names or column order Updated 03/2023: switch from parameter files to argparse arguments @@ -145,21 +146,26 @@ # PURPOSE: keep track of threads def info(args): - logging.info(pathlib.Path(sys.argv[0]).name) - logging.info(args) - logging.info(f'module name: {__name__}') + # get logger + logger = logging.getLogger(__name__) + logger.info(pathlib.Path(sys.argv[0]).name) + logger.info(args) + logger.info(f'module name: {__name__}') if hasattr(os, 'getppid'): - logging.info(f'parent process: {os.getppid():d}') - logging.info(f'process id: {os.getpid():d}') + logger.info(f'parent process: {os.getppid():d}') + logger.info(f'process id: {os.getpid():d}') # PURPOSE: plot Rignot 2012 drainage basin polylines def plot_rignot_basins(ax, base_dir): + # get logger + logger = logging.getLogger(__name__) region_directory = base_dir.joinpath(*region_dir) # for each region for reg in region_title: # read the regional polylines region_file = region_directory.joinpath(region_filename.format(reg)) + logger.debug(str(region_file)) region_ll = np.loadtxt(region_file, dtype=region_dtype) # converting region lat/lon into plot coordinates points = projection.transform_points( @@ -170,9 +176,11 @@ def plot_rignot_basins(ax, base_dir): # PURPOSE: plot Greenland drainage basins from IMBIE2 (Mouginot) def plot_IMBIE2_basins(ax, base_dir): + # get logger + logger = logging.getLogger(__name__) # read drainage basin polylines from shapefile (using splat operator) basin_shapefile = base_dir.joinpath(*IMBIE_basin_file) - logging.debug(str(basin_shapefile)) + logger.debug(str(basin_shapefile)) shape_input = shapefile.Reader(str(basin_shapefile)) shape_entities = shape_input.shapes() shape_attributes = shape_input.records() @@ -198,8 +206,10 @@ def plot_IMBIE2_basins(ax, base_dir): # PURPOSE: plot Greenland grounded ice delineation from GIMP def plot_grounded_ice(ax, base_dir, START=1, END=300, LINEWIDTH=0.6): + # get logger + logger = logging.getLogger(__name__) coast_shapefile = base_dir.joinpath(*coast_file) - logging.debug(str(coast_shapefile)) + logger.debug(str(coast_shapefile)) shape_input = shapefile.Reader(str(coast_shapefile)) shape_entities = shape_input.shapes() shape_attributes = shape_input.records() @@ -218,6 +228,9 @@ def plot_grounded_ice(ax, base_dir, START=1, END=300, LINEWIDTH=0.6): # PURPOSE: plot glaciated regions from Randolph Glacier Inventory def plot_glacier_inventory(ax, base_dir, START=0, END=30, LINEWIDTH=0.6): + # get logger + logger = logging.getLogger(__name__) + # RGI shapefiles for Arctic glaciers within plot region RGI_files = [] RGI_files.append('03_rgi60_ArcticCanadaNorth') RGI_files.append('04_rgi60_ArcticCanadaSouth') @@ -225,7 +238,7 @@ def plot_glacier_inventory(ax, base_dir, START=0, END=30, LINEWIDTH=0.6): RGI_files.append('07_rgi60_Svalbard') for f in RGI_files: RGI_shapefile = base_dir.joinpath('RGI', f, f'{f}_plot.shp') - logging.debug(str(RGI_shapefile)) + logger.debug(str(RGI_shapefile)) shape_input = shapefile.Reader(str(RGI_shapefile)) shape_entities = shape_input.shapes() shape_attributes = shape_input.records() @@ -244,6 +257,8 @@ def plot_glacier_inventory(ax, base_dir, START=0, END=30, LINEWIDTH=0.6): # PURPOSE plot coastlines and islands (GSHHS with G250 Greenland) def plot_coastline(ax, base_dir): + # get logger + logger = logging.getLogger(__name__) # read the coastline shape file coastline_dir = base_dir.joinpath('masks', 'G250') coastline_shape_files = [] @@ -251,7 +266,7 @@ def plot_coastline(ax, base_dir): coastline_shape_files.append('greenland_coastline_islands.shp') for fi, S in zip(coastline_shape_files, [1000, 200]): coast_shapefile = coastline_dir.joinpath(fi) - logging.debug(str(coast_shapefile)) + logger.debug(str(coast_shapefile)) shape_input = shapefile.Reader(str(coast_shapefile)) shape_entities = shape_input.shapes() # for each entity within the shapefile @@ -263,9 +278,11 @@ def plot_coastline(ax, base_dir): # plot the MODIS Mosaic of Greenland as a background image def plot_image_mosaic(ax, base_dir, MASKED=True): + # get logger + logger = logging.getLogger(__name__) # read MODIS mosaic of Greenland image_geotiff_file = base_dir.joinpath(*image_file) - logging.debug(str(image_geotiff_file)) + logger.debug(str(image_geotiff_file)) ds = osgeo.gdal.Open(str(image_geotiff_file)) # get dimensions xsize = ds.RasterXSize @@ -315,11 +332,31 @@ def add_plot_scale(ax, X, Y, dx, dy, masked, fc1='w', fc2='k'): Y - 2.5 * dy, Y + 3.2 * dy, ] - ax.fill([x1, x2, x2, x1, x1], [y1, y1, y2, y2, y1], fc1, zorder=4) + ax.fill( + [x1, x2, x2, x1, x1], + [y1, y1, y2, y2, y1], + fc1, + zorder=4, + ) for i, c in enumerate([fc1, fc2, fc1, fc2]): - x1, x2, y1, y2 = [X + 0.25 * i * dx, X + 0.25 * (i + 1) * dx, Y, Y + dy] - ax.fill([x1, x2, x2, x1, x1], [y1, y1, y2, y2, y1], c, zorder=5) - ax.plot([X, X + dx, X + dx, X, X], [Y, Y, Y + dy, Y + dy, Y], fc2, zorder=6) + x1, x2, y1, y2 = [ + X + 0.25 * i * dx, + X + 0.25 * (i + 1) * dx, + Y, + Y + dy, + ] + ax.fill( + [x1, x2, x2, x1, x1], + [y1, y1, y2, y2, y1], + c, + zorder=5, + ) + ax.plot( + [X, X + dx, X + dx, X, X], + [Y, Y, Y + dy, Y + dy, Y], + fc2, + zorder=6, + ) for i in range(3): ax.plot( [X + 0.5 * i * dx, X + 0.5 * i * dx], @@ -385,6 +422,8 @@ def plot_grid( FIGURE_DPI=None, MODE=0o775, ): + # get logger + logger = logging.getLogger(__name__) # read CPT or use color map if CPT_FILE is not None: # cpt file @@ -729,7 +768,7 @@ def plot_grid( # create output directory if non-existent FIGURE_FILE.parent.mkdir(mode=MODE, parents=True, exist_ok=True) # save to file - logging.info(str(FIGURE_FILE)) + logger.info(str(FIGURE_FILE)) plt.savefig( FIGURE_FILE, metadata={'Title': pathlib.Path(sys.argv[0]).name}, @@ -988,7 +1027,9 @@ def main(): # create logger loglevels = [logging.CRITICAL, logging.INFO, logging.DEBUG] - logging.basicConfig(level=loglevels[args.verbose]) + logger = gravtk.utilities.build_logger( + __name__, level=loglevels[args.verbose] + ) # try to run the analysis with listed parameters try: @@ -1032,8 +1073,8 @@ def main(): # if there has been an error exception # print the type, value, and stack trace of the # current exception being handled - logging.critical(f'process id {os.getpid():d} failed') - logging.error(traceback.format_exc()) + logger.critical(f'process id {os.getpid():d} failed') + logger.error(traceback.format_exc()) # run main program diff --git a/gravity_toolkit/mapping/plot_GrIS_grid_movie.py b/gravity_toolkit/mapping/plot_GrIS_grid_movie.py index 3f13ee5..17b04b9 100644 --- a/gravity_toolkit/mapping/plot_GrIS_grid_movie.py +++ b/gravity_toolkit/mapping/plot_GrIS_grid_movie.py @@ -1,7 +1,7 @@ #!/usr/bin/env python """ plot_GrIS_grid_movie.py -Written by Tyler Sutterley (05/2023) +Written by Tyler Sutterley (08/2026) Creates GMT-like animations for the Greenland Ice Sheet on a NSIDC polar stereographic north (3413) projection @@ -34,6 +34,7 @@ https://pypi.python.org/pypi/GDAL/ UPDATE HISTORY: + Updated 08/2026: use upstream file logger for verbose output Updated 05/2023: use pathlib to define and operate on paths Updated 03/2023: switch from parameter files to argparse arguments updated inputs to spatial from_file function @@ -158,21 +159,26 @@ # PURPOSE: keep track of threads def info(args): - logging.info(pathlib.Path(sys.argv[0]).name) - logging.info(args) - logging.info(f'module name: {__name__}') + # get logger + logger = logging.getLogger(__name__) + logger.info(pathlib.Path(sys.argv[0]).name) + logger.info(args) + logger.info(f'module name: {__name__}') if hasattr(os, 'getppid'): - logging.info(f'parent process: {os.getppid():d}') - logging.info(f'process id: {os.getpid():d}') + logger.info(f'parent process: {os.getppid():d}') + logger.info(f'process id: {os.getpid():d}') # PURPOSE: plot Rignot 2012 drainage basin polylines def plot_rignot_basins(ax, base_dir): + # get logger + logger = logging.getLogger(__name__) region_directory = base_dir.joinpath(*region_dir) # for each region for reg in region_title: # read the regional polylines region_file = region_directory.joinpath(region_filename.format(reg)) + logger.debug(str(region_file)) region_ll = np.loadtxt(region_file, dtype=region_dtype) # converting region lat/lon into plot coordinates points = projection.transform_points( @@ -183,9 +189,11 @@ def plot_rignot_basins(ax, base_dir): # PURPOSE: plot Greenland drainage basins from IMBIE2 (Mouginot) def plot_IMBIE2_basins(ax, base_dir): + # get logger + logger = logging.getLogger(__name__) # read drainage basin polylines from shapefile (using splat operator) basin_shapefile = base_dir.joinpath(*IMBIE_basin_file) - logging.debug(str(basin_shapefile)) + logger.debug(str(basin_shapefile)) shape_input = shapefile.Reader(str(basin_shapefile)) shape_entities = shape_input.shapes() shape_attributes = shape_input.records() @@ -211,8 +219,10 @@ def plot_IMBIE2_basins(ax, base_dir): # PURPOSE: plot Greenland grounded ice delineation from GIMP def plot_grounded_ice(ax, base_dir, START=1, END=300, LINEWIDTH=0.6): + # get logger + logger = logging.getLogger(__name__) coast_shapefile = base_dir.joinpath(*coast_file) - logging.debug(str(coast_shapefile)) + logger.debug(str(coast_shapefile)) shape_input = shapefile.Reader(str(coast_shapefile)) shape_entities = shape_input.shapes() shape_attributes = shape_input.records() @@ -231,6 +241,9 @@ def plot_grounded_ice(ax, base_dir, START=1, END=300, LINEWIDTH=0.6): # PURPOSE: plot glaciated regions from Randolph Glacier Inventory def plot_glacier_inventory(ax, base_dir, START=0, END=30, LINEWIDTH=0.6): + # get logger + logger = logging.getLogger(__name__) + # RGI shapefiles for Arctic glaciers within plot region RGI_files = [] RGI_files.append('03_rgi60_ArcticCanadaNorth') RGI_files.append('04_rgi60_ArcticCanadaSouth') @@ -238,7 +251,7 @@ def plot_glacier_inventory(ax, base_dir, START=0, END=30, LINEWIDTH=0.6): RGI_files.append('07_rgi60_Svalbard') for f in RGI_files: RGI_shapefile = base_dir.joinpath('RGI', f, f'{f}_plot.shp') - logging.debug(str(RGI_shapefile)) + logger.debug(str(RGI_shapefile)) shape_input = shapefile.Reader(str(RGI_shapefile)) shape_entities = shape_input.shapes() shape_attributes = shape_input.records() @@ -257,6 +270,8 @@ def plot_glacier_inventory(ax, base_dir, START=0, END=30, LINEWIDTH=0.6): # PURPOSE plot coastlines and islands (GSHHS with G250 Greenland) def plot_coastline(ax, base_dir): + # get logger + logger = logging.getLogger(__name__) # read the coastline shape file coastline_dir = base_dir.joinpath('masks', 'G250') coastline_shape_files = [] @@ -264,7 +279,7 @@ def plot_coastline(ax, base_dir): coastline_shape_files.append('greenland_coastline_islands.shp') for fi, S in zip(coastline_shape_files, [1000, 200]): coast_shapefile = coastline_dir.joinpath(fi) - logging.debug(str(coast_shapefile)) + logger.debug(str(coast_shapefile)) shape_input = shapefile.Reader(str(coast_shapefile)) shape_entities = shape_input.shapes() # for each entity within the shapefile @@ -276,9 +291,11 @@ def plot_coastline(ax, base_dir): # plot the MODIS Mosaic of Greenland as a background image def plot_image_mosaic(ax, base_dir, MASKED=True): + # get logger + logger = logging.getLogger(__name__) # read MODIS mosaic of Greenland image_geotiff_file = base_dir.joinpath(*image_file) - logging.debug(str(image_geotiff_file)) + logger.debug(str(image_geotiff_file)) ds = osgeo.gdal.Open(str(image_geotiff_file)) # get dimensions xsize = ds.RasterXSize @@ -328,11 +345,31 @@ def add_plot_scale(ax, X, Y, dx, dy, masked, fc1='w', fc2='k'): Y - 2.5 * dy, Y + 3.2 * dy, ] - ax.fill([x1, x2, x2, x1, x1], [y1, y1, y2, y2, y1], fc1, zorder=4) + ax.fill( + [x1, x2, x2, x1, x1], + [y1, y1, y2, y2, y1], + fc1, + zorder=4, + ) for i, c in enumerate([fc1, fc2, fc1, fc2]): - x1, x2, y1, y2 = [X + 0.25 * i * dx, X + 0.25 * (i + 1) * dx, Y, Y + dy] - ax.fill([x1, x2, x2, x1, x1], [y1, y1, y2, y2, y1], c, zorder=5) - ax.plot([X, X + dx, X + dx, X, X], [Y, Y, Y + dy, Y + dy, Y], fc2, zorder=6) + x1, x2, y1, y2 = [ + X + 0.25 * i * dx, + X + 0.25 * (i + 1) * dx, + Y, + Y + dy, + ] + ax.fill( + [x1, x2, x2, x1, x1], + [y1, y1, y2, y2, y1], + c, + zorder=5, + ) + ax.plot( + [X, X + dx, X + dx, X, X], + [Y, Y, Y + dy, Y + dy, Y], + fc2, + zorder=6, + ) for i in range(3): ax.plot( [X + 0.5 * i * dx, X + 0.5 * i * dx], @@ -397,6 +434,8 @@ def animate_grid( FIGURE_DPI=None, MODE=0o775, ): + # get logger + logger = logging.getLogger(__name__) # read CPT or use color map if CPT_FILE is not None: # cpt file @@ -1045,7 +1084,9 @@ def main(): # create logger loglevels = [logging.CRITICAL, logging.INFO, logging.DEBUG] - logging.basicConfig(level=loglevels[args.verbose]) + logger = gravtk.utilities.build_logger( + __name__, level=loglevels[args.verbose] + ) # try to run the analysis with listed parameters try: @@ -1087,8 +1128,8 @@ def main(): # if there has been an error exception # print the type, value, and stack trace of the # current exception being handled - logging.critical(f'process id {os.getpid():d} failed') - logging.error(traceback.format_exc()) + logger.critical(f'process id {os.getpid():d} failed') + logger.error(traceback.format_exc()) # run main program diff --git a/gravity_toolkit/mapping/plot_QML_grid_3maps.py b/gravity_toolkit/mapping/plot_QML_grid_3maps.py index 25eaad3..1478fc1 100644 --- a/gravity_toolkit/mapping/plot_QML_grid_3maps.py +++ b/gravity_toolkit/mapping/plot_QML_grid_3maps.py @@ -1,7 +1,7 @@ #!/usr/bin/env python """ plot_AIS_grid_3maps.py -Written by Tyler Sutterley (05/2023) +Written by Tyler Sutterley (08/2026) Creates 3 GMT-like plots for Queen Maud Land (QML) in Antarctica on a polar stereographic south (3031) projection @@ -31,6 +31,7 @@ https://pypi.python.org/pypi/GDAL/ UPDATE HISTORY: + Updated 08/2026: use upstream file logger for verbose output Updated 05/2023: use pathlib to define and operate on paths added option to set the input variable names or column order Updated 03/2023: switch from parameter files to argparse arguments @@ -176,21 +177,26 @@ # PURPOSE: keep track of threads def info(args): - logging.info(pathlib.Path(sys.argv[0]).name) - logging.info(args) - logging.info(f'module name: {__name__}') + # get logger + logger = logging.getLogger(__name__) + logger.info(pathlib.Path(sys.argv[0]).name) + logger.info(args) + logger.info(f'module name: {__name__}') if hasattr(os, 'getppid'): - logging.info(f'parent process: {os.getppid():d}') - logging.info(f'process id: {os.getpid():d}') + logger.info(f'parent process: {os.getppid():d}') + logger.info(f'process id: {os.getpid():d}') # PURPOSE: plot Rignot 2012 drainage basin polylines def plot_rignot_basins(ax, base_dir): + # get logger + logger = logging.getLogger(__name__) region_directory = base_dir.joinpath(*region_dir) # for each region for reg in region_title: # read the regional polylines region_file = region_directory.joinpath(region_filename.format(reg)) + logger.debug(str(region_file)) region_ll = np.loadtxt(region_file, dtype=region_dtype) # converting region lat/lon into plot coordinates points = projection.transform_points( @@ -201,9 +207,11 @@ def plot_rignot_basins(ax, base_dir): # PURPOSE: plot Antarctic drainage basins from IMBIE2 (Mouginot) def plot_IMBIE2_basins(ax, base_dir): + # get logger + logger = logging.getLogger(__name__) # read drainage basin polylines from shapefile (using splat operator) basin_shapefile = base_dir.joinpath(*IMBIE_basin_file) - logging.debug(str(basin_shapefile)) + logger.debug(str(basin_shapefile)) shape_input = shapefile.Reader(str(basin_shapefile)) shape_entities = shape_input.shapes() shape_attributes = shape_input.records() @@ -225,10 +233,12 @@ def plot_IMBIE2_basins(ax, base_dir): # PURPOSE: plot Antarctic drainage sub-basins from IMBIE-2 (Mouginot) def plot_IMBIE2_subbasins(ax, base_dir): + # get logger + logger = logging.getLogger(__name__) # read drainage basin polylines from shapefile (using splat operator) IMBIE_subbasin_file = ['Basins_20Oct2016_v1.7', 'Basins_v1.7.shp'] basin_shapefile = base_dir.joinpath('masks', *IMBIE_subbasin_file) - logging.debug(str(basin_shapefile)) + logger.debug(str(basin_shapefile)) shape_input = shapefile.Reader(str(basin_shapefile)) shape_entities = shape_input.shapes() shape_attributes = shape_input.records() @@ -248,8 +258,10 @@ def plot_IMBIE2_subbasins(ax, base_dir): # PURPOSE: plot Antarctic grounded ice delineation def plot_grounded_ice(ax, base_dir, START=1): + # get logger + logger = logging.getLogger(__name__) grounded_ice_shapefile = base_dir.joinpath(*coast_file) - logging.debug(str(grounded_ice_shapefile)) + logger.debug(str(grounded_ice_shapefile)) shape_input = shapefile.Reader(str(grounded_ice_shapefile)) shape_entities = shape_input.shapes() shape_attributes = shape_input.records() @@ -262,9 +274,11 @@ def plot_grounded_ice(ax, base_dir, START=1): # PURPOSE: plot MODIS mosaic of Antarctica as background image def plot_image_mosaic(ax, base_dir, MASKED=True): + # get logger + logger = logging.getLogger(__name__) # read MODIS mosaic of Antarctica image_geotiff_file = base_dir.joinpath(*image_file) - logging.debug(str(image_geotiff_file)) + logger.debug(str(image_geotiff_file)) ds = osgeo.gdal.Open(str(image_geotiff_file)) # get geotiff info info_geotiff = ds.GetGeoTransform() @@ -323,11 +337,31 @@ def add_plot_scale(ax, X, Y, dx, dy, masked, fc1='w', fc2='k'): Y - 3.5 * dy, Y + 3.4 * dy, ] - ax.fill([x1, x2, x2, x1, x1], [y1, y1, y2, y2, y1], fc1, zorder=4) + ax.fill( + [x1, x2, x2, x1, x1], + [y1, y1, y2, y2, y1], + fc1, + zorder=4, + ) for i, c in enumerate([fc1, fc2, fc1, fc2]): - x1, x2, y1, y2 = [X + 0.25 * i * dx, X + 0.25 * (i + 1) * dx, Y, Y + dy] - ax.fill([x1, x2, x2, x1, x1], [y1, y1, y2, y2, y1], c, zorder=5) - ax.plot([X, X + dx, X + dx, X, X], [Y, Y, Y + dy, Y + dy, Y], fc2, zorder=6) + x1, x2, y1, y2 = [ + X + 0.25 * i * dx, + X + 0.25 * (i + 1) * dx, + Y, + Y + dy, + ] + ax.fill( + [x1, x2, x2, x1, x1], + [y1, y1, y2, y2, y1], + c, + zorder=5, + ) + ax.plot( + [X, X + dx, X + dx, X, X], + [Y, Y, Y + dy, Y + dy, Y], + fc2, + zorder=6, + ) for i in range(3): ax.plot( [X + 0.5 * i * dx, X + 0.5 * i * dx], @@ -393,6 +427,8 @@ def plot_grid( FIGURE_DPI=None, MODE=0o775, ): + # get logger + logger = logging.getLogger(__name__) # extend list if a single format was entered for all files if len(DATAFORM) < len(FILENAMES): DATAFORM = DATAFORM * len(FILENAMES) @@ -740,7 +776,7 @@ def plot_grid( # create output directory if non-existent FIGURE_FILE.parent.mkdir(mode=MODE, parents=True, exist_ok=True) # save to file - logging.info(str(FIGURE_FILE)) + logger.info(str(FIGURE_FILE)) plt.savefig( FIGURE_FILE, metadata={'Title': pathlib.Path(sys.argv[0]).name}, @@ -996,7 +1032,9 @@ def main(): # create logger loglevels = [logging.CRITICAL, logging.INFO, logging.DEBUG] - logging.basicConfig(level=loglevels[args.verbose]) + logger = gravtk.utilities.build_logger( + __name__, level=loglevels[args.verbose] + ) # try to run the analysis with listed parameters try: @@ -1039,8 +1077,8 @@ def main(): # if there has been an error exception # print the type, value, and stack trace of the # current exception being handled - logging.critical(f'process id {os.getpid():d} failed') - logging.error(traceback.format_exc()) + logger.critical(f'process id {os.getpid():d} failed') + logger.error(traceback.format_exc()) # run main program diff --git a/gravity_toolkit/mapping/plot_global_grid_3maps.py b/gravity_toolkit/mapping/plot_global_grid_3maps.py index 5311085..127929c 100644 --- a/gravity_toolkit/mapping/plot_global_grid_3maps.py +++ b/gravity_toolkit/mapping/plot_global_grid_3maps.py @@ -1,7 +1,7 @@ #!/usr/bin/env python """ plot_global_grid_3maps.py -Written by Tyler Sutterley (05/2023) +Written by Tyler Sutterley (08/2026) Creates 3 GMT-like plots in a Plate Carree (Equirectangular) projection PYTHON DEPENDENCIES: @@ -23,6 +23,7 @@ https://github.com/GeospatialPython/pyshp UPDATE HISTORY: + Updated 08/2026: use upstream file logger for verbose output Updated 05/2023: use pathlib to define and operate on paths added option to set the input variable names or column order Updated 03/2023: switch from parameter files to argparse arguments @@ -95,16 +96,20 @@ # PURPOSE: keep track of threads def info(args): - logging.info(pathlib.Path(sys.argv[0]).name) - logging.info(args) - logging.info(f'module name: {__name__}') + # get logger + logger = logging.getLogger(__name__) + logger.info(pathlib.Path(sys.argv[0]).name) + logger.info(args) + logger.info(f'module name: {__name__}') if hasattr(os, 'getppid'): - logging.info(f'parent process: {os.getppid():d}') - logging.info(f'process id: {os.getpid():d}') + logger.info(f'parent process: {os.getppid():d}') + logger.info(f'process id: {os.getpid():d}') # PURPOSE plot coastlines and islands (GSHHS with G250 Greenland) def plot_coastline(ax, base_dir, LINEWIDTH=0.5): + # get logger + logger = logging.getLogger(__name__) # read the coastline shape file coastline_dir = base_dir.joinpath('masks', 'G250') coastline_shape_files = [] @@ -112,7 +117,7 @@ def plot_coastline(ax, base_dir, LINEWIDTH=0.5): coastline_shape_files.append('greenland_coastline_islands.shp') for fi, S in zip(coastline_shape_files, [1000, 200]): coast_shapefile = coastline_dir.joinpath(fi) - logging.debug(str(coast_shapefile)) + logger.debug(str(coast_shapefile)) shape_input = shapefile.Reader(str(coast_shapefile)) shape_entities = shape_input.shapes() # for each entity within the shapefile @@ -124,13 +129,16 @@ def plot_coastline(ax, base_dir, LINEWIDTH=0.5): # PURPOSE: plot Antarctic grounded ice delineation def plot_grounded_ice(ax, base_dir, LINEWIDTH=0.5): + # get logger + logger = logging.getLogger(__name__) + # path to shapefile for grounded ice delineation grounded_ice_file = [ 'masks', 'IceBoundaries_Antarctica_v02', 'ant_ice_sheet_islands_v2.shp', ] grounded_ice_shapefile = base_dir.joinpath(*grounded_ice_file) - logging.debug(str(grounded_ice_shapefile)) + logger.debug(str(grounded_ice_shapefile)) shape_input = shapefile.Reader(str(grounded_ice_shapefile)) shape_entities = shape_input.shapes() shape_attributes = shape_input.records() @@ -178,6 +186,8 @@ def plot_grid( FIGURE_DPI=None, MODE=0o775, ): + # get logger + logger = logging.getLogger(__name__) # extend list if a single format was entered for all files if len(DATAFORM) < len(FILENAMES): DATAFORM = DATAFORM * len(FILENAMES) @@ -502,7 +512,7 @@ def plot_grid( # create output directory if non-existent FIGURE_FILE.parent.mkdir(mode=MODE, parents=True, exist_ok=True) # save to file - logging.info(str(FIGURE_FILE)) + logger.info(str(FIGURE_FILE)) plt.savefig( FIGURE_FILE, metadata={'Title': pathlib.Path(sys.argv[0]).name}, @@ -740,7 +750,9 @@ def main(): # create logger loglevels = [logging.CRITICAL, logging.INFO, logging.DEBUG] - logging.basicConfig(level=loglevels[args.verbose]) + logger = gravtk.utilities.build_logger( + __name__, level=loglevels[args.verbose] + ) # try to run the analysis with listed parameters try: @@ -780,8 +792,8 @@ def main(): # if there has been an error exception # print the type, value, and stack trace of the # current exception being handled - logging.critical(f'process id {os.getpid():d} failed') - logging.error(traceback.format_exc()) + logger.critical(f'process id {os.getpid():d} failed') + logger.error(traceback.format_exc()) # run main program diff --git a/gravity_toolkit/mapping/plot_global_grid_4maps.py b/gravity_toolkit/mapping/plot_global_grid_4maps.py index 5fb99dd..0af72af 100644 --- a/gravity_toolkit/mapping/plot_global_grid_4maps.py +++ b/gravity_toolkit/mapping/plot_global_grid_4maps.py @@ -1,7 +1,7 @@ #!/usr/bin/env python """ plot_global_grid_4maps.py -Written by Tyler Sutterley (05/2023) +Written by Tyler Sutterley (08/2026) Creates 4 GMT-like plots in a Plate Carree (Equirectangular) projection PYTHON DEPENDENCIES: @@ -23,6 +23,7 @@ https://github.com/GeospatialPython/pyshp UPDATE HISTORY: + Updated 08/2026: use upstream file logger for verbose output Updated 05/2023: use pathlib to define and operate on paths added option to set the input variable names or column order Updated 03/2023: switch from parameter files to argparse arguments @@ -95,16 +96,20 @@ # PURPOSE: keep track of threads def info(args): - logging.info(pathlib.Path(sys.argv[0]).name) - logging.info(args) - logging.info(f'module name: {__name__}') + # get logger + logger = logging.getLogger(__name__) + logger.info(pathlib.Path(sys.argv[0]).name) + logger.info(args) + logger.info(f'module name: {__name__}') if hasattr(os, 'getppid'): - logging.info(f'parent process: {os.getppid():d}') - logging.info(f'process id: {os.getpid():d}') + logger.info(f'parent process: {os.getppid():d}') + logger.info(f'process id: {os.getpid():d}') # PURPOSE plot coastlines and islands (GSHHS with G250 Greenland) def plot_coastline(ax, base_dir, LINEWIDTH=0.5): + # get logger + logger = logging.getLogger(__name__) # read the coastline shape file coastline_dir = base_dir.joinpath('masks', 'G250') coastline_shape_files = [] @@ -112,7 +117,7 @@ def plot_coastline(ax, base_dir, LINEWIDTH=0.5): coastline_shape_files.append('greenland_coastline_islands.shp') for fi, S in zip(coastline_shape_files, [1000, 200]): coast_shapefile = coastline_dir.joinpath(fi) - logging.debug(str(coast_shapefile)) + logger.debug(str(coast_shapefile)) shape_input = shapefile.Reader(str(coast_shapefile)) shape_entities = shape_input.shapes() # for each entity within the shapefile @@ -124,13 +129,16 @@ def plot_coastline(ax, base_dir, LINEWIDTH=0.5): # PURPOSE: plot Antarctic grounded ice delineation def plot_grounded_ice(ax, base_dir, LINEWIDTH=0.5): + # get logger + logger = logging.getLogger(__name__) + # path to shapefile for grounded ice delineation grounded_ice_file = [ 'masks', 'IceBoundaries_Antarctica_v02', 'ant_ice_sheet_islands_v2.shp', ] grounded_ice_shapefile = base_dir.joinpath(*grounded_ice_file) - logging.debug(str(grounded_ice_shapefile)) + logger.debug(str(grounded_ice_shapefile)) shape_input = shapefile.Reader(str(grounded_ice_shapefile)) shape_entities = shape_input.shapes() shape_attributes = shape_input.records() @@ -178,6 +186,8 @@ def plot_grid( FIGURE_DPI=None, MODE=0o775, ): + # get logger + logger = logging.getLogger(__name__) # extend list if a single format was entered for all files if len(DATAFORM) < len(FILENAMES): DATAFORM = DATAFORM * len(FILENAMES) @@ -509,7 +519,7 @@ def plot_grid( # create output directory if non-existent FIGURE_FILE.parent.mkdir(mode=MODE, parents=True, exist_ok=True) # save to file - logging.info(str(FIGURE_FILE)) + logger.info(str(FIGURE_FILE)) plt.savefig( FIGURE_FILE, metadata={'Title': pathlib.Path(sys.argv[0]).name}, @@ -747,7 +757,9 @@ def main(): # create logger loglevels = [logging.CRITICAL, logging.INFO, logging.DEBUG] - logging.basicConfig(level=loglevels[args.verbose]) + logger = gravtk.utilities.build_logger( + __name__, level=loglevels[args.verbose] + ) # try to run the analysis with listed parameters try: @@ -787,8 +799,8 @@ def main(): # if there has been an error exception # print the type, value, and stack trace of the # current exception being handled - logging.critical(f'process id {os.getpid():d} failed') - logging.error(traceback.format_exc()) + logger.critical(f'process id {os.getpid():d} failed') + logger.error(traceback.format_exc()) # run main program diff --git a/gravity_toolkit/mapping/plot_global_grid_5maps.py b/gravity_toolkit/mapping/plot_global_grid_5maps.py index 9e770b6..6d0b9a7 100644 --- a/gravity_toolkit/mapping/plot_global_grid_5maps.py +++ b/gravity_toolkit/mapping/plot_global_grid_5maps.py @@ -1,7 +1,7 @@ #!/usr/bin/env python """ plot_global_grid_5maps.py -Written by Tyler Sutterley (05/2023) +Written by Tyler Sutterley (08/2026) Creates 5 GMT-like plots in a Plate Carree (Equirectangular) projection PYTHON DEPENDENCIES: @@ -23,6 +23,7 @@ https://github.com/GeospatialPython/pyshp UPDATE HISTORY: + Updated 08/2026: use upstream file logger for verbose output Updated 05/2023: use pathlib to define and operate on paths added option to set the input variable names or column order Updated 03/2023: switch from parameter files to argparse arguments @@ -96,16 +97,20 @@ # PURPOSE: keep track of threads def info(args): - logging.info(pathlib.Path(sys.argv[0]).name) - logging.info(args) - logging.info(f'module name: {__name__}') + # get logger + logger = logging.getLogger(__name__) + logger.info(pathlib.Path(sys.argv[0]).name) + logger.info(args) + logger.info(f'module name: {__name__}') if hasattr(os, 'getppid'): - logging.info(f'parent process: {os.getppid():d}') - logging.info(f'process id: {os.getpid():d}') + logger.info(f'parent process: {os.getppid():d}') + logger.info(f'process id: {os.getpid():d}') # PURPOSE plot coastlines and islands (GSHHS with G250 Greenland) def plot_coastline(ax, base_dir, LINEWIDTH=0.5): + # get logger + logger = logging.getLogger(__name__) # read the coastline shape file coastline_dir = base_dir.joinpath('masks', 'G250') coastline_shape_files = [] @@ -113,7 +118,7 @@ def plot_coastline(ax, base_dir, LINEWIDTH=0.5): coastline_shape_files.append('greenland_coastline_islands.shp') for fi, S in zip(coastline_shape_files, [1000, 200]): coast_shapefile = coastline_dir.joinpath(fi) - logging.debug(str(coast_shapefile)) + logger.debug(str(coast_shapefile)) shape_input = shapefile.Reader(str(coast_shapefile)) shape_entities = shape_input.shapes() # for each entity within the shapefile @@ -125,13 +130,16 @@ def plot_coastline(ax, base_dir, LINEWIDTH=0.5): # PURPOSE: plot Antarctic grounded ice delineation def plot_grounded_ice(ax, base_dir, LINEWIDTH=0.5): + # get logger + logger = logging.getLogger(__name__) + # path to shapefile for grounded ice delineation grounded_ice_file = [ 'masks', 'IceBoundaries_Antarctica_v02', 'ant_ice_sheet_islands_v2.shp', ] grounded_ice_shapefile = base_dir.joinpath(*grounded_ice_file) - logging.debug(str(grounded_ice_shapefile)) + logger.debug(str(grounded_ice_shapefile)) shape_input = shapefile.Reader(str(grounded_ice_shapefile)) shape_entities = shape_input.shapes() shape_attributes = shape_input.records() @@ -179,6 +187,8 @@ def plot_grid( FIGURE_DPI=None, MODE=0o775, ): + # get logger + logger = logging.getLogger(__name__) # extend list if a single format was entered for all files if len(DATAFORM) < len(FILENAMES): DATAFORM = DATAFORM * len(FILENAMES) @@ -510,7 +520,7 @@ def plot_grid( # create output directory if non-existent FIGURE_FILE.parent.mkdir(mode=MODE, parents=True, exist_ok=True) # save to file - logging.info(str(FIGURE_FILE)) + logger.info(str(FIGURE_FILE)) plt.savefig( FIGURE_FILE, metadata={'Title': pathlib.Path(sys.argv[0]).name}, @@ -748,7 +758,9 @@ def main(): # create logger loglevels = [logging.CRITICAL, logging.INFO, logging.DEBUG] - logging.basicConfig(level=loglevels[args.verbose]) + logger = gravtk.utilities.build_logger( + __name__, level=loglevels[args.verbose] + ) # try to run the analysis with listed parameters try: @@ -788,8 +800,8 @@ def main(): # if there has been an error exception # print the type, value, and stack trace of the # current exception being handled - logging.critical(f'process id {os.getpid():d} failed') - logging.error(traceback.format_exc()) + logger.critical(f'process id {os.getpid():d} failed') + logger.error(traceback.format_exc()) # run main program diff --git a/gravity_toolkit/mapping/plot_global_grid_9maps.py b/gravity_toolkit/mapping/plot_global_grid_9maps.py index 771ffce..f9a1f62 100644 --- a/gravity_toolkit/mapping/plot_global_grid_9maps.py +++ b/gravity_toolkit/mapping/plot_global_grid_9maps.py @@ -1,7 +1,7 @@ #!/usr/bin/env python """ plot_global_grid_9maps.py -Written by Tyler Sutterley (05/2023) +Written by Tyler Sutterley (08/2026) Creates 9 GMT-like plots in a Plate Carree (Equirectangular) projection PYTHON DEPENDENCIES: @@ -23,6 +23,7 @@ https://github.com/GeospatialPython/pyshp UPDATE HISTORY: + Updated 08/2026: use upstream file logger for verbose output Updated 05/2023: use pathlib to define and operate on paths added option to set the input variable names or column order Updated 03/2023: switch from parameter files to argparse arguments @@ -95,16 +96,20 @@ # PURPOSE: keep track of threads def info(args): - logging.info(pathlib.Path(sys.argv[0]).name) - logging.info(args) - logging.info(f'module name: {__name__}') + # get logger + logger = logging.getLogger(__name__) + logger.info(pathlib.Path(sys.argv[0]).name) + logger.info(args) + logger.info(f'module name: {__name__}') if hasattr(os, 'getppid'): - logging.info(f'parent process: {os.getppid():d}') - logging.info(f'process id: {os.getpid():d}') + logger.info(f'parent process: {os.getppid():d}') + logger.info(f'process id: {os.getpid():d}') # PURPOSE plot coastlines and islands (GSHHS with G250 Greenland) def plot_coastline(ax, base_dir, LINEWIDTH=0.5): + # get logger + logger = logging.getLogger(__name__) # read the coastline shape file coastline_dir = base_dir.joinpath('masks', 'G250') coastline_shape_files = [] @@ -112,7 +117,7 @@ def plot_coastline(ax, base_dir, LINEWIDTH=0.5): coastline_shape_files.append('greenland_coastline_islands.shp') for fi, S in zip(coastline_shape_files, [1000, 200]): coast_shapefile = coastline_dir.joinpath(fi) - logging.debug(str(coast_shapefile)) + logger.debug(str(coast_shapefile)) shape_input = shapefile.Reader(str(coast_shapefile)) shape_entities = shape_input.shapes() # for each entity within the shapefile @@ -124,13 +129,16 @@ def plot_coastline(ax, base_dir, LINEWIDTH=0.5): # PURPOSE: plot Antarctic grounded ice delineation def plot_grounded_ice(ax, base_dir, LINEWIDTH=0.5): + # get logger + logger = logging.getLogger(__name__) + # path to shapefile for grounded ice delineation grounded_ice_file = [ 'masks', 'IceBoundaries_Antarctica_v02', 'ant_ice_sheet_islands_v2.shp', ] grounded_ice_shapefile = base_dir.joinpath(*grounded_ice_file) - logging.debug(str(grounded_ice_shapefile)) + logger.debug(str(grounded_ice_shapefile)) shape_input = shapefile.Reader(str(grounded_ice_shapefile)) shape_entities = shape_input.shapes() shape_attributes = shape_input.records() @@ -178,6 +186,8 @@ def plot_grid( FIGURE_DPI=None, MODE=0o775, ): + # get logger + logger = logging.getLogger(__name__) # extend list if a single format was entered for all files if len(DATAFORM) < len(FILENAMES): DATAFORM = DATAFORM * len(FILENAMES) @@ -512,7 +522,7 @@ def plot_grid( # create output directory if non-existent FIGURE_FILE.parent.mkdir(mode=MODE, parents=True, exist_ok=True) # save to file - logging.info(str(FIGURE_FILE)) + logger.info(str(FIGURE_FILE)) plt.savefig( FIGURE_FILE, metadata={'Title': pathlib.Path(sys.argv[0]).name}, @@ -750,7 +760,9 @@ def main(): # create logger loglevels = [logging.CRITICAL, logging.INFO, logging.DEBUG] - logging.basicConfig(level=loglevels[args.verbose]) + logger = gravtk.utilities.build_logger( + __name__, level=loglevels[args.verbose] + ) # try to run the analysis with listed parameters try: @@ -790,8 +802,8 @@ def main(): # if there has been an error exception # print the type, value, and stack trace of the # current exception being handled - logging.critical(f'process id {os.getpid():d} failed') - logging.error(traceback.format_exc()) + logger.critical(f'process id {os.getpid():d} failed') + logger.error(traceback.format_exc()) # run main program diff --git a/gravity_toolkit/mapping/plot_global_grid_maps.py b/gravity_toolkit/mapping/plot_global_grid_maps.py index 8fe8f96..f9cfc9d 100644 --- a/gravity_toolkit/mapping/plot_global_grid_maps.py +++ b/gravity_toolkit/mapping/plot_global_grid_maps.py @@ -1,7 +1,7 @@ #!/usr/bin/env python """ plot_global_grid_maps.py -Written by Tyler Sutterley (05/2023) +Written by Tyler Sutterley (08/2026) Creates GMT-like plots in a Plate Carree (Equirectangular) projection PYTHON DEPENDENCIES: @@ -23,6 +23,7 @@ https://github.com/GeospatialPython/pyshp UPDATE HISTORY: + Updated 08/2026: use upstream file logger for verbose output Updated 05/2023: use pathlib to define and operate on paths added option to set the input variable names or column order Updated 03/2023: switch from parameter files to argparse arguments @@ -96,16 +97,20 @@ # PURPOSE: keep track of threads def info(args): - logging.info(pathlib.Path(sys.argv[0]).name) - logging.info(args) - logging.info(f'module name: {__name__}') + # get logger + logger = logging.getLogger(__name__) + logger.info(pathlib.Path(sys.argv[0]).name) + logger.info(args) + logger.info(f'module name: {__name__}') if hasattr(os, 'getppid'): - logging.info(f'parent process: {os.getppid():d}') - logging.info(f'process id: {os.getpid():d}') + logger.info(f'parent process: {os.getppid():d}') + logger.info(f'process id: {os.getpid():d}') # PURPOSE plot coastlines and islands (GSHHS with G250 Greenland) def plot_coastline(ax, base_dir, LINEWIDTH=0.5): + # get logger + logger = logging.getLogger(__name__) # read the coastline shape file coastline_dir = base_dir.joinpath('masks', 'G250') coastline_shape_files = [] @@ -113,7 +118,7 @@ def plot_coastline(ax, base_dir, LINEWIDTH=0.5): coastline_shape_files.append('greenland_coastline_islands.shp') for fi, S in zip(coastline_shape_files, [1000, 200]): coast_shapefile = coastline_dir.joinpath(fi) - logging.debug(str(coast_shapefile)) + logger.debug(str(coast_shapefile)) shape_input = shapefile.Reader(str(coast_shapefile)) shape_entities = shape_input.shapes() # for each entity within the shapefile @@ -125,13 +130,16 @@ def plot_coastline(ax, base_dir, LINEWIDTH=0.5): # PURPOSE: plot Antarctic grounded ice delineation def plot_grounded_ice(ax, base_dir, LINEWIDTH=0.5): + # get logger + logger = logging.getLogger(__name__) + # path to shapefile for grounded ice delineation grounded_ice_file = [ 'masks', 'IceBoundaries_Antarctica_v02', 'ant_ice_sheet_islands_v2.shp', ] grounded_ice_shapefile = base_dir.joinpath(*grounded_ice_file) - logging.debug(str(grounded_ice_shapefile)) + logger.debug(str(grounded_ice_shapefile)) shape_input = shapefile.Reader(str(grounded_ice_shapefile)) shape_entities = shape_input.shapes() shape_attributes = shape_input.records() @@ -180,6 +188,8 @@ def plot_grid( FIGURE_DPI=None, MODE=0o775, ): + # get logger + logger = logging.getLogger(__name__) # read CPT or use color map if CPT_FILE is not None: # cpt file @@ -514,7 +524,7 @@ def plot_grid( # create output directory if non-existent FIGURE_FILE.parent.mkdir(mode=MODE, parents=True, exist_ok=True) # save to file - logging.info(str(FIGURE_FILE)) + logger.info(str(FIGURE_FILE)) plt.savefig( FIGURE_FILE, metadata={'Title': pathlib.Path(sys.argv[0]).name}, @@ -761,7 +771,9 @@ def main(): # create logger loglevels = [logging.CRITICAL, logging.INFO, logging.DEBUG] - logging.basicConfig(level=loglevels[args.verbose]) + logger = gravtk.utilities.build_logger( + __name__, level=loglevels[args.verbose] + ) # try to run the analysis with listed parameters try: @@ -803,8 +815,8 @@ def main(): # if there has been an error exception # print the type, value, and stack trace of the # current exception being handled - logging.critical(f'process id {os.getpid():d} failed') - logging.error(traceback.format_exc()) + logger.critical(f'process id {os.getpid():d} failed') + logger.error(traceback.format_exc()) # run main program diff --git a/gravity_toolkit/mapping/plot_global_grid_movie.py b/gravity_toolkit/mapping/plot_global_grid_movie.py index 0b212bd..f3f27eb 100644 --- a/gravity_toolkit/mapping/plot_global_grid_movie.py +++ b/gravity_toolkit/mapping/plot_global_grid_movie.py @@ -1,7 +1,7 @@ #!/usr/bin/env python """ plot_global_grid_maps.py -Written by Tyler Sutterley (05/2023) +Written by Tyler Sutterley (08/2026) Creates GMT-like animations in a Plate Carree (Equirectangular) projection PYTHON DEPENDENCIES: @@ -23,6 +23,7 @@ https://github.com/GeospatialPython/pyshp UPDATE HISTORY: + Updated 08/2026: use upstream file logger for verbose output Updated 05/2023: use pathlib to define and operate on paths Updated 03/2023: switch from parameter files to argparse arguments updated inputs to spatial from_file function @@ -98,16 +99,20 @@ # PURPOSE: keep track of threads def info(args): - logging.info(pathlib.Path(sys.argv[0]).name) - logging.info(args) - logging.info(f'module name: {__name__}') + # get logger + logger = logging.getLogger(__name__) + logger.info(pathlib.Path(sys.argv[0]).name) + logger.info(args) + logger.info(f'module name: {__name__}') if hasattr(os, 'getppid'): - logging.info(f'parent process: {os.getppid():d}') - logging.info(f'process id: {os.getpid():d}') + logger.info(f'parent process: {os.getppid():d}') + logger.info(f'process id: {os.getpid():d}') # PURPOSE plot coastlines and islands (GSHHS with G250 Greenland) def plot_coastline(ax, base_dir, LINEWIDTH=0.5): + # get logger + logger = logging.getLogger(__name__) # read the coastline shape file coastline_dir = base_dir.joinpath('masks', 'G250') coastline_shape_files = [] @@ -115,7 +120,7 @@ def plot_coastline(ax, base_dir, LINEWIDTH=0.5): coastline_shape_files.append('greenland_coastline_islands.shp') for fi, S in zip(coastline_shape_files, [1000, 200]): coast_shapefile = coastline_dir.joinpath(fi) - logging.debug(str(coast_shapefile)) + logger.debug(str(coast_shapefile)) shape_input = shapefile.Reader(str(coast_shapefile)) shape_entities = shape_input.shapes() # for each entity within the shapefile @@ -127,13 +132,16 @@ def plot_coastline(ax, base_dir, LINEWIDTH=0.5): # PURPOSE: plot Antarctic grounded ice delineation def plot_grounded_ice(ax, base_dir, LINEWIDTH=0.5): + # get logger + logger = logging.getLogger(__name__) + # path to shapefile for grounded ice delineation grounded_ice_file = [ 'masks', 'IceBoundaries_Antarctica_v02', 'ant_ice_sheet_islands_v2.shp', ] grounded_ice_shapefile = base_dir.joinpath(*grounded_ice_file) - logging.debug(str(grounded_ice_shapefile)) + logger.debug(str(grounded_ice_shapefile)) shape_input = shapefile.Reader(str(grounded_ice_shapefile)) shape_entities = shape_input.shapes() shape_attributes = shape_input.records() @@ -179,6 +187,8 @@ def animate_grid( FIGURE_DPI=None, MODE=0o775, ): + # get logger + logger = logging.getLogger(__name__) # read CPT or use color map if CPT_FILE is not None: # cpt file @@ -775,7 +785,9 @@ def main(): # create logger loglevels = [logging.CRITICAL, logging.INFO, logging.DEBUG] - logging.basicConfig(level=loglevels[args.verbose]) + logger = gravtk.utilities.build_logger( + __name__, level=loglevels[args.verbose] + ) # try to run the analysis with listed parameters try: @@ -813,8 +825,8 @@ def main(): # if there has been an error exception # print the type, value, and stack trace of the # current exception being handled - logging.critical(f'process id {os.getpid():d} failed') - logging.error(traceback.format_exc()) + logger.critical(f'process id {os.getpid():d} failed') + logger.error(traceback.format_exc()) # run main program diff --git a/gravity_toolkit/read_GRACE_harmonics.py b/gravity_toolkit/read_GRACE_harmonics.py index 70868ad..dd210eb 100644 --- a/gravity_toolkit/read_GRACE_harmonics.py +++ b/gravity_toolkit/read_GRACE_harmonics.py @@ -1,7 +1,7 @@ #!/usr/bin/env python """ read_GRACE_harmonics.py -Written by Tyler Sutterley (11/2024) +Written by Tyler Sutterley (08/2026) Contributions by Hugo Lecomte Reads GRACE files and extracts spherical harmonic data and drift rates (RL04) @@ -42,6 +42,7 @@ time.py: utilities for calculating time operations UPDATE HISTORY: + Updated 08/2026: use python datetime to calculate start and end dates Updated 11/2024: check if the GRACE/GRACE-FO files are gfc format Updated 05/2023: use pathlib to define and operate on paths Updated 03/2023: added regex formatting for CNES GRGS harmonics @@ -75,6 +76,7 @@ import pathlib import numpy as np import gravity_toolkit.time +from datetime import datetime, timedelta # PURPOSE: read Level-2 GRACE and GRACE-FO spherical harmonic files @@ -153,12 +155,19 @@ def read_GRACE_harmonics(input_file, LMAX, **kwargs): FLAG = r'GRCOF2' # output python dictionary with GRACE/GRACE-FO data and metadata - grace_L2_input = {} + # spherical harmonic model (SHM) data + SHM = {} + # extract GRACE/GRACE-FO date information from input file name - start_yr = np.float64(SY) - end_yr = np.float64(EY) - start_day = np.float64(SD) - end_day = np.float64(ED) + start_date = datetime(int(SY), 1, 1) + timedelta(days=int(SD) - 1) + start_struct = start_date.timetuple() + end_date = datetime(int(EY), 1, 1) + timedelta(days=int(ED) - 1) + end_struct = end_date.timetuple() + # start and end day of the year + start_yr = start_struct.tm_year + start_day = start_struct.tm_yday + end_yr = end_struct.tm_year + end_day = end_struct.tm_yday # calculate mid-month date taking into account if measurements are # on different years dpy = gravity_toolkit.time.calendar_days(start_yr).sum() @@ -168,40 +177,44 @@ def read_GRACE_harmonics(input_file, LMAX, **kwargs): # Calculate mid-month value mid_day = np.mean([start_day, end_cyclic]) # Calculating the mid-month date in decimal form - grace_L2_input['time'] = start_yr + mid_day / dpy + SHM['time'] = start_yr + mid_day / dpy + # Calculating the Julian dates of the start and end date - grace_L2_input['start'] = ( - 2400000.5 - + gravity_toolkit.time.convert_calendar_dates( - start_yr, 1.0, start_day, epoch=(1858, 11, 17, 0, 0, 0) - ) + MJD1 = gravity_toolkit.time.convert_calendar_dates( + start_yr, + start_struct.tm_mon, + start_struct.tm_mday, + epoch=(1858, 11, 17, 0, 0, 0), ) - grace_L2_input['end'] = ( - 2400000.5 - + gravity_toolkit.time.convert_calendar_dates( - end_yr, 1.0, end_day, epoch=(1858, 11, 17, 0, 0, 0) - ) + MJD2 = gravity_toolkit.time.convert_calendar_dates( + end_yr, + end_struct.tm_mon, + end_struct.tm_mday, + epoch=(1858, 11, 17, 0, 0, 0), ) + SHM['start'] = 2400000.5 + MJD1 + SHM['end'] = 2400000.5 + MJD2 # set maximum spherical harmonic order - MMAX = ( - np.copy(LMAX) if (kwargs['MMAX'] is None) else np.copy(kwargs['MMAX']) - ) + MMAX = kwargs.get('MMAX', None) + # only replace if None (allow MMAX to be zero, which is typically falsy) + if MMAX is None: + MMAX = np.copy(LMAX) # output dimensions - grace_L2_input['l'] = np.arange(LMAX + 1) - grace_L2_input['m'] = np.arange(MMAX + 1) + SHM['l'] = np.arange(LMAX + 1) + SHM['m'] = np.arange(MMAX + 1) # Spherical harmonic coefficient matrices to be filled from data file - grace_L2_input['clm'] = np.zeros((LMAX + 1, MMAX + 1)) - grace_L2_input['slm'] = np.zeros((LMAX + 1, MMAX + 1)) + SHM['clm'] = np.zeros((LMAX + 1, MMAX + 1)) + SHM['slm'] = np.zeros((LMAX + 1, MMAX + 1)) # spherical harmonic uncalibrated standard deviations - grace_L2_input['eclm'] = np.zeros((LMAX + 1, MMAX + 1)) - grace_L2_input['eslm'] = np.zeros((LMAX + 1, MMAX + 1)) + SHM['eclm'] = np.zeros((LMAX + 1, MMAX + 1)) + SHM['eslm'] = np.zeros((LMAX + 1, MMAX + 1)) if (DREL == 4) and (DSET == 'GSM'): # clm and slm drift rates for RL04 drift_c = np.zeros((LMAX + 1, MMAX + 1)) drift_s = np.zeros((LMAX + 1, MMAX + 1)) # set default degree 0 harmonics for intercomparability between centers - grace_L2_input['clm'][0, 0] = 1.0 + SHM['clm'][0, 0] = 1.0 # extract GRACE and GRACE-FO file headers # replace colons in header if within quotations @@ -223,15 +236,13 @@ def read_GRACE_harmonics(input_file, LMAX, **kwargs): ] header_regex = re.compile(r'(' + r'|'.join(header_parameters) + r')') header = [l.split(maxsplit=1) for l in head if header_regex.match(l)] - grace_L2_input['header'] = {i[0]: i[1] for i in header} + SHM['header'] = {i[0]: i[1] for i in header} elif ((N == 'GRAC') and (DREL >= 6)) or (N == 'GRFO'): # parse the YAML header for RL06 or GRACE-FO (specifying yaml loader) - grace_L2_input.update( - yaml.load('\n'.join(head), Loader=yaml.BaseLoader) - ) + SHM.update(yaml.load('\n'.join(head), Loader=yaml.BaseLoader)) else: # save lines of the GRACE file header removing empty lines - grace_L2_input['header'] = [l.rstrip() for l in head if l] + SHM['header'] = [l.rstrip() for l in head if l] # for each line in the GRACE/GRACE-FO file for line in file_contents: @@ -244,10 +255,10 @@ def read_GRACE_harmonics(input_file, LMAX, **kwargs): m1 = np.int64(line_contents[2]) # if degree and order are below the truncation limits if (l1 <= LMAX) and (m1 <= MMAX): - grace_L2_input['clm'][l1, m1] = np.float64(line_contents[3]) - grace_L2_input['slm'][l1, m1] = np.float64(line_contents[4]) - grace_L2_input['eclm'][l1, m1] = np.float64(line_contents[5]) - grace_L2_input['eslm'][l1, m1] = np.float64(line_contents[6]) + SHM['clm'][l1, m1] = np.float64(line_contents[3]) + SHM['slm'][l1, m1] = np.float64(line_contents[4]) + SHM['eclm'][l1, m1] = np.float64(line_contents[5]) + SHM['eslm'][l1, m1] = np.float64(line_contents[6]) # find if line starts with drift rate flag elif bool(re.match(r'GRDOTA', line)): # split the line into individual components @@ -264,21 +275,21 @@ def read_GRACE_harmonics(input_file, LMAX, **kwargs): # Currently removes 2003.3 to get the temporal average close to 0. if (DREL == 4) and (DSET == 'GSM'): # time since 2003.3 - dt = grace_L2_input['time'] - 2003.3 - grace_L2_input['clm'][:, :] += dt * drift_c[:, :] - grace_L2_input['slm'][:, :] += dt * drift_s[:, :] + dt = SHM['time'] - 2003.3 + SHM['clm'][:, :] += dt * drift_c[:, :] + SHM['slm'][:, :] += dt * drift_s[:, :] # Correct Pole Tide following Wahr et al. (2015) 10.1002/2015JB011986 if kwargs['POLE_TIDE'] and (DSET == 'GSM'): # time since 2000.0 - dt = grace_L2_input['time'] - 2000.0 + dt = SHM['time'] - 2000.0 # CSR and JPL Pole Tide Correction if PRC in ('UTCSR', 'JPLEM', 'JPLMSC'): # values for IERS mean pole [2010] - if grace_L2_input['time'] < 2010.0: + if SHM['time'] < 2010.0: a = np.array([0.055974, 1.8243e-3, 1.8413e-4, 7.024e-6]) b = np.array([-0.346346, -1.7896e-3, 1.0729e-4, 0.908e-6]) - elif grace_L2_input['time'] >= 2010.0: + elif SHM['time'] >= 2010.0: a = np.array([0.023513, 7.6141e-3, 0.0, 0.0]) b = np.array([-0.358891, 0.6287e-3, 0.0, 0.0]) # calculate m1 and m2 values @@ -297,8 +308,8 @@ def read_GRACE_harmonics(input_file, LMAX, **kwargs): m2 + 3.48e-3 * dt ) # correct GRACE/GRACE-FO spherical harmonics for pole tide - grace_L2_input['clm'][2, 1] -= C21_PT - grace_L2_input['slm'][2, 1] -= S21_PT + SHM['clm'][2, 1] -= C21_PT + SHM['slm'][2, 1] -= S21_PT # GFZ Pole Tide Correction elif PRC in ('EIGEN', 'GFZOP'): # pole tide values for GFZ @@ -306,13 +317,13 @@ def read_GRACE_harmonics(input_file, LMAX, **kwargs): C21_PT = -1.551e-9 * (-0.62e-3 * dt) - 0.012e-9 * (3.48e-3 * dt) S21_PT = 0.021e-9 * (-0.62e-3 * dt) - 1.505e-9 * (3.48e-3 * dt) # correct GRACE/GRACE-FO spherical harmonics for pole tide - grace_L2_input['clm'][2, 1] -= C21_PT - grace_L2_input['slm'][2, 1] -= S21_PT + SHM['clm'][2, 1] -= C21_PT + SHM['slm'][2, 1] -= S21_PT # return the header data, GRACE/GRACE-FO data # GRACE/GRACE-FO date (mid-month in decimal) # and the start and end days as Julian dates - return grace_L2_input + return SHM # PURPOSE: extract parameters from filename diff --git a/gravity_toolkit/read_gfc_harmonics.py b/gravity_toolkit/read_gfc_harmonics.py index 314a0da..bd21cd2 100644 --- a/gravity_toolkit/read_gfc_harmonics.py +++ b/gravity_toolkit/read_gfc_harmonics.py @@ -1,7 +1,7 @@ #!/usr/bin/env python """ read_gfc_harmonics.py -Written by Tyler Sutterley (06/2023) +Written by Tyler Sutterley (08/2026) Contributions by Hugo Lecomte Reads gfc files and extracts spherical harmonics for Swarm and @@ -55,6 +55,7 @@ calculate_tidal_offset.py: calculates the C20 offset for a tidal system UPDATE HISTORY: + Updated 08/2026: rename output dictionary to gfc for consistency Updated 06/2024: use wrapper to importlib for optional dependencies Updated 05/2023: use pathlib to define and operate on paths Updated 03/2023: improve typing for variables in docstrings @@ -74,6 +75,7 @@ import numpy as np import gravity_toolkit.time from gravity_toolkit.utilities import import_dependency +from datetime import datetime # attempt imports geoidtk = import_dependency('geoid_toolkit') @@ -165,10 +167,12 @@ def read_gfc_harmonics(input_file, TIDE=None, FLAG='gfc'): # extract parameters from input filename PFX, PRD, trunc, year, month, SFX = rx.findall(input_file.name).pop() # number of days in each month for the calendar year - dpm = gravity_toolkit.time.calendar_days(int(year)) + dpm = gravity_toolkit.time.calendar_days(int(year), astype=int) # create start and end date lists - start_date = [int(year), int(month), 1, 0, 0, 0] - end_date = [int(year), int(month), dpm[int(month) - 1], 23, 59, 59] + start_date = datetime(int(year), int(month), 1, 0, 0, 0) + end_date = datetime( + int(year), int(month), dpm[int(month) - 1], 23, 59, 59 + ) elif re.match(swarm_data, input_file.name): # compile numerical expression operator for parameters from files # Swarm: data from Swarm satellite @@ -177,10 +181,10 @@ def read_gfc_harmonics(input_file, TIDE=None, FLAG='gfc'): SAT, tmp, PROD, starttime, endtime, RL, SFX = rx.findall( input_file.name ).pop() - start_date, _ = gravity_toolkit.time.parse_date_string(starttime) - end_date, _ = gravity_toolkit.time.parse_date_string(endtime) + start_date = gravity_toolkit.time.parse(starttime) + end_date = gravity_toolkit.time.parse(endtime) # number of days in each month for the calendar year - dpm = gravity_toolkit.time.calendar_days(start_date[0]) + dpm = gravity_toolkit.time.calendar_days(start_date.year, astype=int) elif re.match(swarm_model, input_file.name): # compile numerical expression operator for parameters from files # Swarm: dealiasing products for Swarm data @@ -188,63 +192,56 @@ def read_gfc_harmonics(input_file, TIDE=None, FLAG='gfc'): # extract parameters from input filename PROD, trunc, month, year, SFX = rx.findall(input_file.name).pop() # number of days in each month for the calendar year - dpm = gravity_toolkit.time.calendar_days(int(year)) + dpm = gravity_toolkit.time.calendar_days(int(year), astype=int) # create start and end date lists - start_date = [int(year), int(month), 1, 0, 0, 0] - end_date = [int(year), int(month), dpm[int(month) - 1], 23, 59, 59] + start_date = datetime(int(year), int(month), 1, 0, 0, 0) + end_date = datetime( + int(year), int(month), dpm[int(month) - 1], 23, 59, 59 + ) # python dictionary with model input and headers ZIP = bool(re.search('ZIP', SFX, re.IGNORECASE)) - model_input = geoidtk.read_ICGEM_harmonics( + # read gravity field coefficients (gfc) + gfc = geoidtk.read_ICGEM_harmonics( input_file, TIDE=TIDE, FLAG=FLAG, ZIP=ZIP ) # start and end day of the year - start_day = ( - np.sum(dpm[: start_date[1] - 1]) - + start_date[2] - + start_date[3] / 24.0 - + start_date[4] / 1440.0 - + start_date[5] / 86400.0 - ) - end_day = ( - np.sum(dpm[: end_date[1] - 1]) - + end_date[2] - + end_date[3] / 24.0 - + end_date[4] / 1440.0 - + end_date[5] / 86400.0 - ) + start_struct = start_date.timetuple() + start_yr = start_struct.tm_year + start_day = start_struct.tm_yday + end_struct = end_date.timetuple() + end_yr = end_struct.tm_year + end_day = end_struct.tm_yday + # end date taking into account measurements taken on different years - end_cyclic = (end_date[0] - start_date[0]) * np.sum(dpm) + end_day + end_cyclic = (end_yr - start_yr) * np.sum(dpm) + end_day # calculate mid-month value mid_day = np.mean([start_day, end_cyclic]) # Calculating the mid-month date in decimal form - model_input['time'] = start_date[0] + mid_day / np.sum(dpm) + gfc['time'] = start_yr + mid_day / np.sum(dpm) + # Calculating the Julian dates of the start and end date - model_input['start'] = ( - 2400000.5 - + gravity_toolkit.time.convert_calendar_dates( - start_date[0], - start_date[1], - start_date[2], - hour=start_date[3], - minute=start_date[4], - second=start_date[5], - epoch=(1858, 11, 17, 0, 0, 0), - ) + MJD1 = gravity_toolkit.time.convert_calendar_dates( + start_yr, + start_struct.tm_mon, + start_struct.tm_mday, + hour=start_struct.tm_hour, + minute=start_struct.tm_min, + second=start_struct.tm_sec, + epoch=(1858, 11, 17, 0, 0, 0), ) - model_input['end'] = ( - 2400000.5 - + gravity_toolkit.time.convert_calendar_dates( - end_date[0], - end_date[1], - end_date[2], - hour=end_date[3], - minute=end_date[4], - second=end_date[5], - epoch=(1858, 11, 17, 0, 0, 0), - ) + MJD2 = gravity_toolkit.time.convert_calendar_dates( + end_yr, + end_struct.tm_mon, + end_struct.tm_mday, + hour=end_struct.tm_hour, + minute=end_struct.tm_min, + second=end_struct.tm_sec, + epoch=(1858, 11, 17, 0, 0, 0), ) + gfc['start'] = 2400000.5 + MJD1 + gfc['end'] = 2400000.5 + MJD2 # return the spherical harmonics and parameters - return model_input + return gfc diff --git a/gravity_toolkit/scripts/calc_degree_one.py b/gravity_toolkit/scripts/calc_degree_one.py index 4dd0794..df98628 100755 --- a/gravity_toolkit/scripts/calc_degree_one.py +++ b/gravity_toolkit/scripts/calc_degree_one.py @@ -883,8 +883,8 @@ def calc_degree_one( 'lmh...,lm...->mh...', plmout[l2, :, :], Ylms.ilm[l2, :] ) # Multiplying by c/s(phi#m) to get surface density in cmwe (lon,lat) - # ccos/ssin are mXphi, pcos/psin are mXtheta: resultant matrices are phiXtheta - # The summation over spherical harmonic order is in this multiplication + # resultant matrices are phiXtheta + # (summation over spherical harmonic order is in this multiplication) rmass = np.einsum('mp...,mh...->ph...', m_phi, pconv).real # calculate G matrix parameters through a summation of each latitude # summation of integration factors, Legendre polynomials, @@ -932,8 +932,8 @@ def calc_degree_one( ) # Multiplying by c/s(phi#m) to get surface density in cm w.e. (lonxlat) - # ccos/ssin are mXphi, pcos/psin are mXtheta: resultant matrices are phiXtheta - # The summation over spherical harmonic order is in this multiplication + # resultant matrices are phiXtheta + # (summation over spherical harmonic order is in this multiplication) lmass = np.einsum('mp...,mh...->ph...', m_phi, pconv).real # use sea level fingerprints or eustatic from GRACE land components diff --git a/gravity_toolkit/scripts/grace_raster_grids.py b/gravity_toolkit/scripts/grace_raster_grids.py index 08caa6f..fb3c079 100644 --- a/gravity_toolkit/scripts/grace_raster_grids.py +++ b/gravity_toolkit/scripts/grace_raster_grids.py @@ -38,9 +38,6 @@ CF: Center of Surface Figure (default) CM: Center of Mass of Earth System CE: Center of Mass of Solid Earth - -F X, --format X: input/output data format - netCDF4 - HDF5 -G X, --gia X: GIA model type to read IJ05-R2: Ivins R2 GIA Models W12a: Whitehouse GIA Models @@ -148,6 +145,9 @@ UPDATE HISTORY: Updated 08/2026: use default file logger for valid and failed program runs + use new structured netCDF4 output function from geoid-toolkit + reorder dimensions to be time, y, x for output netCDF4 files + include additional attributes to output netCDF4 files for CF compliance Updated 06/2024: use wrapper to importlib for optional dependencies Updated 03/2024: increase mask buffer to twice the smoothing radius Written 08/2023 @@ -237,7 +237,6 @@ def grace_raster_grids( SLR_C30=None, SLR_C40=None, SLR_C50=None, - DATAFORM=None, MEAN_FILE=None, MEANFORM=None, REMOVE_FILES=None, @@ -257,20 +256,22 @@ def grace_raster_grids( # output attributes for raster files attributes = dict(ROOT=collections.OrderedDict()) + # add attributes for software information + attributes['ROOT']['software_reference'] = gravtk.version.project_name + attributes['ROOT']['software_version'] = gravtk.version.full_version + # add attributes for GRACE/GRACE-FO information attributes['ROOT']['generating_institute'] = PROC attributes['ROOT']['product_release'] = DREL attributes['ROOT']['product_name'] = DSET attributes['ROOT']['product_type'] = 'gravity_field' attributes['ROOT']['title'] = 'GRACE/GRACE-FO Spatial Data' - attributes['ROOT']['reference'] = ( - f'Output from {pathlib.Path(sys.argv[0]).name}' - ) + # add citation to John's 1998 paper + attributes['ROOT']['citation'] = 'https://doi.org/10.1029/98jb02844' + reference = f'Output from {pathlib.Path(sys.argv[0]).name}' + attributes['ROOT']['reference'] = reference # list object of output files for file logs (full path) output_files = [] - # file information - suffix = dict(netCDF4='nc', HDF5='H5')[DATAFORM] - # read arrays of kl, hl, and ll Love Numbers LOVE = gravtk.load_love_numbers( LMAX, LOVE_NUMBERS=LOVE_NUMBERS, REFERENCE=REFERENCE, FORMAT='class' @@ -327,6 +328,9 @@ def grace_raster_grids( # add attributes for input GRACE/GRACE-FO spherical harmonics for att_name, att_val in Ylms['attributes'].items(): attributes['ROOT'][att_name] = att_val + # get the start and end dates for the GRACE/GRACE-FO data + SD = Ylms.attrs['start_date'].min().astype('datetime64[D]') + ED = Ylms.attrs['end_date'].max().astype('datetime64[D]') # use a mean file for the static field to remove if MEAN_FILE: @@ -435,10 +439,12 @@ def grace_raster_grids( crs1 = get_projection(PROJECTION) crs2 = pyproj.CRS.from_epsg(4326) transformer = pyproj.Transformer.from_crs(crs1, crs2, always_xy=True) - # dictionary of coordinate reference system variables + # dictionaries of coordinate reference system variables crs_to_dict = crs1.to_dict() + crs_to_cf = crs1.to_cf() + standard_name = crs_to_cf['grid_mapping_name'].title() # Climate and Forecast (CF) Metadata Conventions - if crs1.to_epsg() == 4326: + if crs1.to_epsg() == 4326 or crs1.is_geographic: y_cf, x_cf = crs1.cs_to_cf() else: x_cf, y_cf = crs1.cs_to_cf() @@ -460,62 +466,84 @@ def grace_raster_grids( attributes['ROOT']['earth_density'] = f'{factors.rho_e:0.3f} g/cm^3' attributes['ROOT']['earth_gravity_constant'] = f'{factors.GM:0.3f} cm^3/s^2' + # dictionary describing the output netCDF4 structure + struct = dict( + dimensions=('time', 'y', 'x'), + variables={ + 'z': ('time', 'y', 'x'), + 'crs': (), + }, + ) + # projection attributes attributes['crs'] = {} - # add projection attributes - attributes['crs']['standard_name'] = crs1.to_cf()[ - 'grid_mapping_name' - ].title() + attributes['crs']['standard_name'] = standard_name attributes['crs']['spatial_epsg'] = crs1.to_epsg() attributes['crs']['spatial_ref'] = crs1.to_wkt() attributes['crs']['proj4_params'] = crs1.to_proj4() - for att_name, att_val in crs1.to_cf().items(): + for att_name, att_val in crs_to_cf.items(): attributes['crs'][att_name] = att_val - if 'lat_0' in crs_to_dict.keys() and (crs1.to_epsg() != 4326): - attributes['crs']['latitude_of_projection_origin'] = crs_to_dict[ - 'lat_0' - ] - # x and y + if 'lat_0' in crs_to_dict.keys() and not crs1.is_geographic: + lat_0 = crs_to_dict.get('lat_0', None) + attributes['crs']['latitude_of_projection_origin'] = lat_0 + # x and y coordinate attributes attributes['x'], attributes['y'] = ({}, {}) for att_name in ['long_name', 'standard_name', 'units']: attributes['x'][att_name] = x_cf[att_name] attributes['y'][att_name] = y_cf[att_name] - # time + # time attributes attributes['time'] = {} attributes['time']['units'] = 'years' attributes['time']['long_name'] = 'Date_in_Decimal_Years' - # output gridded data - fill_value = -9999.0 + attributes['time']['standard_name'] = 'time' + attributes['time']['calendar'] = 'standard' + # data attributes attributes['z'] = {} attributes['z']['units'] = units_name attributes['z']['long_name'] = units_longname + attributes['z']['short_name'] = units attributes['z']['degree_of_truncation'] = LMAX - attributes['z']['_FillValue'] = fill_value - # set grid mapping attribute attributes['z']['grid_mapping'] = 'crs' + attributes['z']['coordinates'] = ' '.join(struct['dimensions']) # output data variables output = {} # projection variable - output['crs'] = np.array((), dtype=np.byte) + output['crs'] = np.byte() # spacing and bounds of output grid dx, dy = np.broadcast_to(np.atleast_1d(SPACING), (2,)) xmin, xmax, ymin, ymax = np.copy(BOUNDS) # create x and y from spacing and bounds output['x'] = np.arange(xmin + dx / 2.0, xmax + dx, dx) - output['y'] = np.arange(ymin + dx / 2.0, ymax + dy, dy) + output['y'] = np.arange(ymin + dy / 2.0, ymax + dy, dy) ny, nx = (len(output['y']), len(output['x'])) + # output time variable + output['time'] = np.zeros((nt)) + # output gridded raster data + fill_value = -9999.0 + output['z'] = np.ma.zeros((nt, ny, nx), fill_value=fill_value) + output['z'].mask = np.ones((nt, ny, nx), dtype=bool) + + # create meshgrid of x and y gridx, gridy = np.meshgrid(output['x'], output['y']) gridlon, gridlat = transformer.transform(gridx, gridy) - - # semimajor axis of ellipsoid [m] - a_axis = crs1.ellipsoid.semi_major_metre # ellipsoidal flattening flat = 1.0 / crs1.ellipsoid.inverse_flattening # calculate geocentric latitude and convert to degrees latitude_geocentric = geoidtk.spatial.geocentric_latitude( - gridlon, gridlat, a_axis=a_axis, flat=flat + gridlat, flat=flat ) + # add geospatial attributes + attributes['ROOT']['geospatial_lat_min'] = gridlat.min() + attributes['ROOT']['geospatial_lat_max'] = gridlat.max() + attributes['ROOT']['geospatial_lon_min'] = gridlon.min() + attributes['ROOT']['geospatial_lon_max'] = gridlon.max() + attributes['ROOT']['geospatial_lat_units'] = 'degrees_north' + attributes['ROOT']['geospatial_lon_units'] = 'degrees_east' + # add temporal attributes + attributes['ROOT']['time_coverage_start'] = np.datetime_as_string(SD) + attributes['ROOT']['time_coverage_end'] = np.datetime_as_string(ED) + attributes['ROOT']['time_coverage_duration'] = str(ED - SD) # calculate spatial mask with an extended radius THRESHOLD = 0.025 @@ -531,26 +559,22 @@ def grace_raster_grids( ) ii, jj = np.nonzero(mask.reshape(ny, nx) > THRESHOLD) - # output gridded raster data - output['z'] = np.ma.zeros((ny, nx, nt), fill_value=fill_value) - output['z'].mask = np.ones((ny, nx, nt), dtype=bool) - output['time'] = np.zeros((nt)) - # converting harmonics to truncated, smoothed coefficients in units # combining harmonics to calculate output raster grids - for i, grace_month in enumerate(GRACE_Ylms.month): + for t, grace_month in enumerate(GRACE_Ylms.month): + # keep track of the current month being processed logger.debug(grace_month) # GRACE/GRACE-FO harmonics for time t - Ylms = GRACE_Ylms.index(i) + Ylms = GRACE_Ylms.index(t) # Remove GIA rate for time - Ylms.subtract(GIA_Ylms.index(i)) + Ylms.subtract(GIA_Ylms.index(t)) # Remove monthly files to be removed - Ylms.subtract(remove_Ylms.index(i)) + Ylms.subtract(remove_Ylms.index(t)) # truncate to degree and order LMAX and MMAX # truncate minimum degree to LMIN Ylms.truncate(LMAX, lmin=LMIN, mmax=MMAX) # convert spherical harmonics to output raster grid - output['z'].data[ii, jj, i] = gravtk.clenshaw_summation( + output['z'].data[t, ii, jj] = gravtk.clenshaw_summation( Ylms.clm, Ylms.slm, gridlon[ii, jj], @@ -560,27 +584,25 @@ def grace_raster_grids( LMAX=LMAX, LOVE=LOVE, ) - output['z'].mask[ii, jj, i] = False + output['z'].mask[t, ii, jj] = False # copy time variables for month - output['time'][i] = np.copy(Ylms.time) + output['time'][t] = np.copy(Ylms.time) # convert masked values to fill value output['z'].data[output['z'].mask] = output['z'].fill_value - # output raster files to netCDF4 or HDF5 - FILE = ( + # output netCDF4 raster files + output_file = OUTPUT_DIRECTORY.joinpath( f'{FILE_PREFIX}{units}_L{LMAX:d}{order_str}{gw_str}{ds_str}_' - f'{START:03d}-{END:03d}.{suffix}' + f'{START:03d}-{END:03d}.nc' ) - output_file = OUTPUT_DIRECTORY.joinpath(FILE) # use spatial functions from geoid toolkit to write rasters - if DATAFORM == 'netCDF4': - geoidtk.spatial.to_netCDF4( - output, attributes, output_file, data_type='grid' - ) - elif DATAFORM == 'HDF5': - geoidtk.spatial.to_HDF5( - output, attributes, output_file, data_type='grid' - ) + geoidtk.spatial.to_netCDF4( + output, + attributes, + output_file, + data_type='structured', + structure=struct, + ) # set the permissions mode of the output files output_file.chmod(mode=MODE) # add file to list @@ -907,15 +929,6 @@ def arguments(): choices=['CSR', 'GSFC', 'LARES'], help='Replace C50 coefficients with SLR values', ) - # input data format (netCDF4, HDF5) - parser.add_argument( - '--format', - '-F', - type=str, - default='netCDF4', - choices=['netCDF4', 'HDF5'], - help='Input/output data format', - ) # mean file to remove parser.add_argument( '--mean-file', @@ -1038,7 +1051,6 @@ def main(): SLR_C30=args.slr_c30, SLR_C40=args.slr_c40, SLR_C50=args.slr_c50, - DATAFORM=args.format, MEAN_FILE=args.mean_file, MEANFORM=args.mean_format, REMOVE_FILES=args.remove_file, diff --git a/gravity_toolkit/scripts/grace_spatial_error.py b/gravity_toolkit/scripts/grace_spatial_error.py index ab44657..573d1be 100755 --- a/gravity_toolkit/scripts/grace_spatial_error.py +++ b/gravity_toolkit/scripts/grace_spatial_error.py @@ -120,6 +120,8 @@ UPDATE HISTORY: Updated 08/2026: use default file logger for valid and failed program runs + include additional attributes to output files for CF compliance + use np.einsum and Euler's formula for spherical harmonic summations Updated 05/2023: use pathlib to define and operate on paths Updated 03/2023: add root attributes to output netCDF4 and HDF5 files use attributes from units class for writing to netCDF4/HDF5 files @@ -243,6 +245,8 @@ def grace_spatial_error( attributes['product_name'] = DSET attributes['product_type'] = 'gravity_field' attributes['title'] = 'GRACE/GRACE-FO Spatial Error' + # add citation to John's 2006 paper + attributes['citation'] = 'http://doi.org/10.1029/2005GL025305' # list object of output files for file logs (full path) output_files = [] @@ -308,6 +312,9 @@ def grace_spatial_error( # add attributes for input GRACE/GRACE-FO spherical harmonics for att_name, att_val in Ylms['attributes'].items(): attributes[att_name] = att_val + # get the start and end dates for the GRACE/GRACE-FO data + SD = Ylms.attrs['start_date'].min().astype('datetime64[D]') + ED = Ylms.attrs['end_date'].max().astype('datetime64[D]') # use a mean file for the static field to remove if MEAN_FILE: @@ -416,7 +423,7 @@ def grace_spatial_error( nsmth = np.int64(delta_Ylms.month) # Output spatial data object - delta = gravtk.spatial() + grid = gravtk.spatial() # Output Degree Spacing dlon, dlat = (DDEG[0], DDEG[0]) if (len(DDEG) == 1) else (DDEG[0], DDEG[1]) # Output Degree Interval @@ -424,21 +431,21 @@ def grace_spatial_error( # (-180:180,90:-90) nlon = np.int64((360.0 / dlon) + 1.0) nlat = np.int64((180.0 / dlat) + 1.0) - delta.lon = -180 + dlon * np.arange(0, nlon) - delta.lat = 90.0 - dlat * np.arange(0, nlat) + grid.lon = -180 + dlon * np.arange(0, nlon) + grid.lat = 90.0 - dlat * np.arange(0, nlat) elif INTERVAL == 2: # (Degree spacing)/2 - delta.lon = np.arange(-180 + dlon / 2.0, 180 + dlon / 2.0, dlon) - delta.lat = np.arange(90.0 - dlat / 2.0, -90.0 - dlat / 2.0, -dlat) - nlon = len(delta.lon) - nlat = len(delta.lat) + grid.lon = np.arange(-180 + dlon / 2.0, 180 + dlon / 2.0, dlon) + grid.lat = np.arange(90.0 - dlat / 2.0, -90.0 - dlat / 2.0, -dlat) + nlon = len(grid.lon) + nlat = len(grid.lat) elif INTERVAL == 3: # non-global grid set with BOUNDS parameter minlon, maxlon, minlat, maxlat = BOUNDS.copy() - delta.lon = np.arange(minlon + dlon / 2.0, maxlon + dlon / 2.0, dlon) - delta.lat = np.arange(maxlat - dlat / 2.0, minlat - dlat / 2.0, -dlat) - nlon = len(delta.lon) - nlat = len(delta.lat) + grid.lon = np.arange(minlon + dlon / 2.0, maxlon + dlon / 2.0, dlon) + grid.lat = np.arange(maxlat - dlat / 2.0, minlat - dlat / 2.0, -dlat) + nlon = len(grid.lon) + nlat = len(grid.lat) # output spatial units # dfactor is the degree dependent coefficients @@ -457,22 +464,32 @@ def grace_spatial_error( attributes['earth_radius'] = f'{factors.rad_e:0.3f} cm' attributes['earth_density'] = f'{factors.rho_e:0.3f} g/cm^3' attributes['earth_gravity_constant'] = f'{factors.GM:0.3f} cm^3/s^2' + # add geospatial attributes + attributes['geospatial_lat_min'] = grid.lat.min() + attributes['geospatial_lat_max'] = grid.lat.max() + attributes['geospatial_lon_min'] = grid.lon.min() + attributes['geospatial_lon_max'] = grid.lon.max() + attributes['geospatial_lat_units'] = 'degrees_north' + attributes['geospatial_lon_units'] = 'degrees_east' + # add temporal attributes + attributes['time_coverage_start'] = np.datetime_as_string(SD) + attributes['time_coverage_end'] = np.datetime_as_string(ED) + attributes['time_coverage_duration'] = str(ED - SD) # add attributes to output spatial object attributes['reference'] = f'Output from {pathlib.Path(sys.argv[0]).name}' - delta.attributes['ROOT'] = attributes + grid.attributes['ROOT'] = attributes # Computing plms for converting to spatial domain - phi = np.radians(delta.lon[np.newaxis, :]) - theta = np.radians(90.0 - delta.lat) + phi = np.radians(grid.lon) + theta = np.radians(90.0 - grid.lat) PLM, dPLM = gravtk.plm_holmes(LMAX, np.cos(theta)) # square of legendre polynomials truncated to order MMAX mm = np.arange(0, MMAX + 1) PLM2 = PLM[:, mm, :] ** 2 # Calculating cos(m*phi)^2 and sin(m*phi)^2 - m = delta_Ylms.m[:, np.newaxis] - ccos = np.cos(np.dot(m, phi)) ** 2 - ssin = np.sin(np.dot(m, phi)) ** 2 + mp = np.einsum('m...,p...->mp...', mm, phi) + m_phi2 = (1.0 - 1j) * (np.exp(-2j * mp) + np.exp(2j * mp) + 2j) / 4.0 # truncate delta harmonics to spherical harmonic range Ylms = delta_Ylms.truncate(LMAX, lmin=LMIN, mmax=MMAX) @@ -480,20 +497,16 @@ def grace_spatial_error( # smooth harmonics and convert to output units Ylms = Ylms.convolve(dfactor * wt).power(2.0).scale(1.0 / nsmth) # Calculate fourier coefficients - d_cos = np.zeros((MMAX + 1, nlat)) # [m,th] - d_sin = np.zeros((MMAX + 1, nlat)) # [m,th] - # Calculating delta spatial values - for k in range(0, nlat): - # summation over all spherical harmonic degrees - d_cos[:, k] = np.sum(PLM2[:, :, k] * Ylms.clm, axis=0) - d_sin[:, k] = np.sum(PLM2[:, :, k] * Ylms.slm, axis=0) + # summation over all spherical harmonic degrees + pconv2 = np.einsum('lmh...,lm...->mh...', PLM2, Ylms.ilm) - # Multiplying by c/s(phi#m) to get spatial maps (lon,lat) - delta.data = np.sqrt(np.dot(ccos.T, d_cos) + np.dot(ssin.T, d_sin)).T + # Multiplying by c/s(phi#m) to get spatial error map + # take the square root and drop imaginary component + grid.data = np.sqrt(np.einsum('mp...,mh...->hp...', m_phi2, pconv2)).real # output file format file_format = '{0}{1}_L{2:d}{3}{4}{5}_ERR_{6:03d}-{7:03d}.{8}' - # output error file to ascii, netCDF4 or HDF5 + # build output filename fargs = ( FILE_PREFIX, units, @@ -505,8 +518,10 @@ def grace_spatial_error( GRACE_Ylms.month[-1], suffix[DATAFORM], ) - OUTPUT_FILE = OUTPUT_DIRECTORY.joinpath(file_format.format(*fargs)) - delta.to_file( + filename = file_format.format(*fargs) + OUTPUT_FILE = OUTPUT_DIRECTORY.joinpath(filename) + # write spatial data to output file + grid.to_file( OUTPUT_FILE, format=DATAFORM, date=False, diff --git a/gravity_toolkit/scripts/grace_spatial_maps.py b/gravity_toolkit/scripts/grace_spatial_maps.py index ccb1d59..5801c05 100755 --- a/gravity_toolkit/scripts/grace_spatial_maps.py +++ b/gravity_toolkit/scripts/grace_spatial_maps.py @@ -151,6 +151,7 @@ UPDATE HISTORY: Updated 08/2026: use default file logger for valid and failed program runs + include additional attributes to output files for CF compliance Updated 05/2023: use pathlib to define and operate on paths Updated 03/2023: add root attributes to output netCDF4 and HDF5 files use attributes from units class for writing to netCDF4/HDF5 files @@ -270,6 +271,8 @@ def grace_spatial_maps( attributes['product_name'] = DSET attributes['product_type'] = 'gravity_field' attributes['title'] = 'GRACE/GRACE-FO Spatial Data' + # add citation to John's 1998 paper + attributes['citation'] = 'https://doi.org/10.1029/98jb02844' # list object of output files for file logs (full path) output_files = [] @@ -333,6 +336,10 @@ def grace_spatial_maps( # add attributes for input GRACE/GRACE-FO spherical harmonics for att_name, att_val in Ylms['attributes'].items(): attributes[att_name] = att_val + # get the start and end dates for the GRACE/GRACE-FO data + SD = Ylms.attrs['start_date'].min().astype('datetime64[D]') + ED = Ylms.attrs['end_date'].max().astype('datetime64[D]') + # use a mean file for the static field to remove if MEAN_FILE: # read data form for input mean file (ascii, netCDF4, HDF5, gfc) @@ -478,6 +485,17 @@ def grace_spatial_maps( attributes['earth_radius'] = f'{factors.rad_e:0.3f} cm' attributes['earth_density'] = f'{factors.rho_e:0.3f} g/cm^3' attributes['earth_gravity_constant'] = f'{factors.GM:0.3f} cm^3/s^2' + # add geospatial attributes + attributes['ROOT']['geospatial_lat_min'] = grid.lat.min() + attributes['ROOT']['geospatial_lat_max'] = grid.lat.max() + attributes['ROOT']['geospatial_lon_min'] = grid.lon.min() + attributes['ROOT']['geospatial_lon_max'] = grid.lon.max() + attributes['ROOT']['geospatial_lat_units'] = 'degrees_north' + attributes['ROOT']['geospatial_lon_units'] = 'degrees_east' + # add temporal attributes + attributes['ROOT']['time_coverage_start'] = np.datetime_as_string(SD) + attributes['ROOT']['time_coverage_end'] = np.datetime_as_string(ED) + attributes['ROOT']['time_coverage_duration'] = str(ED - SD) # add attributes to output spatial object attributes['reference'] = f'Output from {pathlib.Path(sys.argv[0]).name}' grid.attributes['ROOT'] = attributes @@ -486,13 +504,13 @@ def grace_spatial_maps( file_format = '{0}{1}_L{2:d}{3}{4}{5}_{6:03d}.{7}' # converting harmonics to truncated, smoothed coefficients in units # combining harmonics to calculate output spatial fields - for i, grace_month in enumerate(GRACE_Ylms.month): + for t, grace_month in enumerate(GRACE_Ylms.month): # GRACE/GRACE-FO harmonics for time t - Ylms = GRACE_Ylms.index(i) + Ylms = GRACE_Ylms.index(t) # Remove GIA rate for time - Ylms.subtract(GIA_Ylms.index(i)) + Ylms.subtract(GIA_Ylms.index(t)) # Remove monthly files to be removed - Ylms.subtract(remove_Ylms.index(i)) + Ylms.subtract(remove_Ylms.index(t)) # smooth harmonics and convert to output units Ylms.convolve(dfactor * wt) # convert spherical harmonics to output spatial grid @@ -511,7 +529,7 @@ def grace_spatial_maps( grid.time = np.copy(Ylms.time) grid.month = np.copy(Ylms.month) - # output monthly files to ascii, netCDF4 or HDF5 + # build output filename fargs = ( FILE_PREFIX, units, @@ -522,7 +540,9 @@ def grace_spatial_maps( grace_month, suffix[DATAFORM], ) - OUTPUT_FILE = OUTPUT_DIRECTORY.joinpath(file_format.format(*fargs)) + filename = file_format.format(*fargs) + OUTPUT_FILE = OUTPUT_DIRECTORY.joinpath(filename) + # write spatial data to output file grid.to_file( OUTPUT_FILE, format=DATAFORM, diff --git a/gravity_toolkit/scripts/monte_carlo_degree_one.py b/gravity_toolkit/scripts/monte_carlo_degree_one.py index b32f720..7ea8aac 100644 --- a/gravity_toolkit/scripts/monte_carlo_degree_one.py +++ b/gravity_toolkit/scripts/monte_carlo_degree_one.py @@ -800,8 +800,8 @@ def monte_carlo_degree_one( 'lmh...,lm...->mh...', plmout[l2, :, :], GRACE_Ylms.ilm[l2, :] ) # Multiplying by c/s(phi#m) to get surface density in cmwe (lon,lat) - # ccos/ssin are mXphi, pcos/psin are mXtheta: resultant matrices are phiXtheta - # The summation over spherical harmonic order is in this multiplication + # resultant matrices are phiXtheta + # (summation over spherical harmonic order is in this multiplication) rmass = np.einsum('mp...,mh...->ph...', m_phi, pconv).real # calculate G matrix parameters through a summation of each latitude # summation of integration factors, Legendre polynomials, @@ -849,8 +849,8 @@ def monte_carlo_degree_one( ) # Multiplying by c/s(phi#m) to get surface density in cm w.e. (lonxlat) - # ccos/ssin are mXphi, pcos/psin are mXtheta: resultant matrices are phiXtheta - # The summation over spherical harmonic order is in this multiplication + # resultant matrices are phiXtheta + # (summation over spherical harmonic order is in this multiplication) lmass = np.einsum('mp...,mh...->ph...', m_phi, pconv).real # use sea level fingerprints or eustatic from GRACE land components diff --git a/gravity_toolkit/scripts/scale_grace_maps.py b/gravity_toolkit/scripts/scale_grace_maps.py index acc402f..9e6fa20 100644 --- a/gravity_toolkit/scripts/scale_grace_maps.py +++ b/gravity_toolkit/scripts/scale_grace_maps.py @@ -156,6 +156,8 @@ UPDATE HISTORY: Updated 08/2026: use default file logger for valid and failed program runs + use np.einsum and Euler's formula for spherical harmonic summations + include additional attributes to output files for CF compliance Updated 05/2023: use pathlib to define and operate on paths Updated 03/2023: use new scaling_factors inheritance of spatial class single input file with scaling factor variables @@ -193,8 +195,9 @@ import logging import pathlib import argparse -import numpy as np import traceback +import collections +import numpy as np import gravity_toolkit as gravtk @@ -261,6 +264,14 @@ def scale_grace_maps( if not OUTPUT_DIRECTORY.exists(): OUTPUT_DIRECTORY.mkdir(mode=MODE, parents=True, exist_ok=True) + # output attributes for spatial files + attributes = collections.OrderedDict() + attributes['generating_institute'] = PROC + attributes['product_release'] = DREL + attributes['product_name'] = DSET + attributes['product_type'] = 'gravity_field' + # add citation to Felix and Sean's 2012 paper + attributes['citation'] = 'https://doi.org/10.1029/2011WR011453' # list object of output files for file logs (full path) output_files = [] @@ -273,15 +284,33 @@ def scale_grace_maps( LOVE = gravtk.load_love_numbers( LMAX, LOVE_NUMBERS=LOVE_NUMBERS, REFERENCE=REFERENCE, FORMAT='class' ) + # add attributes for earth model and love numbers + attributes['earth_model'] = LOVE.model + attributes['earth_love_numbers'] = LOVE.citation + attributes['reference_frame'] = LOVE.reference # atmospheric ECMWF "jump" flag (if ATM) atm_str = '_wATM' if ATM else '' # output string for both LMAX==MMAX and LMAX != MMAX cases MMAX = np.copy(LMAX) if not MMAX else MMAX order_str = f'M{MMAX:d}' if (MMAX != LMAX) else '' + # add attributes for LMAX and MMAX + attributes['max_degree'] = LMAX + attributes['max_order'] = MMAX + # output spatial units units = 'cmwe' + # dfactor is the degree dependent coefficients + # for converting to centimeters water equivalent (cmwe) + factors = gravtk.units(lmax=LMAX).harmonic(*LOVE) + dfactor = factors.get(units) + # units attributes units_name, units_longname = gravtk.units.get_attributes(units) + # add attributes for earth parameters + attributes['earth_radius'] = f'{factors.rad_e:0.3f} cm' + attributes['earth_density'] = f'{factors.rho_e:0.3f} g/cm^3' + attributes['earth_gravity_constant'] = f'{factors.GM:0.3f} cm^3/s^2' + # invalid value fill_value = -9999.0 @@ -357,6 +386,13 @@ def scale_grace_maps( ) # create harmonics object from GRACE/GRACE-FO data GRACE_Ylms = gravtk.harmonics().from_dict(Ylms) + # add attributes for input GRACE/GRACE-FO spherical harmonics + for att_name, att_val in Ylms['attributes'].items(): + attributes[att_name] = att_val + # get the start and end dates for the GRACE/GRACE-FO data + SD = Ylms.attrs['start_date'].min().astype('datetime64[D]') + ED = Ylms.attrs['end_date'].max().astype('datetime64[D]') + # use a mean file for the static field to remove if MEAN_FILE: # read data form for input mean file (ascii, netCDF4, HDF5, gfc) @@ -527,18 +563,30 @@ def scale_grace_maps( grid.data = np.zeros((nlat, nlon, nfiles)) grid.mask = np.zeros((nlat, nlon, nfiles), dtype=bool) + # add geospatial attributes + attributes['geospatial_lat_min'] = grid.lat.min() + attributes['geospatial_lat_max'] = grid.lat.max() + attributes['geospatial_lon_min'] = grid.lon.min() + attributes['geospatial_lon_max'] = grid.lon.max() + attributes['geospatial_lat_units'] = 'degrees_north' + attributes['geospatial_lon_units'] = 'degrees_east' + # add temporal attributes + attributes['time_coverage_start'] = np.datetime_as_string(SD) + attributes['time_coverage_end'] = np.datetime_as_string(ED) + attributes['time_coverage_duration'] = str(ED - SD) + # add attributes to output spatial object + attributes['title'] = 'GRACE/GRACE-FO Scaled Spatial Data' + attributes['reference'] = f'Output from {pathlib.Path(sys.argv[0]).name}' + grid.attributes['ROOT'] = attributes + # Computing plms for converting to spatial domain - phi = np.radians(grid.lon[np.newaxis, :]) + phi = np.radians(grid.lon) theta = np.radians(90.0 - grid.lat) PLM, dPLM = gravtk.plm_holmes(LMAX, np.cos(theta)) # square of legendre polynomials truncated to order MMAX mm = np.arange(0, MMAX + 1) PLM2 = PLM[:, mm, :] ** 2 - # dfactor is the degree dependent coefficients - # for converting to centimeters water equivalent (cmwe) - dfactor = gravtk.units(lmax=LMAX).harmonic(*LOVE).cmwe - # converting harmonics to truncated, smoothed coefficients in units # combining harmonics to calculate output spatial fields for i, gm in enumerate(GRACE_Ylms.month): @@ -581,22 +629,30 @@ def scale_grace_maps( grid.month[-1], suffix[DATAFORM], ) - FILE = OUTPUT_DIRECTORY.joinpath(file_format.format(*fargs)) - # attributes for output files - attributes = {} - attributes['units'] = copy.copy(units_name) - attributes['longname'] = copy.copy(units_longname) - attributes['title'] = 'GRACE/GRACE-FO Spatial Data' - attributes['reference'] = f'Output from {pathlib.Path(sys.argv[0]).name}' + filename = file_format.format(*fargs) + FILE = OUTPUT_DIRECTORY.joinpath(filename) + # write to output file if DATAFORM == 'ascii': # ascii (.txt) grid.to_ascii(FILE, date=True, verbose=VERBOSE) elif DATAFORM == 'netCDF4': # netCDF4 - grid.to_netCDF4(FILE, date=True, verbose=VERBOSE, **attributes) + grid.to_netCDF4( + FILE, + date=True, + units=units_name, + longname=units_longname, + verbose=VERBOSE, + ) elif DATAFORM == 'HDF5': # HDF5 - grid.to_HDF5(FILE, date=True, verbose=VERBOSE, **attributes) + grid.to_HDF5( + FILE, + date=True, + units=units_name, + longname=units_longname, + verbose=VERBOSE, + ) # set the permissions mode of the output files FILE.chmod(mode=MODE) # add file to list @@ -613,6 +669,9 @@ def scale_grace_maps( error.data = kfactor.error * ratio error.mask = np.copy(kfactor.mask) error.update_mask() + # add attributes to output error object + attributes['title'] = 'GRACE/GRACE-FO Scaling Error' + error.attributes['ROOT'] = attributes # output monthly error files to ascii, netCDF4 or HDF5 fargs = ( @@ -627,18 +686,30 @@ def scale_grace_maps( grid.month[-1], suffix[DATAFORM], ) - FILE = OUTPUT_DIRECTORY.joinpath(file_format.format(*fargs)) - # attributes for output files - attributes['title'] = 'GRACE/GRACE-FO Scaling Error' + filename = file_format.format(*fargs) + FILE = OUTPUT_DIRECTORY.joinpath(filename) + # write to output file if DATAFORM == 'ascii': # ascii (.txt) error.to_ascii(FILE, date=False, verbose=VERBOSE) elif DATAFORM == 'netCDF4': # netCDF4 - error.to_netCDF4(FILE, date=False, verbose=VERBOSE, **attributes) + error.to_netCDF4( + FILE, + date=False, + units=units_name, + longname=units_longname, + verbose=VERBOSE, + ) elif DATAFORM == 'HDF5': # HDF5 - error.to_HDF5(FILE, date=False, verbose=VERBOSE, **attributes) + error.to_HDF5( + FILE, + date=False, + units=units_name, + longname=units_longname, + verbose=VERBOSE, + ) # set the permissions mode of the output files FILE.chmod(mode=MODE) # add file to list @@ -652,11 +723,11 @@ def scale_grace_maps( delta.month = np.copy(nsmth) delta.data = np.zeros((nlat, nlon)) delta.mask = np.zeros((nlat, nlon), dtype=bool) + # calculate scaled spatial error # Calculating cos(m*phi)^2 and sin(m*phi)^2 - m = delta_Ylms.m[:, np.newaxis] - ccos = np.cos(np.dot(m, phi)) ** 2 - ssin = np.sin(np.dot(m, phi)) ** 2 + mp = np.einsum('m...,p...->mp...', mm, phi) + m_phi2 = (1.0 - 1j) * (np.exp(-2j * mp) + np.exp(2j * mp) + 2j) / 4.0 # truncate delta harmonics to spherical harmonic range Ylms = delta_Ylms.truncate(LMAX, lmin=LMIN, mmax=MMAX) @@ -664,19 +735,19 @@ def scale_grace_maps( # smooth harmonics and convert to output units Ylms = Ylms.convolve(dfactor * wt).power(2.0).scale(1.0 / nsmth) # Calculate fourier coefficients - d_cos = np.zeros((MMAX + 1, nlat)) # [m,th] - d_sin = np.zeros((MMAX + 1, nlat)) # [m,th] - # Calculating delta spatial values - for k in range(0, nlat): - # summation over all spherical harmonic degrees - d_cos[:, k] = np.sum(PLM2[:, :, k] * Ylms.clm, axis=0) - d_sin[:, k] = np.sum(PLM2[:, :, k] * Ylms.slm, axis=0) + # summation over all spherical harmonic degrees + pconv2 = np.einsum('lmh...,lm...->mh...', PLM2, Ylms.ilm) + # Multiplying by c/s(phi#m) to get spatial error map - delta.data[:] = np.sqrt(np.dot(ccos.T, d_cos) + np.dot(ssin.T, d_sin)).T + # take the square root and drop imaginary component + delta.data = np.sqrt(np.einsum('mp...,mh...->hp...', m_phi2, pconv2)).real # scale output harmonic errors with kfactor delta = delta.scale(kfactor.data) delta.replace_invalid(fill_value, mask=kfactor.mask) + # add attributes to output error object + attributes['title'] = 'GRACE/GRACE-FO Scaled Spatial Error' + delta.attributes['ROOT'] = attributes # output monthly files to ascii, netCDF4 or HDF5 fargs = ( @@ -691,18 +762,30 @@ def scale_grace_maps( grid.month[-1], suffix[DATAFORM], ) - FILE = OUTPUT_DIRECTORY.joinpath(file_format.format(*fargs)) + filename = file_format.format(*fargs) + FILE = OUTPUT_DIRECTORY.joinpath(filename) # attributes for output files - attributes['title'] = 'GRACE/GRACE-FO Spatial Error' if DATAFORM == 'ascii': # ascii (.txt) delta.to_ascii(FILE, date=True, verbose=VERBOSE) elif DATAFORM == 'netCDF4': # netCDF4 - delta.to_netCDF4(FILE, date=True, verbose=VERBOSE, **attributes) + delta.to_netCDF4( + FILE, + date=True, + units=units_name, + longname=units_longname, + verbose=VERBOSE, + ) elif DATAFORM == 'HDF5': # HDF5 - delta.to_HDF5(FILE, date=True, verbose=VERBOSE, **attributes) + delta.to_HDF5( + FILE, + date=True, + units=units_name, + longname=units_longname, + verbose=VERBOSE, + ) # set the permissions mode of the output files FILE.chmod(mode=MODE) # add file to list diff --git a/gravity_toolkit/time.py b/gravity_toolkit/time.py index 8b89c9f..6eef200 100644 --- a/gravity_toolkit/time.py +++ b/gravity_toolkit/time.py @@ -1,7 +1,7 @@ #!/usr/bin/env python """ time.py -Written by Tyler Sutterley (07/2026) +Written by Tyler Sutterley (08/2026) Contributions by Hugo Lecomte Utilities for calculating time operations @@ -13,6 +13,8 @@ https://dateutil.readthedocs.io/en/stable/ UPDATE HISTORY: + Updated 08/2026: output numpy.datetime64 objects from file parsers + added astype option to calendar_days function (can now be int) Updated 07/2026: added function to convert times to numpy datetimes Updated 06/2024: remove timescale class and leap second calculations Updated 06/2023: add timescale class for converting between time scales @@ -46,11 +48,11 @@ import logging import pathlib import warnings -import datetime import traceback import numpy as np import dateutil.parser import gravity_toolkit.utilities +from datetime import datetime, timedelta # conversion factors between time units and seconds _to_sec = { @@ -228,7 +230,7 @@ def to_datetime(timedelta: np.ndarray, attributes: str, unit: str = 's'): # get the epoch and units from the attributes epoch, to_secs = parse_date_string(attributes) # convert epoch to datetime variable - epoch = np.datetime64(datetime.datetime(*epoch)) + epoch = np.datetime64(datetime(*epoch)) # calculate the delta time in seconds delta_time = np.atleast_1d(timedelta * to_secs).astype(np.int64) # return the datetime array @@ -266,7 +268,7 @@ def to_string( # return the string representations of the datetime objects if strftime is not None: # convert to datetime objects and use strftime formatting - dtime = dtime.astype(datetime.datetime) + dtime = dtime.astype(datetime) return np.array([d.strftime(strftime) for d in dtime]) else: return np.datetime_as_string(dtime, unit=unit) @@ -281,9 +283,16 @@ def parse_grace_file(granule): ---------- granule: str GRACE/GRACE-FO Level-2 spherical harmonic data file + + Returns + ------- + start_date: datetime + start date of GRACE/GRACE-FO data file + end_date: datetime + end date of GRACE/GRACE-FO data file """ # verify that filename is reduced to basename - file_basename = pathlib.Path(granule).name + name = pathlib.Path(granule).name # compile numerical expression operator for parameters from files # UTCSR: The University of Texas at Austin Center for Space Research # EIGEN: GFZ German Research Center for Geosciences (RL01-RL05) @@ -301,11 +310,11 @@ def parse_grace_file(granule): ) rx = re.compile(regex_pattern, re.VERBOSE) # extract parameters from input filename - PFX, SY, SD, EY, ED, AUX, PRC, F1, DRL, F2, SFX = rx.findall( - file_basename - ).pop() - # return the start and end date lists - return ((SY, SD), (EY, ED)) + PFX, SY, SD, EY, ED, AUX, PRC, F1, DRL, F2, SFX = rx.findall(name).pop() + # return the start and end dates + start_date = datetime(int(SY), 1, 1) + timedelta(days=int(SD) - 1) + end_date = datetime(int(EY), 1, 1) + timedelta(days=int(ED) - 1) + return (start_date, end_date) # PURPOSE: extract dates from GRAZ or Swarm files with regular expressions @@ -330,6 +339,13 @@ def parse_gfc_file(granule, PROC, DSET): - ``'GAC'``: combined non-tidal atmospheric and oceanic correction - ``'GAD'``: ocean bottom pressure product - ``'GSM'``: corrected monthly static gravity field product + + Returns + ------- + start_date: datetime + start date of gfc data file + end_date: datetime + end date of gfc data file """ # verify that filename is reduced to basename file_basename = pathlib.Path(granule).name @@ -354,10 +370,12 @@ def parse_gfc_file(granule, PROC, DSET): # extract parameters from input filename PFX, PRD, trunc, year, month, SFX = rx.findall(file_basename).pop() # number of days in each month for the calendar year - dpm = calendar_days(int(year)) + dpm = calendar_days(int(year), astype=int) # create start and end date lists - start_date = [int(year), int(month), 1, 0, 0, 0] - end_date = [int(year), int(month), dpm[int(month) - 1], 23, 59, 59] + start_date = datetime(int(year), int(month), 1, 0, 0, 0) + end_date = datetime( + int(year), int(month), dpm[int(month) - 1], 23, 59, 59 + ) elif (PROC == 'Swarm') and (DSET == 'GSM'): # regular expression operators for Swarm data regex_pattern = ( @@ -369,8 +387,8 @@ def parse_gfc_file(granule, PROC, DSET): SAT, tmp, PROD, starttime, endtime, RL, SFX = rx.findall( file_basename ).pop() - start_date, _ = parse_date_string(starttime) - end_date, _ = parse_date_string(endtime) + start_date = parse(starttime) + end_date = parse(endtime) elif (PROC == 'Swarm') and (DSET != 'GSM'): # regular expression operators for Swarm models regex_pattern = ( @@ -382,10 +400,12 @@ def parse_gfc_file(granule, PROC, DSET): # extract parameters from input filename PROD, trunc, month, year, SFX = rx.findall(file_basename).pop() # number of days in each month for the calendar year - dpm = calendar_days(int(year)) + dpm = calendar_days(int(year), astype=int) # create start and end date lists - start_date = [int(year), int(month), 1, 0, 0, 0] - end_date = [int(year), int(month), dpm[int(month) - 1], 23, 59, 59] + start_date = datetime(int(year), int(month), 1, 0, 0, 0) + end_date = datetime( + int(year), int(month), dpm[int(month) - 1], 23, 59, 59 + ) # return the start and end date lists return (start_date, end_date) @@ -589,7 +609,7 @@ def calendar_to_julian(year_decimal): # PURPOSE: gets the number of days per month for a given year -def calendar_days(year): +def calendar_days(year, astype=np.float64): """ Calculates the number of days per month for a given year @@ -597,6 +617,8 @@ def calendar_days(year): ---------- year: np.ndarray calendar year + astype: obj, default np.float64 + data type for output array Returns ------- @@ -616,9 +638,9 @@ def calendar_days(year): m4000 = year % 4000 # find indices for standard years and leap years using criteria if (m4 == 0) & (m100 != 0) | (m400 == 0) & (m4000 != 0): - return np.array(_dpm_leap, dtype=np.float64) + return np.array(_dpm_leap, dtype=astype) elif (m4 != 0) | (m100 == 0) & (m400 != 0) | (m4000 == 0): - return np.array(_dpm_stnd, dtype=np.float64) + return np.array(_dpm_stnd, dtype=astype) # PURPOSE: convert a numpy datetime array to delta times since an epoch @@ -640,7 +662,7 @@ def convert_datetime(date, epoch=_unix_epoch): """ # convert epoch to datetime variables if isinstance(epoch, (tuple, list)): - epoch = np.datetime64(datetime.datetime(*epoch)) + epoch = np.datetime64(datetime(*epoch)) elif isinstance(epoch, str): epoch = np.datetime64(parse(epoch)) # convert to delta time @@ -665,11 +687,11 @@ def convert_delta_time(delta_time, epoch1=None, epoch2=None, scale=1.0): """ # convert epochs to datetime variables if isinstance(epoch1, (tuple, list)): - epoch1 = np.datetime64(datetime.datetime(*epoch1)) + epoch1 = np.datetime64(datetime(*epoch1)) elif isinstance(epoch1, str): epoch1 = np.datetime64(parse(epoch1)) if isinstance(epoch2, (tuple, list)): - epoch2 = np.datetime64(datetime.datetime(*epoch2)) + epoch2 = np.datetime64(datetime(*epoch2)) elif isinstance(epoch2, str): epoch2 = np.datetime64(parse(epoch2)) # calculate the total difference in time in seconds @@ -734,9 +756,9 @@ def convert_calendar_dates( - 2400000.5 ) # convert epochs to datetime variables - epoch1 = np.datetime64(datetime.datetime(*_mjd_epoch)) + epoch1 = np.datetime64(datetime(*_mjd_epoch)) if isinstance(epoch, (tuple, list)): - epoch = np.datetime64(datetime.datetime(*epoch)) + epoch = np.datetime64(datetime(*epoch)) elif isinstance(epoch, str): epoch = np.datetime64(parse(epoch)) # calculate the total difference in time in days diff --git a/test/conftest.py b/test/conftest.py index be72234..ed8e1f8 100644 --- a/test/conftest.py +++ b/test/conftest.py @@ -1,15 +1,22 @@ import pytest + def pytest_addoption(parser): - parser.addoption("--username", action="store", help="NASA Earthdata username") - parser.addoption("--password", action="store", help="NASA Earthdata password") + parser.addoption( + '--username', action='store', help='NASA Earthdata username' + ) + parser.addoption( + '--password', action='store', help='NASA Earthdata password' + ) + @pytest.fixture def username(request): - """ Returns NASA Earthdata username """ - return request.config.getoption("--username") + """Returns NASA Earthdata username""" + return request.config.getoption('--username') + @pytest.fixture def password(request): - """ Returns NASA Earthdata password """ - return request.config.getoption("--password") + """Returns NASA Earthdata password""" + return request.config.getoption('--password') diff --git a/test/test_download_and_read.py b/test/test_download_and_read.py index 7352904..c905985 100644 --- a/test/test_download_and_read.py +++ b/test/test_download_and_read.py @@ -1,49 +1,69 @@ #!/usr/bin/env python -u""" +""" test_download_and_read.py (11/2021) Tests the read program to verify that coefficients are being extracted """ + import pytest import pathlib import posixpath import gravity_toolkit as gravtk + # PURPOSE: Download a GRACE file from PO.DAAC Cumulus and check that read program runs -def test_podaac_cumulus_download_and_read(username,password): +def test_podaac_cumulus_download_and_read(username, password): # find the path to the data files ids, urls, mtimes = gravtk.utilities.cmr( - mission='grace', center='CSR', release='RL06', level='L2', - product='GSM', start_date='2002-04-01', end_date='2002-04-30', - provider='POCLOUD', endpoint='data') + mission='grace', + center='CSR', + release='RL06', + level='L2', + product='GSM', + start_date='2002-04-01', + end_date='2002-04-30', + provider='POCLOUD', + endpoint='data', + ) # attempt to download the GRACE file try: # build opener for data client access URS = 'urs.earthdata.nasa.gov' - opener = gravtk.utilities.attempt_login(URS, - username=username, password=password, - authorization_header=False, verbose=True) + opener = gravtk.utilities.attempt_login( + URS, + username=username, + password=password, + authorization_header=False, + verbose=True, + ) # download and read as virtual file object FILE = gravtk.utilities.from_http(urls[0], context=None, verbose=True) except gravtk.utilities.urllib2.HTTPError as exc: pytest.xfail(exc.reason) - except EOFError as exc: - pytest.xfail("NASA Earthdata Login Error") + except (EOFError, gravtk.utilities.urllib2.URLError) as exc: + pytest.xfail('NASA Earthdata Login Error') # read as virtual file object Ylms = gravtk.read_GRACE_harmonics(FILE, 60) keys = ['time', 'start', 'end', 'clm', 'slm', 'eclm', 'eslm', 'header'] test = dict(start=2452369.5, end=2452394.5) assert all((key in Ylms.keys()) for key in keys) - assert all((Ylms[key] == val) for key,val in test.items()) - assert (Ylms['clm'][2,0] == -0.484169355584e-03) + assert all((Ylms[key] == val) for key, val in test.items()) + assert Ylms['clm'][2, 0] == -0.484169355584e-03 + # PURPOSE: Download a GRACE file from GFZ and check that read program runs def test_gfz_http_download_and_read(): - HOST=['https://isdc-data.gfz.de','grace','Level-2','CSR','RL06', - 'GSM-2_2002095-2002120_GRAC_UTCSR_BA01_0600.gz'] + HOST = [ + 'https://isdc-data.gfz.de', + 'grace', + 'Level-2', + 'CSR', + 'RL06', + 'GSM-2_2002095-2002120_GRAC_UTCSR_BA01_0600.gz', + ] # attempt to download the GRACE file try: # download and read as virtual file object - FILE = gravtk.utilities.from_http(HOST,verbose=True) + FILE = gravtk.utilities.from_http(HOST, verbose=True) except gravtk.utilities.urllib2.HTTPError as exc: pytest.xfail(exc.reason) # read as virtual file object @@ -51,14 +71,21 @@ def test_gfz_http_download_and_read(): keys = ['time', 'start', 'end', 'clm', 'slm', 'eclm', 'eslm', 'header'] test = dict(start=2452369.5, end=2452394.5) assert all((key in Ylms.keys()) for key in keys) - assert all((Ylms[key] == val) for key,val in test.items()) - assert (Ylms['clm'][2,0] == -0.484169355584e-03) + assert all((Ylms[key] == val) for key, val in test.items()) + assert Ylms['clm'][2, 0] == -0.484169355584e-03 + # PURPOSE: Download a GRACE file from GFZ and check that read program runs -@pytest.mark.skip(reason="Deprecated GFZ FTP server") +@pytest.mark.skip(reason='Deprecated GFZ FTP server') def test_gfz_ftp_download_and_read(): - HOST=['isdcftp.gfz-potsdam.de','grace','Level-2','CSR','RL06', - 'GSM-2_2002095-2002120_GRAC_UTCSR_BA01_0600.gz'] + HOST = [ + 'isdcftp.gfz-potsdam.de', + 'grace', + 'Level-2', + 'CSR', + 'RL06', + 'GSM-2_2002095-2002120_GRAC_UTCSR_BA01_0600.gz', + ] # attempt to download the GRACE file try: # download and read as virtual file object @@ -70,17 +97,23 @@ def test_gfz_ftp_download_and_read(): keys = ['time', 'start', 'end', 'clm', 'slm', 'eclm', 'eslm', 'header'] test = dict(start=2452369.5, end=2452394.5) assert all((key in Ylms.keys()) for key in keys) - assert all((Ylms[key] == val) for key,val in test.items()) - assert (Ylms['clm'][2,0] == -0.484169355584e-03) + assert all((Ylms[key] == val) for key, val in test.items()) + assert Ylms['clm'][2, 0] == -0.484169355584e-03 + # PURPOSE: Download a GRACE-FO COST-G file from the GFZ ICGEM def test_gfz_icgem_costg_download_and_read(): - HOST=['https://icgem.gfz.de','getseries','02_COST-G_', - 'Grace-FO_RL02','GSM-2_2018152-2018181_GRFO_COSTG_BF01_0200.gfc'] + HOST = [ + 'https://icgem.gfz.de', + 'getseries', + '02_COST-G_', + 'Grace-FO_RL02', + 'GSM-2_2018152-2018181_GRFO_COSTG_BF01_0200.gfc', + ] # attempt to download the GRACE file try: # download and read as virtual file object - FILE = gravtk.utilities.from_http(HOST,verbose=True) + FILE = gravtk.utilities.from_http(HOST, verbose=True) except gravtk.utilities.urllib2.HTTPError as exc: pytest.xfail(exc.reason) # read as virtual file object @@ -88,22 +121,23 @@ def test_gfz_icgem_costg_download_and_read(): keys = ['time', 'start', 'end', 'clm', 'slm', 'eclm', 'eslm', 'header'] test = dict(start=2458270.5, end=2458299.5) assert all((key in Ylms.keys()) for key in keys) - assert all((Ylms[key] == val) for key,val in test.items()) - assert (Ylms['clm'][2,0] == -0.484165368910e-03) + assert all((Ylms[key] == val) for key, val in test.items()) + assert Ylms['clm'][2, 0] == -0.484165368910e-03 + # PURPOSE: Download a Swarm file from ESA and check that read program runs def test_esa_swarm_download_and_read(): # build url for Swarm file - HOST='https://swarm-diss.eo.esa.int' - swarm_file='SW_OPER_EGF_SHA_2__20131201T000000_20131231T235959_0101.ZIP' - parameters = gravtk.utilities.urlencode({'file': - posixpath.join('swarm','Level2longterm','EGF',swarm_file)}) - remote_file = [HOST,'?do=download&{0}'.format(parameters)] + HOST = 'https://swarm-diss.eo.esa.int' + swarm_file = 'SW_OPER_EGF_SHA_2__20131201T000000_20131231T235959_0101.ZIP' + parameters = gravtk.utilities.urlencode( + {'file': posixpath.join('swarm', 'Level2longterm', 'EGF', swarm_file)} + ) + remote_file = [HOST, '?do=download&{0}'.format(parameters)] # attempt to download the Swarm file try: # download as local file object - gravtk.utilities.from_http(remote_file, - local=swarm_file, verbose=True) + gravtk.utilities.from_http(remote_file, local=swarm_file, verbose=True) except gravtk.utilities.urllib2.HTTPError as exc: pytest.xfail(exc.reason) # read the local file @@ -112,16 +146,24 @@ def test_esa_swarm_download_and_read(): keys = ['time', 'start', 'end', 'clm', 'slm', 'eclm', 'eslm'] test = dict(start=2456627.5, end=2456658.499988426) assert all((key in Ylms.keys()) for key in keys) - assert all((Ylms[key] == val) for key,val in test.items()) - assert (Ylms['clm'][2,0] == -0.48416530506600003e-03) + assert all((Ylms[key] == val) for key, val in test.items()) + assert Ylms['clm'][2, 0] == -0.48416530506600003e-03 # clean up swarm_file.unlink() + # PURPOSE: Download a GRACE ITSG GRAZ file and check that read program runs def test_itsg_graz_download_and_read(): - HOST=['http://ftp.tugraz.at','pub','ITSG','GRACE', - 'ITSG-Grace_operational','monthly','monthly_n60', - 'ITSG-Grace_operational_n60_2018-06.gfc'] + HOST = [ + 'http://ftp.tugraz.at', + 'pub', + 'ITSG', + 'GRACE', + 'ITSG-Grace_operational', + 'monthly', + 'monthly_n60', + 'ITSG-Grace_operational_n60_2018-06.gfc', + ] # attempt to download the GRAZ file try: # download as local file object @@ -134,7 +176,7 @@ def test_itsg_graz_download_and_read(): keys = ['time', 'start', 'end', 'clm', 'slm', 'eclm', 'eslm'] test = dict(start=2458270.5, end=2458300.499988426) assert all((key in Ylms.keys()) for key in keys) - assert all((Ylms[key] == val) for key,val in test.items()) - assert (Ylms['clm'][2,0] == -0.4841694727612e-03) + assert all((Ylms[key] == val) for key, val in test.items()) + assert Ylms['clm'][2, 0] == -0.4841694727612e-03 # clean up itsg_file.unlink() diff --git a/test/test_gia.py b/test/test_gia.py index 9abe91a..06e2f0b 100644 --- a/test/test_gia.py +++ b/test/test_gia.py @@ -1,8 +1,9 @@ #!/usr/bin/env python -u""" +""" test_gia.py (12/2022) Tests that the GIA model readers are equivalent """ + import gzip import time import pytest @@ -11,15 +12,27 @@ import numpy as np import gravity_toolkit as gravtk + # PURPOSE: Download ICE-6G GIA model -@pytest.fixture(scope="module", autouse=True) +@pytest.fixture(scope='module', autouse=True) def download_GIA_model(): # output GIA file GIA_FILE = pathlib.Path('Stokes_trend_High_Res.txt') # download GIA model - HOST = ['https://www.atmosp.physics.utoronto.ca','~peltier','datasets', - 'Ice6G_C_VM5a','ICE-6G_High_Res_Stokes_trend.txt.gz'] - fid = gravtk.utilities.from_http(HOST, verbose=True) + HOST = [ + 'https://www.atmosp.physics.utoronto.ca', + '~peltier', + 'datasets', + 'Ice6G_C_VM5a', + 'ICE-6G_High_Res_Stokes_trend.txt.gz', + ] + # attempt to download the GIA file + try: + fid = gravtk.utilities.from_http(HOST, verbose=True) + except gravtk.utilities.urllib2.HTTPError as exc: + pytest.xfail(exc.reason) + except gravtk.utilities.urllib2.URLError as exc: + pytest.xfail('Connection Refused') # decompress GIA model from virtual BytesIO object with gzip.open(fid, 'rb') as f_in, open(GIA_FILE, 'wb') as f_out: shutil.copyfileobj(f_in, f_out) @@ -28,6 +41,7 @@ def download_GIA_model(): # clean up GIA_FILE.unlink() + # PURPOSE: read ICE-6G GIA test outputs def test_GIA_model_read(): # output GIA file and type @@ -36,25 +50,26 @@ def test_GIA_model_read(): # read GIA model Ylms = gravtk.read_GIA_model(GIA_FILE, GIA=GIA) # assert input GIA values - assert Ylms['clm'][2,0] == 1.43961238E-11 - assert Ylms['clm'][3,0] == 1.52009079E-12 - assert Ylms['slm'][3,1] == -8.05198489E-12 + assert Ylms['clm'][2, 0] == 1.43961238e-11 + assert Ylms['clm'][3, 0] == 1.52009079e-12 + assert Ylms['slm'][3, 1] == -8.05198489e-12 # assert parameters assert Ylms['title'] == 'ICE6G-D_High_Res' + # PURPOSE: read ICE-6G GIA model and test harmonic outputs def test_GIA_model_harmonics(): # output GIA file and type GIA_FILE = 'Stokes_trend_High_Res.txt' GIA = 'ICE6G-D' # degree of truncation - LMAX,MMAX = (60, 30) + LMAX, MMAX = (60, 30) # read GIA model Ylms = gravtk.gia(lmax=LMAX).from_GIA(GIA_FILE, GIA=GIA, mmax=MMAX) # assert input GIA values - assert Ylms.clm[2,0] == 1.43961238E-11 - assert Ylms.clm[3,0] == 1.52009079E-12 - assert Ylms.slm[3,1] == -8.05198489E-12 + assert Ylms.clm[2, 0] == 1.43961238e-11 + assert Ylms.clm[3, 0] == 1.52009079e-12 + assert Ylms.slm[3, 1] == -8.05198489e-12 # assert parameters assert Ylms.title == 'ICE6G-D_High_Res' # assert truncation @@ -63,36 +78,42 @@ def test_GIA_model_harmonics(): assert Ylms.mmax == MMAX assert Ylms.m[-1] == MMAX + # PURPOSE: read ICE-6G GIA model and compare drift estimates def test_GIA_model_drift_estimate(): # output GIA file and type GIA_FILE = 'Stokes_trend_High_Res.txt' GIA = 'ICE6G-D' # degree and order of truncation - LMAX,MMAX = (60, 30) + LMAX, MMAX = (60, 30) # synthetic time estimate now = time.gmtime() - tdec = np.arange(2002, now.tm_year+1, 1.0/12.0) + tdec = np.arange(2002, now.tm_year + 1, 1.0 / 12.0) epoch = 2003.3 # read GIA model - GIA_Ylms_rate = gravtk.read_GIA_model(GIA_FILE, GIA=GIA, LMAX=LMAX, MMAX=MMAX) + GIA_Ylms_rate = gravtk.read_GIA_model( + GIA_FILE, GIA=GIA, LMAX=LMAX, MMAX=MMAX + ) # calculate the monthly mass change from GIA GIA_Ylms = gravtk.harmonics(lmax=LMAX, mmax=MMAX) GIA_Ylms.time = np.copy(tdec) GIA_Ylms.month = gravtk.time.calendar_to_grace(tdec) GIA_Ylms.month = gravtk.time.adjust_months(GIA_Ylms.month) # allocate for output harmonics - GIA_Ylms.clm = np.zeros((GIA_Ylms.lmax+1, GIA_Ylms.mmax+1, len(tdec))) - GIA_Ylms.slm = np.zeros((GIA_Ylms.lmax+1, GIA_Ylms.mmax+1, len(tdec))) + GIA_Ylms.clm = np.zeros((GIA_Ylms.lmax + 1, GIA_Ylms.mmax + 1, len(tdec))) + GIA_Ylms.slm = np.zeros((GIA_Ylms.lmax + 1, GIA_Ylms.mmax + 1, len(tdec))) # assert input GIA values # monthly GIA calculated by gia_rate*time elapsed # finding change in GIA each month - for i,t in enumerate(tdec): - GIA_Ylms.clm[:,:,i] = GIA_Ylms_rate['clm']*(t - epoch) - GIA_Ylms.slm[:,:,i] = GIA_Ylms_rate['slm']*(t - epoch) + for i, t in enumerate(tdec): + GIA_Ylms.clm[:, :, i] = GIA_Ylms_rate['clm'] * (t - epoch) + GIA_Ylms.slm[:, :, i] = GIA_Ylms_rate['slm'] * (t - epoch) # read GIA model and calculate drift from harmonics class - Ylms = gravtk.gia(lmax=LMAX).from_GIA( - GIA_FILE, GIA=GIA, mmax=MMAX).drift(tdec, epoch=epoch) + Ylms = ( + gravtk.gia(lmax=LMAX) + .from_GIA(GIA_FILE, GIA=GIA, mmax=MMAX) + .drift(tdec, epoch=epoch) + ) # assert that spherical harmonics are equal assert np.all(GIA_Ylms.clm == Ylms.clm) assert np.all(GIA_Ylms.slm == Ylms.slm) diff --git a/test/test_harmonics.py b/test/test_harmonics.py index 34de7e0..dcb8def 100755 --- a/test/test_harmonics.py +++ b/test/test_harmonics.py @@ -1,5 +1,5 @@ #!/usr/bin/env python -u""" +""" test_harmonics.py (08/2020) Tests harmonic programs using the Velicogna and Wahr (2013) Greenland synthetic 1. Converts synthetic spatial distribution to spherical harmonics @@ -8,6 +8,7 @@ 4. Compares output smoothed spatial distribution with validation dataset Tests harmonic objects flatten, expansion and iteration routines """ + import pytest import inspect import pathlib @@ -18,19 +19,27 @@ filename = inspect.getframeinfo(inspect.currentframe()).filename filepath = pathlib.Path(filename).absolute().parent + # PURPOSE: test harmonic conversion programs def test_harmonics(): # path to load Love numbers file - love_numbers_file = gravtk.utilities.get_data_path( - ['data','love_numbers']) + love_numbers_file = gravtk.utilities.get_data_path(['data', 'love_numbers']) # read load Love numbers LOVE = gravtk.read_love_numbers(love_numbers_file, FORMAT='class') # read input spatial distribution file - distribution_file = filepath.joinpath('out.green_ice.grid.0.5.2008.cmh20.gz') - input_distribution = gravtk.spatial().from_ascii(distribution_file, - date=False, spacing=[0.5,0.5], nlat=361, nlon=721, - extent=[0,360.0,-90,90], compression='gzip') + distribution_file = filepath.joinpath( + 'out.green_ice.grid.0.5.2008.cmh20.gz' + ) + input_distribution = gravtk.spatial().from_ascii( + distribution_file, + date=False, + spacing=[0.5, 0.5], + nlat=361, + nlon=721, + extent=[0, 360.0, -90, 90], + compression='gzip', + ) # spherical harmonic parameters # maximum spherical harmonic degree @@ -44,17 +53,25 @@ def test_harmonics(): theta[theta > np.arccos(-0.9999999)] = np.arccos(-0.9999999) theta[theta < np.arccos(0.9999999)] = np.arccos(0.9999999) # calculate Legendre polynomials with Martin Mohlenkamp's relation - PLM,dPLM = gravtk.associated_legendre(LMAX, np.cos(theta), - method='mohlenkamp') + PLM, dPLM = gravtk.associated_legendre( + LMAX, np.cos(theta), method='mohlenkamp' + ) # convert to spherical harmonics - test_Ylms = gravtk.gen_stokes(input_distribution.data, - input_distribution.lon, input_distribution.lat, UNITS=1, LMAX=LMAX, - PLM=PLM, LOVE=LOVE) + test_Ylms = gravtk.gen_stokes( + input_distribution.data, + input_distribution.lon, + input_distribution.lat, + UNITS=1, + LMAX=LMAX, + PLM=PLM, + LOVE=LOVE, + ) # read harmonics from file harmonics_file = filepath.joinpath('out.geoid.green_ice.0.5.2008.60.gz') valid_Ylms = gravtk.harmonics(lmax=LMAX, mmax=LMAX).from_ascii( - harmonics_file, date=False, compression='gzip') + harmonics_file, date=False, compression='gzip' + ) # check that harmonic data is equal to machine precision difference_Ylms = test_Ylms.copy() @@ -66,27 +83,51 @@ def test_harmonics(): # cmwe, centimeters water equivalent smooth_Ylms = test_Ylms.copy() dfactor = gravtk.units(lmax=LMAX).harmonic(*LOVE) - wt = 2.0*np.pi*gravtk.gauss_weights(RAD,LMAX) - smooth_Ylms.convolve(dfactor.cmwe*wt) + wt = 2.0 * np.pi * gravtk.gauss_weights(RAD, LMAX) + smooth_Ylms.convolve(dfactor.cmwe * wt) transform_Ylms = smooth_Ylms.copy() # convert harmonics back to spatial domain at same grid spacing - test_distribution = gravtk.harmonic_summation(smooth_Ylms.clm, - smooth_Ylms.slm, input_distribution.lon, input_distribution.lat, - LMAX=LMAX, PLM=PLM).T + test_distribution = gravtk.harmonic_summation( + smooth_Ylms.clm, + smooth_Ylms.slm, + input_distribution.lon, + input_distribution.lat, + LMAX=LMAX, + PLM=PLM, + ).T # convert harmonics using fast-fourier transform method - test_transform = gravtk.harmonic_transform(transform_Ylms.clm, - transform_Ylms.slm, input_distribution.lon, input_distribution.lat, - LMAX=LMAX, PLM=PLM).T + test_transform = gravtk.harmonic_transform( + transform_Ylms.clm, + transform_Ylms.slm, + input_distribution.lon, + input_distribution.lat, + LMAX=LMAX, + PLM=PLM, + ).T # convert harmonics back to spatial domain using wrapper function - test_combine = gravtk.stokes_summation(test_Ylms.clm, - test_Ylms.slm, input_distribution.lon, input_distribution.lat, - LMAX=LMAX, PLM=PLM, UNITS=1, RAD=RAD, LOVE=LOVE).T + test_combine = gravtk.stokes_summation( + test_Ylms.clm, + test_Ylms.slm, + input_distribution.lon, + input_distribution.lat, + LMAX=LMAX, + PLM=PLM, + UNITS=1, + RAD=RAD, + LOVE=LOVE, + ).T # read input and output spatial distribution files spatial_file = filepath.joinpath('out.combine.green_ice.0.5.2008.60.gz') - output_distribution = gravtk.spatial().from_ascii(spatial_file, - date=False, spacing=[0.5,0.5], nlat=361, nlon=721, - extent=[0,360.0,-90,90], compression='gzip') + output_distribution = gravtk.spatial().from_ascii( + spatial_file, + date=False, + spacing=[0.5, 0.5], + nlat=361, + nlon=721, + extent=[0, 360.0, -90, 90], + compression='gzip', + ) # check that data is equal to machine precision difference_distribution = test_distribution - output_distribution.data @@ -96,33 +137,34 @@ def test_harmonics(): assert np.all(np.abs(difference_transform) < distribution_eps) assert np.all(test_distribution == test_combine) + # PURPOSE: test harmonic objects def test_iterate(): # maximum spherical harmonic degree and order LMAX = 60 MMAX = 30 # number of harmonics - n_harm = (LMAX**2 + 3*LMAX - (LMAX-MMAX)**2 - (LMAX-MMAX))//2 + 1 + n_harm = (LMAX**2 + 3 * LMAX - (LMAX - MMAX) ** 2 - (LMAX - MMAX)) // 2 + 1 # number of time points for test nt = 12 # create flattened test harmonics flat_Ylms = gravtk.harmonics(lmax=LMAX, mmax=MMAX) - ll,mm = np.meshgrid(np.arange(LMAX+1), np.arange(MMAX+1)) - ii,jj = np.triu_indices(MMAX+1, k=0, m=LMAX+1) - flat_Ylms.l = ll[ii,jj].astype(np.int64) - flat_Ylms.m = mm[ii,jj].astype(np.int64) + ll, mm = np.meshgrid(np.arange(LMAX + 1), np.arange(MMAX + 1)) + ii, jj = np.triu_indices(MMAX + 1, k=0, m=LMAX + 1) + flat_Ylms.l = ll[ii, jj].astype(np.int64) + flat_Ylms.m = mm[ii, jj].astype(np.int64) # add date variables flat_Ylms.month = 1 + np.arange(nt) - flat_Ylms.time = 2002.0 + (flat_Ylms.month - 0.5)/12.0 + flat_Ylms.time = 2002.0 + (flat_Ylms.month - 0.5) / 12.0 # create random harmonics - flat_Ylms.clm = np.random.rand(n_harm,nt) - flat_Ylms.slm = np.random.rand(n_harm,nt) + flat_Ylms.clm = np.random.rand(n_harm, nt) + flat_Ylms.slm = np.random.rand(n_harm, nt) # reshape harmonics object valid_Ylms = flat_Ylms.expand(date=True) # iterate over harmonics - for i,test in enumerate(valid_Ylms): + for i, test in enumerate(valid_Ylms): valid = valid_Ylms.index(i) assert np.isclose(valid.l, test.l).all() assert np.isclose(valid.m, test.m).all() diff --git a/test/test_legendre.py b/test/test_legendre.py index b5fb5e9..1606dad 100644 --- a/test/test_legendre.py +++ b/test/test_legendre.py @@ -1,51 +1,62 @@ #!/usr/bin/env python -u""" +""" test_legendre.py (11/2021) """ + import numpy as np import gravity_toolkit as gravtk + # PURPOSE: test unnormalized Legendre polynomials def test_unnormalized(l=3, x=[-1.0, -0.9, -0.8]): obs = gravtk.legendre(l, x) - expected = np.array([ - [-1.00000, -0.47250, -0.08000], - [0.00000, -1.99420, -1.98000], - [0.00000, -2.56500, -4.32000], - [0.00000, -1.24229, -3.24000] - ]) + expected = np.array( + [ + [-1.00000, -0.47250, -0.08000], + [0.00000, -1.99420, -1.98000], + [0.00000, -2.56500, -4.32000], + [0.00000, -1.24229, -3.24000], + ] + ) assert np.isclose(obs, expected, atol=1e-05).all() + # PURPOSE: test fully-normalized Legendre polynomials def test_normalized(l=3, x=[-1.0, -0.9, -0.8]): obs = gravtk.legendre(l, x, NORMALIZE=True) - expected = np.array([ - [-2.64575, -1.25012, -0.21166], - [-0.00000, 2.15398, 2.13864], - [0.00000, -0.87611, -1.47556], - [-0.00000, 0.17323, 0.45180] - ]) + expected = np.array( + [ + [-2.64575, -1.25012, -0.21166], + [-0.00000, 2.15398, 2.13864], + [0.00000, -0.87611, -1.47556], + [-0.00000, 0.17323, 0.45180], + ] + ) assert np.isclose(obs, expected, atol=1e-05).all() + # PURPOSE: test fully-normalized zonal Legendre polynomials def test_zonal(l=3, x=[-1.0, -0.9, -0.8]): - obs,_ = gravtk.legendre_polynomials(l, x) - expected = np.array([ - [1.00000, 1.00000, 1.00000], - [-1.73205, -1.55885, -1.38564], - [2.23607, 1.59879, 1.02859], - [-2.64575, -1.25012, -0.21166], - ]) + obs, _ = gravtk.legendre_polynomials(l, x) + expected = np.array( + [ + [1.00000, 1.00000, 1.00000], + [-1.73205, -1.55885, -1.38564], + [2.23607, 1.59879, 1.02859], + [-2.64575, -1.25012, -0.21166], + ] + ) assert np.isclose(obs, expected, atol=1e-05).all() + # PURPOSE: compare fully-normalized Legendre polynomials def test_plms(l=240, x=0.1): obs = gravtk.legendre(l, x, NORMALIZE=True) # calculate associated Legendre polynomials - holmes,_ = gravtk.plm_holmes(l, x) - colombo,_ = gravtk.plm_colombo(l, x) - mohlenkamp,_ = gravtk.plm_mohlenkamp(l, x) + holmes, _ = gravtk.plm_holmes(l, x) + colombo, _ = gravtk.plm_colombo(l, x) + mohlenkamp, _ = gravtk.plm_mohlenkamp(l, x) # compare Legendre polynomials - assert np.isclose(obs, holmes[l,:]).all() + assert np.isclose(obs, holmes[l, :]).all() assert np.isclose(holmes, colombo).all() assert np.isclose(holmes, mohlenkamp).all() diff --git a/test/test_love_numbers.py b/test/test_love_numbers.py index 0d9c558..60a107e 100644 --- a/test/test_love_numbers.py +++ b/test/test_love_numbers.py @@ -1,5 +1,5 @@ #!/usr/bin/env python -u""" +""" test_love_numbers.py (01/2023) UPDATE HISTORY: Updated 01/2023: single implicit import of gravity toolkit @@ -7,10 +7,12 @@ add tests for Gegout and Wang load Love number sets Written 08/2020 """ + import os import pytest import gravity_toolkit as gravtk + # PURPOSE: Define Load Love Numbers for lower degree harmonics def get_love_numbers(): """ @@ -21,38 +23,56 @@ def get_love_numbers(): ll = [0.0, 0.13026466961444, 0.023882296795977, 0.069842389427609] return dict(hl=hl, kl=kl, ll=ll) + # PURPOSE: Check that Load Love Numbers match expected for reference frame def test_love_numbers(): # valid low degree Love numbers for reference frame CF VALID = get_love_numbers() # path to load Love numbers file - love_numbers_file = gravtk.utilities.get_data_path( - ['data','love_numbers']) + love_numbers_file = gravtk.utilities.get_data_path(['data', 'love_numbers']) # read load Love numbers and convert to reference frame CF - TEST = gravtk.read_love_numbers(love_numbers_file, - LMAX=1000, HEADER=2, FORMAT='dict', REFERENCE='CF') - assert all((v==t).all() for key in ['hl','kl','ll'] - for v,t in zip(VALID[key],TEST[key])) - assert (TEST['l'].max() == 1000) + TEST = gravtk.read_love_numbers( + love_numbers_file, LMAX=1000, HEADER=2, FORMAT='dict', REFERENCE='CF' + ) + assert all( + (v == t).all() + for key in ['hl', 'kl', 'll'] + for v, t in zip(VALID[key], TEST[key]) + ) + assert TEST['l'].max() == 1000 + # PURPOSE: Check that Gegout (2005) Load Love Numbers can be read def test_Gegout_love_numbers(): # path to load Love numbers file love_numbers_file = gravtk.utilities.get_data_path( - ['data','Load_Love2_CE.dat']) - COLUMNS = ['l','hl','ll','kl'] + ['data', 'Load_Love2_CE.dat'] + ) + COLUMNS = ['l', 'hl', 'll', 'kl'] # read load Love numbers and convert to reference frame CM - TEST = gravtk.read_love_numbers(love_numbers_file, - HEADER=3, COLUMNS=COLUMNS, FORMAT='class', REFERENCE='CM') - assert (TEST.lmax == 1024) + TEST = gravtk.read_love_numbers( + love_numbers_file, + HEADER=3, + COLUMNS=COLUMNS, + FORMAT='class', + REFERENCE='CM', + ) + assert TEST.lmax == 1024 + # PURPOSE: Check that Wang et al. (2012) Load Love Numbers can be read def test_Wang_love_numbers(): # path to load Love numbers file (truncated from degree 46341) love_numbers_file = gravtk.utilities.get_data_path( - ['data','PREM-LLNs-truncated.dat']) - COLUMNS = ['l','hl','ll','kl','nl','nk'] + ['data', 'PREM-LLNs-truncated.dat'] + ) + COLUMNS = ['l', 'hl', 'll', 'kl', 'nl', 'nk'] # read load Love numbers and convert to reference frame CE - TEST = gravtk.read_love_numbers(love_numbers_file, - HEADER=1, COLUMNS=COLUMNS, FORMAT='class', REFERENCE='CE') - assert (TEST.lmax == 5000) + TEST = gravtk.read_love_numbers( + love_numbers_file, + HEADER=1, + COLUMNS=COLUMNS, + FORMAT='class', + REFERENCE='CE', + ) + assert TEST.lmax == 5000 diff --git a/test/test_masks.py b/test/test_masks.py index d04cb87..5b6d56c 100644 --- a/test/test_masks.py +++ b/test/test_masks.py @@ -1,25 +1,30 @@ #!/usr/bin/env python -u""" +""" test_masks.py (06/2025) Tests that the stored masks can be read """ + import pytest import numpy as np import gravity_toolkit as gravtk # parameterize land-sea masks -_land_sea_masks = ["landsea_1d.nc", "landsea_hd.nc", "landsea_qd.nc"] -@pytest.mark.parametrize("LANDMASK", _land_sea_masks) +_land_sea_masks = ['landsea_1d.nc', 'landsea_hd.nc', 'landsea_qd.nc'] + + +@pytest.mark.parametrize('LANDMASK', _land_sea_masks) def test_lsmask(LANDMASK): - u""" + """ Test that the land mask can be read and that the area of the ocean is close to the expected value """ # Land-Sea Mask with Antarctica from Rignot (2017) and Greenland from GEUS # 0=Ocean, 1=Land, 2=Lake, 3=Small Island, 4=Ice Shelf # Open the land-sea NetCDF file for reading - lsmask = gravtk.utilities.get_data_path(['data',LANDMASK]) - landsea = gravtk.spatial().from_netCDF4(lsmask, date=False, varname='LSMASK') + lsmask = gravtk.utilities.get_data_path(['data', LANDMASK]) + landsea = gravtk.spatial().from_netCDF4( + lsmask, date=False, varname='LSMASK' + ) # degree spacing and grid dimensions dlon, dlat = landsea.spacing nlat, nlon = landsea.shape @@ -32,25 +37,27 @@ def test_lsmask(LANDMASK): # create land function land_function = np.zeros((nlat, nlon), dtype=np.float64) # combine land and island levels for land function - indx,indy = np.nonzero((landsea.data >= 1) & (landsea.data <= 3)) - land_function[indx,indy] = 1.0 + indx, indy = np.nonzero((landsea.data >= 1) & (landsea.data <= 3)) + land_function[indx, indy] = 1.0 # calculate the ocean function ocean_function = 1.0 - land_function # average Radius of the Earth [km] - rad_e = gravtk.units().rad_e/1e5 + rad_e = gravtk.units().rad_e / 1e5 # total area of ocean calculated by integrating the ocean function - area = np.sum(ocean_function*np.sin(th)*dphi*dth*rad_e**2) + area = np.sum(ocean_function * np.sin(th) * dphi * dth * rad_e**2) # assert that the area is close to the expected value assert np.isclose(area, 3.62e8, rtol=0.01) + def test_smoothed(): """ Test that the smoothed land mask can be read """ # Read Smoothed Ocean and Land Functions # will mask out land regions in the final current maps - LANDMASK = gravtk.utilities.get_data_path(['data','land_fcn_300km.nc']) - landsea = gravtk.spatial().from_netCDF4(LANDMASK, - date=False, varname='LSMASK') + LANDMASK = gravtk.utilities.get_data_path(['data', 'land_fcn_300km.nc']) + landsea = gravtk.spatial().from_netCDF4( + LANDMASK, date=False, varname='LSMASK' + ) assert landsea.shape == (180, 360) assert landsea.spacing == (1.0, 1.0) diff --git a/test/test_point_masses.py b/test/test_point_masses.py index 10b1b17..79b0da3 100644 --- a/test/test_point_masses.py +++ b/test/test_point_masses.py @@ -1,45 +1,46 @@ #!/usr/bin/env python -u""" +""" test_point_masses.py (01/2023) UPDATE HISTORY: Updated 01/2023: single implicit import of gravity toolkit Written 02/2021 """ + import pytest import numpy as np import gravity_toolkit as gravtk + # parameterize the number of point masses -@pytest.mark.parametrize("NPTS", np.random.randint(2,2000,size=1)) +@pytest.mark.parametrize('NPTS', np.random.randint(2, 2000, size=1)) def test_point_masses(NPTS): # create spatial grid - dlon,dlat = (1.0,1.0) - lat = np.arange(90.0 - dlat/2.0, -90.0 - dlat/2.0, -dlat) - lon = np.arange(-180.0 + dlon/2.0, 180.0 + dlon/2.0, dlon) - gridlon,gridlat = np.meshgrid(lon,lat) + dlon, dlat = (1.0, 1.0) + lat = np.arange(90.0 - dlat / 2.0, -90.0 - dlat / 2.0, -dlat) + lon = np.arange(-180.0 + dlon / 2.0, 180.0 + dlon / 2.0, dlon) + gridlon, gridlat = np.meshgrid(lon, lat) nlat, nlon = np.shape(gridlon) # parameterize point masses - LAT = lat[0]-dlat*np.random.randint(0,nlat,size=NPTS) - LON = lon[0]+dlon*np.random.randint(0,nlon,size=NPTS) - MASS = 100.0 - 200.0*np.random.randn(NPTS) + LAT = lat[0] - dlat * np.random.randint(0, nlat, size=NPTS) + LON = lon[0] + dlon * np.random.randint(0, nlon, size=NPTS) + MASS = 100.0 - 200.0 * np.random.randn(NPTS) # create test gridded field data = np.zeros((nlat, nlon)) for i in range(NPTS): - indy,indx = np.nonzero((gridlat == LAT[i]) & (gridlon == LON[i])) - data[indy,indx] += MASS[i] + indy, indx = np.nonzero((gridlat == LAT[i]) & (gridlon == LON[i])) + data[indy, indx] += MASS[i] # path to load Love numbers file - love_numbers_file = gravtk.utilities.get_data_path( - ['data','love_numbers']) + love_numbers_file = gravtk.utilities.get_data_path(['data', 'love_numbers']) # read load Love numbers LOVE = gravtk.read_love_numbers(love_numbers_file) # calculate harmonics and degree amplitudes for each case - grid_Ylms = gravtk.gen_stokes(data, lon, lat, - LMAX=60, UNITS=2, LOVE=LOVE) - point_Ylms = gravtk.gen_point_load(MASS, LON, LAT, - LMAX=60, UNITS=2, LOVE=LOVE) + grid_Ylms = gravtk.gen_stokes(data, lon, lat, LMAX=60, UNITS=2, LOVE=LOVE) + point_Ylms = gravtk.gen_point_load( + MASS, LON, LAT, LMAX=60, UNITS=2, LOVE=LOVE + ) # check that harmonic data is equal to machine precision difference_Ylms = grid_Ylms.copy() diff --git a/test/test_sea_level.py b/test/test_sea_level.py index ebe093a..a06c68f 100644 --- a/test/test_sea_level.py +++ b/test/test_sea_level.py @@ -1,7 +1,8 @@ #!/usr/bin/env python -u""" +""" test_sea_level.py (07/2026) """ + import inspect import pathlib import numpy as np @@ -11,19 +12,25 @@ filename = inspect.getframeinfo(inspect.currentframe()).filename filepath = pathlib.Path(filename).absolute().parent + # PURPOSE: test sea level equation programs def test_sea_level(): # path to load Love numbers file - love_numbers_file = gravtk.utilities.get_data_path( - ['data','love_numbers']) + love_numbers_file = gravtk.utilities.get_data_path(['data', 'love_numbers']) # read load Love numbers LOVE = gravtk.read_love_numbers(love_numbers_file, FORMAT='class') # read land function file LANDMASK = filepath.joinpath('land.fcn.1_deg.gz') - landsea = gravtk.spatial().from_ascii(LANDMASK, - date=False, spacing=[1.0, 1.0], nlat=180, nlon=360, - extent=[0.5,359.5,-89.5,89.5], compression='gzip') + landsea = gravtk.spatial().from_ascii( + LANDMASK, + date=False, + spacing=[1.0, 1.0], + nlat=180, + nlon=360, + extent=[0.5, 359.5, -89.5, 89.5], + compression='gzip', + ) # spherical harmonic parameters # maximum spherical harmonic degree @@ -31,44 +38,75 @@ def test_sea_level(): # read harmonics from file harmonics_file = filepath.joinpath('out.geoid.green_ice.0.5.2008.60.gz') Ylms = gravtk.harmonics(lmax=LMAX, mmax=LMAX).from_ascii( - harmonics_file, date=False, compression='gzip') + harmonics_file, date=False, compression='gzip' + ) # calculate the legendre functions using Martin Mohlenkamp's relation th = np.radians(90.0 - landsea.lat) PLM, dPLM = gravtk.plm_mohlenkamp(LMAX, np.cos(th)) # run pseudo-spectral sea level equation solver - sea_level = gravtk.sea_level_equation(Ylms.clm, Ylms.slm, - landsea.lon, landsea.lat, landsea.data.T, LMAX=LMAX, - LOVE=LOVE, BODY_TIDE_LOVE=0, - FLUID_LOVE=0, DENSITY=1.0, POLAR=True, - PLM=PLM, ITERATIONS=2, FILL_VALUE=np.nan).T - + sea_level = gravtk.sea_level_equation( + Ylms.clm, + Ylms.slm, + landsea.lon, + landsea.lat, + landsea.data.T, + LMAX=LMAX, + LOVE=LOVE, + BODY_TIDE_LOVE=0, + FLUID_LOVE=0, + DENSITY=1.0, + POLAR=True, + PLM=PLM, + ITERATIONS=2, + FILL_VALUE=np.nan, + ).T + # check that sea level data is equal to file precision valid_file = filepath.joinpath('out.slf.green_ice.1_deg.2008.60.gz') - validation = gravtk.spatial().from_ascii(valid_file, - date=False, spacing=[1.0, 1.0], nlat=180, nlon=360, - extent=[0.5,359.5,-89.5,89.5], compression='gzip') + validation = gravtk.spatial().from_ascii( + valid_file, + date=False, + spacing=[1.0, 1.0], + nlat=180, + nlon=360, + extent=[0.5, 359.5, -89.5, 89.5], + compression='gzip', + ) # check differences difference = validation.data - sea_level valid_difference = difference[np.isfinite(difference)] assert np.all(np.abs(valid_difference) < 1e-8) + def test_harmonics(): # read land function file LANDMASK = filepath.joinpath('land.fcn.1_deg.gz') - landsea = gravtk.spatial().from_ascii(LANDMASK, - date=False, spacing=[1.0, 1.0], nlat=180, nlon=360, - extent=[0.5,359.5,-89.5,89.5], compression='gzip') + landsea = gravtk.spatial().from_ascii( + LANDMASK, + date=False, + spacing=[1.0, 1.0], + nlat=180, + nlon=360, + extent=[0.5, 359.5, -89.5, 89.5], + compression='gzip', + ) # calculate ocean function from land function land_function = landsea.data.T ocean_function = 1.0 - land_function # maximum spherical harmonic degree LMAX = 60 # calculate spherical harmonics using integration and fourier methods - YlmI = gravtk.gen_harmonics(ocean_function, landsea.lon, landsea.lat, - LMAX=LMAX, METHOD="integration") - YlmF = gravtk.gen_harmonics(ocean_function, landsea.lon, landsea.lat, - LMAX=LMAX, METHOD="fourier") + YlmI = gravtk.gen_harmonics( + ocean_function, + landsea.lon, + landsea.lat, + LMAX=LMAX, + METHOD='integration', + ) + YlmF = gravtk.gen_harmonics( + ocean_function, landsea.lon, landsea.lat, LMAX=LMAX, METHOD='fourier' + ) # check that amplitudes of harmonic data are nearly equal difference = YlmI.amplitude - YlmF.amplitude harmonic_eps = np.finfo(np.float16).eps diff --git a/test/test_time.py b/test/test_time.py index 13de331..c8524e4 100644 --- a/test/test_time.py +++ b/test/test_time.py @@ -1,5 +1,5 @@ #!/usr/bin/env python -u""" +""" test_time.py (01/2023) Verify time conversion and utility functions @@ -11,126 +11,140 @@ test date parser for cases when only a date and no units Written 12/2020 """ + import pytest import numpy as np import gravity_toolkit as gravtk + # parameterize calendar dates -@pytest.mark.parametrize("YEAR", np.random.randint(1992,2020,size=2)) -@pytest.mark.parametrize("MONTH", np.random.randint(1,13,size=2)) +@pytest.mark.parametrize('YEAR', np.random.randint(1992, 2020, size=2)) +@pytest.mark.parametrize('MONTH', np.random.randint(1, 13, size=2)) # PURPOSE: verify forward and backwards time conversions -def test_julian(YEAR,MONTH): +def test_julian(YEAR, MONTH): # days per month in a leap and a standard year # only difference is February (29 vs. 28) - dpm_leap = np.array([31,29,31,30,31,30,31,31,30,31,30,31]) - dpm_stnd = np.array([31,28,31,30,31,30,31,31,30,31,30,31]) - DPM = dpm_stnd if np.mod(YEAR,4) else dpm_leap - assert (np.sum(DPM) == gravtk.time.calendar_days(YEAR).sum()) + dpm_leap = np.array([31, 29, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31]) + dpm_stnd = np.array([31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31]) + DPM = dpm_stnd if np.mod(YEAR, 4) else dpm_leap + assert np.sum(DPM) == gravtk.time.calendar_days(YEAR).sum() # calculate Modified Julian Day (MJD) from calendar date - DAY = np.random.randint(1,DPM[MONTH-1]+1) - HOUR = np.random.randint(0,23+1) - MINUTE = np.random.randint(0,59+1) - SECOND = 60.0*np.random.random_sample(1) - MJD = gravtk.time.convert_calendar_dates(YEAR, MONTH, DAY, - hour=HOUR, minute=MINUTE, second=SECOND, - epoch=(1858,11,17,0,0,0)) + DAY = np.random.randint(1, DPM[MONTH - 1] + 1) + HOUR = np.random.randint(0, 23 + 1) + MINUTE = np.random.randint(0, 59 + 1) + SECOND = 60.0 * np.random.random_sample(1) + MJD = gravtk.time.convert_calendar_dates( + YEAR, + MONTH, + DAY, + hour=HOUR, + minute=MINUTE, + second=SECOND, + epoch=(1858, 11, 17, 0, 0, 0), + ) # convert MJD to calendar date JD = np.squeeze(MJD) + 2400000.5 - YY,MM,DD,HH,MN,SS = gravtk.time.convert_julian(JD, - format='tuple', astype=np.float64) + YY, MM, DD, HH, MN, SS = gravtk.time.convert_julian( + JD, format='tuple', astype=np.float64 + ) # assert dates eps = np.finfo(np.float16).eps - assert (YY == YEAR) - assert (MM == MONTH) - assert (DD == DAY) - assert (HH == HOUR) - assert (MN == MINUTE) - assert (np.abs(SS - SECOND) < eps) + assert YY == YEAR + assert MM == MONTH + assert DD == DAY + assert HH == HOUR + assert MN == MINUTE + assert np.abs(SS - SECOND) < eps + # parameterize calendar dates -@pytest.mark.parametrize("YEAR", np.random.randint(1992,2020,size=2)) -@pytest.mark.parametrize("MONTH", np.random.randint(1,13,size=2)) +@pytest.mark.parametrize('YEAR', np.random.randint(1992, 2020, size=2)) +@pytest.mark.parametrize('MONTH', np.random.randint(1, 13, size=2)) # PURPOSE: verify forward and backwards time conversions -def test_decimal_dates(YEAR,MONTH): +def test_decimal_dates(YEAR, MONTH): # days per month in a leap and a standard year # only difference is February (29 vs. 28) - dpm_leap = np.array([31,29,31,30,31,30,31,31,30,31,30,31]) - dpm_stnd = np.array([31,28,31,30,31,30,31,31,30,31,30,31]) - DPM = dpm_stnd if np.mod(YEAR,4) else dpm_leap - assert (np.sum(DPM) == gravtk.time.calendar_days(YEAR).sum()) + dpm_leap = np.array([31, 29, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31]) + dpm_stnd = np.array([31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31]) + DPM = dpm_stnd if np.mod(YEAR, 4) else dpm_leap + assert np.sum(DPM) == gravtk.time.calendar_days(YEAR).sum() # calculate Modified Julian Day (MJD) from calendar date - DAY = np.random.randint(1,DPM[MONTH-1]+1) - HOUR = np.random.randint(0,23+1) - MINUTE = np.random.randint(0,59+1) - SECOND = 60.0*np.random.random_sample(1) + DAY = np.random.randint(1, DPM[MONTH - 1] + 1) + HOUR = np.random.randint(0, 23 + 1) + MINUTE = np.random.randint(0, 59 + 1) + SECOND = 60.0 * np.random.random_sample(1) # calculate year-decimal time - tdec = gravtk.time.convert_calendar_decimal(YEAR, MONTH, day=DAY, - hour=HOUR, minute=MINUTE, second=SECOND) + tdec = gravtk.time.convert_calendar_decimal( + YEAR, MONTH, day=DAY, hour=HOUR, minute=MINUTE, second=SECOND + ) # day of the year 1 = Jan 1, 365 = Dec 31 (std) - day_temp = np.mod(tdec, 1)*np.sum(DPM) + day_temp = np.mod(tdec, 1) * np.sum(DPM) DofY = np.floor(day_temp) + 1 # cumulative sum of the calendar dates - day_cumulative = np.cumsum(np.concatenate(([0],DPM))) + 1 + day_cumulative = np.cumsum(np.concatenate(([0], DPM))) + 1 # finding which month date is in i = np.nonzero((DofY >= day_cumulative[0:-1]) & (DofY < day_cumulative[1:])) - month_range = np.arange(1,13) + month_range = np.arange(1, 13) month = month_range[i] # finding day of the month day = (DofY - day_cumulative[i]) + 1 # convert residuals into time (hour, minute and second) - hour_temp = np.mod(day_temp,1)*24.0 - minute_temp = np.mod(hour_temp,1)*60.0 - second = np.mod(minute_temp,1)*60.0 + hour_temp = np.mod(day_temp, 1) * 24.0 + minute_temp = np.mod(hour_temp, 1) * 60.0 + second = np.mod(minute_temp, 1) * 60.0 # assert dates eps = np.finfo(np.float16).eps - assert (np.floor(tdec) == YEAR) - assert (month == MONTH) - assert (day == DAY) - assert (np.floor(hour_temp) == HOUR) - assert (np.floor(minute_temp) == MINUTE) - assert (np.abs(second - SECOND) < eps) + assert np.floor(tdec) == YEAR + assert month == MONTH + assert day == DAY + assert np.floor(hour_temp) == HOUR + assert np.floor(minute_temp) == MINUTE + assert np.abs(second - SECOND) < eps + # PURPOSE: test UNIX time def test_unix_time(): # ATLAS Standard Data Epoch UNIX = gravtk.utilities.get_unix_time('2018-01-01 00:00:00') - assert (UNIX == 1514764800) + assert UNIX == 1514764800 # check UNIX time conversion with delta times - output_time = gravtk.time.convert_delta_time(UNIX, - epoch1=gravtk.time._unix_epoch, epoch2=(2018,1,1), - scale=1.0) - assert (output_time == 0) + output_time = gravtk.time.convert_delta_time( + UNIX, epoch1=gravtk.time._unix_epoch, epoch2=(2018, 1, 1), scale=1.0 + ) + assert output_time == 0 + # PURPOSE: test parsing time strings def test_parse_date_string(): # time string for Modified Julian Days time_string = 'days since 1858-11-17T00:00:00' - epoch,to_secs = gravtk.time.parse_date_string(time_string) + epoch, to_secs = gravtk.time.parse_date_string(time_string) # check the epoch and the time unit conversion factors - assert np.all(epoch == [1858,11,17,0,0,0]) - assert (to_secs == 86400.0) + assert np.all(epoch == [1858, 11, 17, 0, 0, 0]) + assert to_secs == 86400.0 # time string for ATLAS Standard Data Epoch time_string = 'seconds since 2018-01-01T00:00:00' - epoch,to_secs = gravtk.time.parse_date_string(time_string) + epoch, to_secs = gravtk.time.parse_date_string(time_string) # check the epoch and the time unit conversion factors - assert np.all(epoch == [2018,1,1,0,0,0]) - assert (to_secs == 1.0) + assert np.all(epoch == [2018, 1, 1, 0, 0, 0]) + assert to_secs == 1.0 # time string for unitless case time_string = '2000-01-01T12:00:00' - epoch,to_secs = gravtk.time.parse_date_string(time_string) + epoch, to_secs = gravtk.time.parse_date_string(time_string) # check the epoch and the time unit conversion factors - assert np.all(epoch == [2000,1,1,12,0,0]) - assert (to_secs == 0.0) + assert np.all(epoch == [2000, 1, 1, 12, 0, 0]) + assert to_secs == 0.0 # time string for unitless case with a time zone time_string = '2000-01-01T12:00:00.000-06:00' - epoch,to_secs = gravtk.time.parse_date_string(time_string) + epoch, to_secs = gravtk.time.parse_date_string(time_string) # check the epoch and the time unit conversion factors - assert np.all(epoch == [2000,1,1,18,0,0]) - assert (to_secs == 0.0) + assert np.all(epoch == [2000, 1, 1, 18, 0, 0]) + assert to_secs == 0.0 + # PURPOSE: test months adjustment for special cases # parameterize calendar dates -@pytest.mark.parametrize("PROC", ['CSR','GFZ','GSFC','JPL']) +@pytest.mark.parametrize('PROC', ['CSR', 'GFZ', 'GSFC', 'JPL']) def test_adjust_months(PROC): # The 'Special Months' (Nov 2011, Dec 2011 and April 2012) with # Accelerometer shutoffs make the relation between month number @@ -141,7 +155,7 @@ def test_adjust_months(PROC): # For GSFC: Oct 2018 (202) is centered in Nov 2018 (203) # dates with special months for each processing center - center_dates = dict(CSR=[],GFZ=[],GSFC=[],JPL=[]) + center_dates = dict(CSR=[], GFZ=[], GSFC=[], JPL=[]) # CSR dates to test (year-decimal, GRACE month) center_dates['CSR'].append([2011.62465753, 116]) center_dates['CSR'].append([2011.70821918, 117]) @@ -196,9 +210,9 @@ def test_adjust_months(PROC): center_dates['JPL'].append([2015.62465753, 164]) # get dates and months for center - tdec,months = np.transpose(center_dates[PROC]) + tdec, months = np.transpose(center_dates[PROC]) # GRACE/GRACE-FO months with duplicates - temp = np.array(12.0*(tdec-2002.0)+1,dtype='i') + temp = np.array(12.0 * (tdec - 2002.0) + 1, dtype='i') assert np.any(temp != months.astype('i')) # run months adjustment to fix special cases temp = gravtk.time.adjust_months(temp) diff --git a/test/test_units.py b/test/test_units.py index efcc805..45b2879 100644 --- a/test/test_units.py +++ b/test/test_units.py @@ -1,5 +1,5 @@ #!/usr/bin/env python -u""" +""" test_units.py (03/2023) Verify spherical harmonic and spatial unit factors @@ -8,25 +8,33 @@ include comparisons without including elastic deformation Written 01/2023 """ + import pytest import numpy as np import gravity_toolkit as gravtk + # PURPOSE: test spherical harmonic units -@pytest.mark.parametrize("LMAX", np.random.randint(60,696,size=1)) +@pytest.mark.parametrize('LMAX', np.random.randint(60, 696, size=1)) def test_harmonic_units(LMAX): # extract arrays of kl, hl, and ll Love Numbers - LOVE = gravtk.load_love_numbers(LMAX, - LOVE_NUMBERS=0, REFERENCE='CF', - FORMAT='class') + LOVE = gravtk.load_love_numbers( + LMAX, LOVE_NUMBERS=0, REFERENCE='CF', FORMAT='class' + ) factors = gravtk.units(lmax=LMAX).harmonic(*LOVE) # cmwe, centimeters water equivalent - cmwe = factors.rho_e*factors.rad_e*(2.0*factors.l+1.0)/(1.0+LOVE.kl)/3.0 + cmwe = ( + factors.rho_e + * factors.rad_e + * (2.0 * factors.l + 1.0) + / (1.0 + LOVE.kl) + / 3.0 + ) assert np.all(factors.get('cmwe') == factors.cmwe) assert np.all(factors.get(gravtk.units.bycode(1)) == factors.cmwe) assert np.all(cmwe == factors.cmwe) # mmGH, mm geoid height - mmGH = np.ones((LMAX + 1))*(10.0*factors.rad_e) + mmGH = np.ones((LMAX + 1)) * (10.0 * factors.rad_e) assert np.all(factors.get('mmGH') == factors.mmGH) assert np.all(factors.get(gravtk.units.bycode(2)) == factors.mmGH) assert np.all(mmGH == factors.mmGH) @@ -44,20 +52,21 @@ def test_harmonic_units(LMAX): assert np.all(factors.get(gravtk.units.bycode(6)) == factors.cmVCU) # cmwe, centimeters water equivalent without elastic deformation factors = gravtk.units(lmax=LMAX).harmonic(*LOVE, include_elastic=False) - cmwe = factors.rho_e*factors.rad_e*(2.0*factors.l+1.0)/3.0 + cmwe = factors.rho_e * factors.rad_e * (2.0 * factors.l + 1.0) / 3.0 assert np.all(factors.get('cmwe') == factors.cmwe) + # PURPOSE: test harmonic units with different love number formats -@pytest.mark.parametrize("LMAX", np.random.randint(60,696,size=1)) +@pytest.mark.parametrize('LMAX', np.random.randint(60, 696, size=1)) def test_harmonic_love_numbers(LMAX): # extract arrays of kl, hl, and ll Love Numbers - hl,kl,ll = gravtk.load_love_numbers(LMAX, - LOVE_NUMBERS=0, REFERENCE='CF', - FORMAT='tuple') - LOVE = gravtk.load_love_numbers(LMAX, - LOVE_NUMBERS=0, REFERENCE='CF', - FORMAT='class') - factors_tuple = gravtk.units(lmax=LMAX).harmonic(hl,kl,ll) + hl, kl, ll = gravtk.load_love_numbers( + LMAX, LOVE_NUMBERS=0, REFERENCE='CF', FORMAT='tuple' + ) + LOVE = gravtk.load_love_numbers( + LMAX, LOVE_NUMBERS=0, REFERENCE='CF', FORMAT='class' + ) + factors_tuple = gravtk.units(lmax=LMAX).harmonic(hl, kl, ll) factors_class = gravtk.units(lmax=LMAX).harmonic(*LOVE) # cmwe, centimeters water equivalent assert np.all(factors_tuple.get('cmwe') == factors_class.get('cmwe')) @@ -66,35 +75,51 @@ def test_harmonic_love_numbers(LMAX): # mbar, equivalent surface pressure assert np.all(factors_tuple.get('mbar') == factors_class.get('mbar')) + # PURPOSE: test spatial units -@pytest.mark.parametrize("LMAX", np.random.randint(60,696,size=1)) +@pytest.mark.parametrize('LMAX', np.random.randint(60, 696, size=1)) def test_spatial_units(LMAX): # extract arrays of kl, hl, and ll Love Numbers - LOVE = gravtk.load_love_numbers(LMAX, - LOVE_NUMBERS=0, REFERENCE='CF', - FORMAT='class') + LOVE = gravtk.load_love_numbers( + LMAX, LOVE_NUMBERS=0, REFERENCE='CF', FORMAT='class' + ) factors = gravtk.units(lmax=LMAX).spatial(*LOVE) # cmwe, centimeters water equivalent - cmwe = 3.0*(1.0+LOVE.kl)/(1.0+2.0*factors.l)/(4.0*np.pi*factors.rad_e*factors.rho_e) + cmwe = ( + 3.0 + * (1.0 + LOVE.kl) + / (1.0 + 2.0 * factors.l) + / (4.0 * np.pi * factors.rad_e * factors.rho_e) + ) assert np.all(factors.get('cmwe') == factors.cmwe) assert np.all(cmwe == factors.cmwe) # mmwe, millimeters water equivalent - mmwe = 3.0*(1.0+LOVE.kl)/(1.0+2.0*factors.l)/(40.0*np.pi*factors.rad_e*factors.rho_e) + mmwe = ( + 3.0 + * (1.0 + LOVE.kl) + / (1.0 + 2.0 * factors.l) + / (40.0 * np.pi * factors.rad_e * factors.rho_e) + ) assert np.all(factors.get('mmwe') == factors.mmwe) assert np.all(mmwe == factors.mmwe) # cmwe, centimeters water equivalent without elastic deformation factors = gravtk.units(lmax=LMAX).spatial(*LOVE, include_elastic=False) - cmwe = 3.0/(1.0+2.0*factors.l)/(4.0*np.pi*factors.rad_e*factors.rho_e) + cmwe = ( + 3.0 + / (1.0 + 2.0 * factors.l) + / (4.0 * np.pi * factors.rad_e * factors.rho_e) + ) assert np.all(factors.get('cmwe') == factors.cmwe) + # PURPOSE: test unit attributes def test_unit_attributes(): units_name, units_longname = gravtk.units.get_attributes('cmwe') - assert (units_name == 'cm') - assert (units_longname == 'Equivalent_Water_Thickness') + assert units_name == 'cm' + assert units_longname == 'Equivalent_Water_Thickness' units_name, units_longname = gravtk.units.get_attributes('microGal') - assert (units_name == u'\u03BCGal') - assert (units_longname == 'Gravitational_Undulation') + assert units_name == '\u03bcGal' + assert units_longname == 'Gravitational_Undulation' units_name, units_longname = gravtk.units.get_attributes('mVCU') - assert (units_name == 'meters') - assert (units_longname == 'Viscoelastic_Crustal_Uplift') + assert units_name == 'meters' + assert units_longname == 'Viscoelastic_Crustal_Uplift'