diff --git a/docs/users_guide/tasks/climatologyMapWoa.rst b/docs/users_guide/tasks/climatologyMapWoa.rst index 2202960dc..3a4127900 100644 --- a/docs/users_guide/tasks/climatologyMapWoa.rst +++ b/docs/users_guide/tasks/climatologyMapWoa.rst @@ -4,7 +4,7 @@ climatologyMapWoa ================= An analysis task for comparing potential temperature and salinity -at various depths against WOA18 climatology. +at various depths against WOA23 climatology. Component and Tags:: @@ -40,7 +40,7 @@ The following configuration options are available for this task:: [climatologyMapWoaTemperature] ## options related to plotting climatology maps of potential temperature ## at various levels, including the sea floor against control model results - ## and WOA18 climatological data + ## and WOA23 climatological data # colormap for model/observations colormapNameResult = RdYlBu_r @@ -66,8 +66,8 @@ The following configuration options are available for this task:: [climatologyMapWoaSalinity] ## options related to plotting climatology maps of salinity - ## at various levels, including the sea floor against control model results - ## and WOA18 climatological data + ## at various levels, including the sea floor against control model results + ## and WOA23 climatological data # colormap for model/observations colormapNameResult = haline @@ -98,14 +98,14 @@ For more details, see: The option ``depths`` is a list of (approximate) depths at which to sample the temperature and salinity fields. A value of ``'top'`` indicates the sea -surface. Note that, for the annual climatology, WOA18 data is available down -to 5500 m, whereas, for the seasonal or monthly climatologies, WOA18 data +surface. Note that, for the annual climatology, WOA23 data is available down +to 5500 m, whereas, for the seasonal or monthly climatologies, WOA23 data is only available down to 1500 m. Observations ------------ -:ref:`woa18_t_s` +:ref:`woa23_t_s` Example Result -------------- diff --git a/docs/users_guide/tasks/regionalTSDiagrams.rst b/docs/users_guide/tasks/regionalTSDiagrams.rst index c2d240c92..f6a2115b3 100644 --- a/docs/users_guide/tasks/regionalTSDiagrams.rst +++ b/docs/users_guide/tasks/regionalTSDiagrams.rst @@ -80,7 +80,7 @@ The following configuration options are available for this task: # volMax = 1e12 # Obserational data sets to compare against - obs = ['SOSE', 'WOA18'] + obs = ['SOSE', 'WOA23'] [TSDiagramsForOceanBasins] ## options related to plotting T/S diagrams of major ocean basins @@ -119,7 +119,7 @@ The following configuration options are available for this task: zmax = 0 # Obserational data sets to compare against - obs = ['WOA18'] + obs = ['WOA23'] Similar config sections are included for other region groups. @@ -193,11 +193,11 @@ default is to read them from geojson file. Observations ------------ The ``obs`` option contains a list of the names of observational data sets. -Currently, "SOSE" and "WOA18" are the only data sets available, but we +Currently, "SOSE" and "WOA23" are the only data sets available, but we anticipate adding several additional data sets in the near future. :ref:`sose` -:ref:`woa18_t_s` +:ref:`woa23_t_s` Other Config Options -------------------- diff --git a/docs/users_guide/tasks/timeSeriesOceanRegions.rst b/docs/users_guide/tasks/timeSeriesOceanRegions.rst index fe4c6b378..344b65085 100644 --- a/docs/users_guide/tasks/timeSeriesOceanRegions.rst +++ b/docs/users_guide/tasks/timeSeriesOceanRegions.rst @@ -58,7 +58,7 @@ The following configuration options are available for this task: # zmax = -400 # Observational data sets to compare against - obs = ['SOSE', 'WOA18'] + obs = ['SOSE', 'WOA23'] Region Groups ------------- @@ -121,12 +121,12 @@ Observations ------------ ``obs`` is a list of the observational data sets to plot as reference lines -(constant in time). Possible values are ``'SOSE'`` and ``'WOA18'``. An empty +(constant in time). Possible values are ``'SOSE'`` and ``'WOA23'``. An empty list can be provided if no observations should be plotted. :ref:`sose` -:ref:`woa18_t_s` +:ref:`woa23_t_s` Example Result -------------- diff --git a/mpas_analysis/default.cfg b/mpas_analysis/default.cfg index b453faec9..c21ca49b3 100755 --- a/mpas_analysis/default.cfg +++ b/mpas_analysis/default.cfg @@ -517,7 +517,7 @@ meltSubdirectory = Melt soseSubdirectory = SOSE sshSubdirectory = SSH argoSubdirectory = ARGO -woa18Subdirectory = WOA18 +woa23Subdirectory = WOA23 schmidtkoSubdirectory = Schmidtko woceSubdirectory = WOCE no3Subdirectory = BGC/NO3 @@ -1716,7 +1716,7 @@ zmin = -1000 zmax = 0 # Observational data sets to compare against -obs = ['WOA18'] +obs = ['WOA23'] [timeSeriesAntarcticRegions] @@ -1758,7 +1758,7 @@ variables = [{'name': 'temperature', # zmax = -400 # Observational data sets to compare against -obs = ['SOSE', 'WOA18'] +obs = ['SOSE', 'WOA23'] [timeSeriesTransport] @@ -1874,7 +1874,7 @@ normType = log # volMax = 1e12 # Obserational data sets to compare against -obs = ['SOSE', 'WOA18'] +obs = ['SOSE', 'WOA23'] [TSDiagramsForGreenlandRegions] ## options related to plotting T/S diagrams of Greenland regions @@ -1917,7 +1917,7 @@ zmax = 0 #volMax = 1e12 # Obserational data sets to compare against -obs = ['WOA18'] +obs = ['WOA23'] [TSDiagramsForArcticOceanRegions] ## options related to plotting T/S diagrams of Arctic regions @@ -1960,7 +1960,7 @@ zmax = 0 #volMax = 1e12 # Obserational data sets to compare against -obs = ['WOA18'] +obs = ['WOA23'] [TSDiagramsForOceanBasins] ## options related to plotting T/S diagrams of major ocean basins @@ -1999,7 +1999,7 @@ zmin = -1000 zmax = 0 # Obserational data sets to compare against -obs = ['WOA18'] +obs = ['WOA23'] [climatologyMapVel] @@ -2463,7 +2463,7 @@ shallowVsDeepColormapDepth = -1000 [climatologyMapWoaTemperatureShallow] ## options related to plotting climatology maps of potential temperature ## at various levels, including the sea floor against control model results -## and WOA18 climatological data +## and WOA23 climatological data # colormap for model/observations colormapNameResult = RdYlBu_r @@ -2490,7 +2490,7 @@ normArgsDifference = {'vmin': -5., 'vmax': 5.} [climatologyMapWoaTemperatureDeep] ## options related to plotting climatology maps of potential temperature ## at various levels, including the sea floor against control model results -## and WOA18 climatological data +## and WOA23 climatological data # colormap for model/observations colormapNameResult = RdYlBu_r @@ -2518,7 +2518,7 @@ normArgsDifference = {'vmin': -2., 'vmax': 2.} [climatologyMapWoaSalinityShallow] ## options related to plotting climatology maps of salinity ## at various levels, including the sea floor against control model results -## and WOA18 climatological data +## and WOA23 climatological data # colormap for model/observations colormapNameResult = haline @@ -2545,7 +2545,7 @@ normArgsDifference = {'vmin': -1.5, 'vmax': 1.5} [climatologyMapWoaSalinityDeep] ## options related to plotting climatology maps of salinity ## at various levels, including the sea floor against control model results -## and WOA18 climatological data +## and WOA23 climatological data # colormap for model/observations colormapNameResult = haline diff --git a/mpas_analysis/obs/observational_datasets.xml b/mpas_analysis/obs/observational_datasets.xml index a987d2e82..b94b32620 100755 --- a/mpas_analysis/obs/observational_datasets.xml +++ b/mpas_analysis/obs/observational_datasets.xml @@ -188,6 +188,73 @@ + + + WOA23 Temperature and Salinity Climatology + + + ocean + + + The World Ocean Atlas 2023 (WOA23) was released on February 14, 2024 and + includes temperature, salinity, oxygen (and parameters Apparent Oxygen + Utilization and percent oxygen saturation) and inorganic nutrients + (phosphate, silicate, and nitrate). WOA23 includes approximately 1.8 + million new oceanographic casts added to the WOD since WOA18’s release, + as well as renewed and updated quality controls. + + + [NOAA National Centers for Environmental Information (NCEI) website] + (https://www.ncei.noaa.gov/products/world-ocean-atlas) + + + These data are openly available to the public. Please acknowledge the + use of these data with the text given in the acknowledgment attribute. + + + - [Locarnini et al. (2024)](https://repository.library.noaa.gov/view/noaa/60599) + - [Reagan et al. (2024)](https://repository.library.noaa.gov/view/noaa/60600) + + + @manual{locarnini_world_nodate, + title = {World Ocean Atlas 2023, Volume 1: Temperature}, + url = {https://repository.library.noaa.gov/view/noaa/60599}, + shorttitle = {World Ocean Atlas 2023, Volume 1}, + year = {2024}, + author = {Locarnini, Ricardo A. and Mishonov, Alexey V. and Baranova, + Olga K. and Reagan, James R. and Boyer, Timothy P. and Seidov, Dan and + Wang, Zhankun and Garcia, Hernan E. and Bouchard, Courtney and Cross, + Scott L. and Paver, Christopher R. and Dukhovskoy, Dmitry}, + } + + @manual{reagan_world_nodate, + title = {World Ocean Atlas 2023, Volume 2: Salinity}, + url = {https://repository.library.noaa.gov/view/noaa/60600}, + shorttitle = {World Ocean Atlas 2023, Volume 2}, + author = {Reagan, James R. and Seidov, Dan and Wang, Zhankun and + Dukhovskoy, Dmitry and Boyer, Timothy P. and Locarnini, Ricardo A. and + Baranova, Olga K. and Mishonov, Alexey V. and Garcia, Hernan E. and + Bouchard, Courtney and Cross, Scott L. and Paver, Christopher R.}, + year = {2024}, + } + + + - https://www.ncei.noaa.gov/access/world-ocean-atlas-2023/ + + + - preprocess_observations/preprocess_woa23.py + + + - climatologyMapWoa + + + Ocean/WOA23 + + + woa23_t_s + + + AVISO Absolute Dynamic Topography diff --git a/mpas_analysis/ocean/climatology_map_woa.py b/mpas_analysis/ocean/climatology_map_woa.py index 6b554a5ba..e7e34c1d1 100644 --- a/mpas_analysis/ocean/climatology_map_woa.py +++ b/mpas_analysis/ocean/climatology_map_woa.py @@ -10,7 +10,7 @@ # https://raw.githubusercontent.com/MPAS-Dev/MPAS-Analysis/main/LICENSE """ Analysis tasks for comparing Polar (and global) climatology maps against -WOA18 climatological data. +WOA23 climatological data. """ # Authors # ------- @@ -34,7 +34,7 @@ class ClimatologyMapWoa(AnalysisTask): """ An analysis task for comparison of polar and global temperature and - salinity against WOA18 climatology fields + salinity against WOA23 climatology fields """ # Authors # ------- @@ -65,7 +65,7 @@ def __init__(self, config, mpasClimatologyTask, 'mpas': 'timeMonthly_avg_activeTracers_temperature', 'units': r'$\degree$C', 'titleName': 'Potential Temperature', - 'obsFieldName': 't_an'}, + 'obsFieldName': 'pt_an'}, {'prefix': 'salinity', 'mpas': 'timeMonthly_avg_activeTracers_salinity', 'units': r'PSU', @@ -130,87 +130,57 @@ def __init__(self, config, mpasClimatologyTask, comparisonGridNames=comparisonGridNames, iselValues=None) - remapAnnObsSubtask = None - remapMonObsSubtask = None + remapObsSubtask = None if controlConfig is None: - # Since we have a WOA18 annual climatology file and - # another file containing the 12 WOA18 monthly climatologies, - # do the remapping for each season separately - observationsDirectory = build_obs_path( - config, 'ocean', 'woa18Subdirectory') + config, 'ocean', 'woa23Subdirectory') refFieldNames = [field['obsFieldName'] for field in fields] - if 'ANN' in seasons: - obsFileName = \ - '{}/woa18_decav_04_TS_ann_20190829.nc'.format( - observationsDirectory) - remapAnnObsSubtask = RemapWoaClimatology( - parentTask=self, seasons=['ANN'], - fileName=obsFileName, - outFilePrefix='woa18_ann', - fieldNames=refFieldNames, - depths=depths, - comparisonGridNames=comparisonGridNames, - subtaskName='remapObservationsAnn') - self.add_subtask(remapAnnObsSubtask) - - seasonsMinusAnn = list(seasons) - if 'ANN' in seasonsMinusAnn: - seasonsMinusAnn.remove('ANN') - if len(seasonsMinusAnn) > 0: - obsFileName = \ - '{}/woa18_decav_04_TS_mon_20190829.nc'.format( - observationsDirectory) - remapMonObsSubtask = RemapWoaClimatology( - parentTask=self, seasons=seasonsMinusAnn, - fileName=obsFileName, - outFilePrefix='woa18_mon', - fieldNames=refFieldNames, - depths=depths, - comparisonGridNames=comparisonGridNames, - subtaskName='remapObservationsMon') - - self.add_subtask(remapMonObsSubtask) + obsFileName = f'{observationsDirectory}/' \ + f'woa23_decav_04_pt_s_mon_ann.20241101.nc' + + remapObsSubtask = RemapWoaClimatology( + parentTask=self, seasons=seasons, + fileName=obsFileName, + outFilePrefix='woa23', + fieldNames=refFieldNames, + depths=depths, + comparisonGridNames=comparisonGridNames, + subtaskName='remapObservationsMon') + + self.add_subtask(remapObsSubtask) for field in fields: fieldPrefix = field['prefix'] upperFieldPrefix = fieldPrefix[0].upper() + fieldPrefix[1:] - sectionName = '{}{}'.format(self.taskName, upperFieldPrefix) + sectionName = f'{self.taskName}{upperFieldPrefix}' if controlConfig is None: - refTitleLabel = 'WOA18 Climatology' + refTitleLabel = 'WOA23 Climatology' refFieldName = field['obsFieldName'] - outFileLabel = '{}WOA18'.format(fieldPrefix) - galleryName = 'WOA18 Climatology' + outFileLabel = f'{fieldPrefix}WOA23' + galleryName = 'WOA23 Climatology' diffTitleLabel = 'Model - Climatology' else: controlRunName = controlConfig.get('runs', 'mainRunName') - galleryName = 'Control: {}'.format(controlRunName) + galleryName = f'Control: {controlRunName}' refTitleLabel = galleryName refFieldName = field['mpas'] - outFileLabel = '{}WOA18'.format(fieldPrefix) + outFileLabel = f'{fieldPrefix}WOA23' diffTitleLabel = 'Main - Control' for comparisonGridName in comparisonGridNames: for depthIndex, depth in enumerate(depths): for season in seasons: - if season == 'ANN': - remapObsSubtask = remapAnnObsSubtask - else: - remapObsSubtask = remapMonObsSubtask - subtaskName = 'plot{}_{}_{}_depth_{}'.format( - upperFieldPrefix, - season, - comparisonGridName, - depth) + subtaskName = f'plot{upperFieldPrefix}_{season}_' \ + f'{comparisonGridName}_depth_{depth}' subtask = PlotClimatologyMapSubtask( parentTask=self, @@ -222,8 +192,8 @@ def __init__(self, config, mpasClimatologyTask, depth=depth, subtaskName=subtaskName) - configSectionName = 'climatologyMapWoa{}'.format( - upperFieldPrefix) + configSectionName = \ + f'climatologyMapWoa{upperFieldPrefix}' # if available, use a separate color map for shallow # and deep @@ -248,7 +218,7 @@ def __init__(self, config, mpasClimatologyTask, imageCaption=field['titleName'], galleryGroup=field['titleName'], groupSubtitle=None, - groupLink='{}_woa'.format(fieldPrefix), + groupLink=f'{fieldPrefix}_woa', galleryName=galleryName, configSectionName=configSectionName) @@ -334,7 +304,7 @@ def get_observation_descriptor(self, fileName): # ------- # Xylar Asay-Davis, Milena Veneziani - # Load WOA18 climatological data + # Load WOA23 climatological data dsObs = self.build_observational_dataset(fileName) # create a descriptor of the observation grid using Lat/Lon @@ -364,15 +334,15 @@ def build_observational_dataset(self, fileName): # ------- # Xylar Asay-Davis, Milena Veneziani - # Load WOA18 climatological data + # Load WOA23 climatological data dsObs = xr.open_dataset(fileName) - # Rename coordinates to be consistent with other datasets - if 'month' in dsObs.coords: - # need to add a dummy year coordinate - dsObs = dsObs.rename({'month': 'Time'}) - dsObs.coords['month'] = dsObs['Time'] - dsObs.coords['year'] = ('Time', np.ones(dsObs.sizes['Time'], int)) + # must have capital Time + dsObs = dsObs.rename({'time': 'Time'}) + # make sure month is a coord + dsObs = dsObs.set_coords('month') + # add a dummy year to the dataset + dsObs.coords['year'] = ('Time', np.ones(dsObs.sizes['Time'], int)) data_vars = {} for fieldName in self.fieldNames: @@ -390,7 +360,8 @@ def build_observational_dataset(self, fileName): field = xr.concat(slices, dim='depthSlice') data_vars[fieldName] = field - dsObs = xr.Dataset(data_vars=data_vars, - coords={'depthSlice': depthNames}) + coords = {'depthSlice': depthNames} + + dsObs = xr.Dataset(data_vars=data_vars, coords=coords) return dsObs diff --git a/mpas_analysis/ocean/regional_ts_diagrams.py b/mpas_analysis/ocean/regional_ts_diagrams.py index e7bab4449..10f516444 100644 --- a/mpas_analysis/ocean/regional_ts_diagrams.py +++ b/mpas_analysis/ocean/regional_ts_diagrams.py @@ -23,6 +23,7 @@ from multiprocessing.pool import ThreadPool from geometric_features import FeatureCollection, read_feature_collection +from mpas_tools.cime.constants import constants as cime_constants from mpas_analysis.shared.analysis_task import AnalysisTask @@ -98,6 +99,7 @@ def __init__(self, config, mpasClimatologyTask, regionMasksTask, self.seasons = config.getexpression(self.taskName, 'seasons') + woaFilename = 'WOA23/woa23_decav_04_pt_s_mon_ann.20241101.nc' obsDicts = { 'SOSE': { 'suffix': 'SOSE', @@ -120,29 +122,29 @@ def __init__(self, config, mpasClimatologyTask, regionMasksTask, 'volVar': 'volume', 'zVar': 'z', 'tVar': 'Time'}, - 'WOA18': { - 'suffix': 'WOA18', + 'WOA23': { + 'suffix': 'WOA23', 'gridName': 'Global_0.25x0.25degree', - 'gridFileName': 'WOA18/woa18_decav_04_TS_mon_20190829.nc', - 'TFileName': 'WOA18/woa18_decav_04_TS_mon_20190829.nc', - 'SFileName': 'WOA18/woa18_decav_04_TS_mon_20190829.nc', + 'gridFileName': woaFilename, + 'TFileName': woaFilename, + 'SFileName': woaFilename, 'volFileName': None, 'preprocessedFileTemplate': - 'WOA18/woa18_{}_T_S_z_vol_20200514.nc', + 'WOA23/woa23_{}_decav_04_pt_s_z_vol.20241101.nc', 'lonVar': 'lon', 'latVar': 'lat', - 'TVar': 't_an', + 'TVar': 'pt_an', 'SVar': 's_an', 'volVar': 'volume', 'zVar': 'depth', - 'tVar': 'month'}} + 'tVar': 'time'}} allObsUsed = [] for regionGroup in regionGroups: sectionSuffix = regionGroup[0].upper() + \ regionGroup[1:].replace(' ', '') - sectionName = 'TSDiagramsFor{}'.format(sectionSuffix) + sectionName = f'TSDiagramsFor{sectionSuffix}' obsList = config.getexpression(sectionName, 'obs') allObsUsed = allObsUsed + obsList allObsUsed = set(allObsUsed) @@ -161,7 +163,7 @@ def __init__(self, config, mpasClimatologyTask, regionMasksTask, for regionGroup in regionGroups: sectionSuffix = regionGroup[0].upper() + \ regionGroup[1:].replace(' ', '') - sectionName = 'TSDiagramsFor{}'.format(sectionSuffix) + sectionName = f'TSDiagramsFor{sectionSuffix}' regionNames = config.getexpression(sectionName, 'regionNames') if len(regionNames) == 0: @@ -290,7 +292,7 @@ def __init__(self, parentTask, obsName, obsDict, season): taskName=parentTask.taskName, componentName=parentTask.componentName, tags=parentTask.tags, - subtaskName='climatolgy{}_{}'.format(obsDict['suffix'], season)) + subtaskName=f'climatolgy{obsDict["suffix"]}_{season}') self.obsName = obsName self.obsDict = obsDict @@ -332,8 +334,8 @@ def run_task(self): with dask.config.set(schedular='threads', pool=ThreadPool(self.daskThreads)): - self.logger.info("\n computing T S climatogy for {}...".format( - self.obsName)) + self.logger.info( + f"\n computing T S climatogy for {self.obsName}...") chunk = {obsDict['tVar']: 6} @@ -347,19 +349,19 @@ def run_task(self): obsFileName = build_obs_path( config, component=self.componentName, relativePath=obsDict['TFileName']) - self.logger.info(' Reading from {}...'.format(obsFileName)) + self.logger.info(f' Reading from {obsFileName}...') ds = xarray.open_dataset(obsFileName, chunks=chunk) if obsDict['SFileName'] != obsDict['TFileName']: obsFileName = build_obs_path( config, component=self.componentName, relativePath=obsDict['SFileName']) - self.logger.info(' Reading from {}...'.format(obsFileName)) + self.logger.info(f' Reading from {obsFileName}...') dsS = xarray.open_dataset(obsFileName, chunks=chunk) ds[SVarName] = dsS[SVarName] if obsDict['volFileName'] is None: # compute volume from lat, lon, depth bounds - self.logger.info(' Computing volume...'.format(obsFileName)) + self.logger.info(' Computing volume...') latBndsName = ds[latVarName].attrs['bounds'] lonBndsName = ds[lonVarName].attrs['bounds'] zBndsName = ds[zVarName].attrs['bounds'] @@ -370,7 +372,7 @@ def run_task(self): dLon = numpy.deg2rad(lonBnds[:, 1] - lonBnds[:, 0]) lat = numpy.deg2rad(ds[latVarName]) dz = zBnds[:, 1] - zBnds[:, 0] - radius = 6378137.0 + radius = cime_constants['SHR_CONST_REARTH'] area = radius**2*numpy.cos(lat)*dLat*dLon ds[volVarName] = dz*area @@ -378,7 +380,7 @@ def run_task(self): obsFileName = build_obs_path( config, component=self.componentName, relativePath=obsDict['volFileName']) - self.logger.info(' Reading from {}...'.format(obsFileName)) + self.logger.info(f' Reading from {obsFileName}...') dsVol = xarray.open_dataset(obsFileName) ds[volVarName] = dsVol[volVarName] @@ -386,13 +388,13 @@ def run_task(self): temp_files = [temp_file_name] self.logger.info(' computing the dataset') ds.compute() - self.logger.info(' writing temp file {}...'.format(temp_file_name)) + self.logger.info(f' writing temp file {temp_file_name}...') write_netcdf_with_fill(ds, temp_file_name) chunk = {obsDict['latVar']: 400, obsDict['lonVar']: 400} - self.logger.info(' Reading back from {}...'.format(temp_file_name)) + self.logger.info(f' Reading back from {temp_file_name}...') ds = xarray.open_dataset(temp_file_name, chunks=chunk) if obsDict['tVar'] in ds.dims: @@ -413,9 +415,9 @@ def run_task(self): temp_files.append(temp_file_name) self.logger.info(' computing the dataset') ds.compute() - self.logger.info(' writing temp file {}...'.format(temp_file_name)) + self.logger.info(f' writing temp file {temp_file_name}...') write_netcdf_with_fill(ds, temp_file_name) - self.logger.info(' Reading back from {}...'.format(temp_file_name)) + self.logger.info(f' Reading back from {temp_file_name}...') ds = xarray.open_dataset(temp_file_name, chunks=chunk) self.logger.info(' Broadcasting z coordinate...') @@ -433,15 +435,15 @@ def run_task(self): self.logger.info(' computing the dataset') ds.compute() - self.logger.info(' writing {}...'.format(self.fileName)) + self.logger.info(f' writing {self.fileName}...') write_netcdf_with_fill(ds, self.fileName) for file_name in temp_files: - self.logger.info(' Deleting temp file {}'.format(file_name)) + self.logger.info(f' Deleting temp file {file_name}') os.remove(file_name) self.logger.info(' Done!') def _get_file_name(self, obsDict, suffix=''): - obsSection = '{}Observations'.format(self.componentName) + obsSection = f'{self.componentName}Observations' climatologyDirectory = build_config_full_path( config=self.config, section='output', relativePathOption='climatologySubdirectory', @@ -449,9 +451,8 @@ def _get_file_name(self, obsDict, suffix=''): make_directories(climatologyDirectory) - fileName = '{}/{}_{}_{}{}.nc'.format( - climatologyDirectory, 'TS_{}'.format(obsDict['suffix']), - obsDict['gridName'], self.season, suffix) + fileName = f'{climatologyDirectory}/TS_{obsDict["suffix"]}_' \ + f'{obsDict["gridName"]}_{self.season}{suffix}.nc' return fileName @@ -751,7 +752,7 @@ def _write_obs_t_s(self, obsDict, zmin, zmax): ds = ds.reset_index('nCells').drop_vars( [obsDict['latVar'], obsDict['lonVar']]) if 'nCells' in ds.data_vars: - ds = ds.drop_vars(['nCells']) + ds = ds.drop_vars(['nCells']) ds = ds.where(cellMask, drop=True) @@ -1077,7 +1078,7 @@ def run_task(self): else: norm = colors.LogNorm(vmin=volMinMpas, vmax=volMaxMpas) else: - raise ValueError('Unsupported normType {}'.format(normType)) + raise ValueError(f'Unsupported normType {normType}') if norm is not None: lastPanel.set_norm(norm) else: diff --git a/mpas_analysis/ocean/time_series_ocean_regions.py b/mpas_analysis/ocean/time_series_ocean_regions.py index fded1f505..5109bfa13 100644 --- a/mpas_analysis/ocean/time_series_ocean_regions.py +++ b/mpas_analysis/ocean/time_series_ocean_regions.py @@ -24,7 +24,7 @@ from mpas_analysis.shared.io import open_mpas_dataset, write_netcdf_with_fill from mpas_analysis.shared.io.utility import build_config_full_path, \ - build_obs_path, get_files_year_month, decode_strings, get_region_mask + build_obs_path, get_files_year_month, decode_strings from mpas_analysis.shared.html import write_image_xml @@ -74,6 +74,7 @@ def __init__(self, config, regionMasksTask, controlConfig=None): regionGroups = config.getexpression(self.taskName, 'regionGroups') + woaFilename = 'WOA23/woa23_decav_04_pt_s_mon_ann.20241101.nc' obsDicts = { 'SOSE': { 'suffix': 'SOSE', @@ -94,21 +95,21 @@ def __init__(self, config, regionMasksTask, controlConfig=None): 'zVar': 'z', 'tDim': 'Time', 'legend': 'SOSE 2005-2010 ANN mean'}, - 'WOA18': { - 'suffix': 'WOA18', + 'WOA23': { + 'suffix': 'WOA23', 'gridName': 'Global_0.25x0.25degree', - 'gridFileName': 'WOA18/woa18_decav_04_TS_mon_20190829.nc', - 'TFileName': 'WOA18/woa18_decav_04_TS_mon_20190829.nc', - 'SFileName': 'WOA18/woa18_decav_04_TS_mon_20190829.nc', + 'gridFileName': woaFilename, + 'TFileName': woaFilename, + 'SFileName': woaFilename, 'volFileName': None, 'lonVar': 'lon', 'latVar': 'lat', - 'TVar': 't_an', + 'TVar': 'pt_an', 'SVar': 's_an', 'volVar': 'volume', 'zVar': 'depth', - 'tDim': 'month', - 'legend': 'WOA18 1955-2017 ANN mean'}} + 'tDim': 'time', + 'legend': 'WOA23 1991-2020 ANN mean'}} for regionGroup in regionGroups: sectionSuffix = regionGroup[0].upper() + \ @@ -185,7 +186,8 @@ def __init__(self, config, regionMasksTask, controlConfig=None): localObsDict = dict(groupObsDicts[obsName]) obsSubtask = ComputeObsRegionalTimeSeriesSubtask( - self, regionGroup, regionName, fullSuffix, localObsDict) + self, regionGroup, regionName, fullSuffix, + localObsDict) obsSubtasks[obsName] = obsSubtask plotRegionSubtask = PlotRegionTimeSeriesSubtask( @@ -828,8 +830,7 @@ def run_task(self): # ------- # Xylar Asay-Davis - self.logger.info("\nAveraging T and S for {}...".format( - self.regionName)) + self.logger.info(f"\nAveraging T and S for {self.regionName}...") obsDict = self.obsDict config = self.config @@ -839,20 +840,19 @@ def run_task(self): sectionSuffix = regionGroup[0].upper() + \ regionGroup[1:].replace(' ', '') - sectionName = 'timeSeries{}'.format(sectionSuffix) + sectionName = f'timeSeries{sectionSuffix}' - outputDirectory = '{}/{}/'.format( - build_config_full_path(self.config, 'output', - 'timeseriesSubdirectory'), - timeSeriesName) + baseDir = build_config_full_path(self.config, 'output', + 'timeseriesSubdirectory') + outputDirectory = f'{baseDir}/{timeSeriesName}/' try: os.makedirs(outputDirectory) except OSError: pass - outFileName = '{}/TS_{}_{}.nc'.format( - outputDirectory, obsDict['suffix'], self.prefix) + outFileName = \ + f"{outputDirectory}/TS_{obsDict['suffix']}_{self.prefix}.nc" if os.path.exists(outFileName): return @@ -868,7 +868,7 @@ def run_task(self): dsRegionMask = dsRegionMask.reset_index('nCells').drop_vars( [obsDict['latVar'], obsDict['lonVar']]) if 'nCells' in dsRegionMask.data_vars: - dsRegionMaks = dsRegionMask.drop_vars(['nCells']) + dsRegionMask = dsRegionMask.drop_vars(['nCells']) maskRegionNames = decode_strings(dsRegionMask.regionNames) regionIndex = maskRegionNames.index(self.regionName) @@ -927,7 +927,7 @@ def run_task(self): if obsDict['volFileName'] is None: # compute volume from lat, lon, depth bounds - self.logger.info(' Computing volume...'.format(obsFileName)) + self.logger.info(' Computing volume...') latBndsName = ds[latVarName].attrs['bounds'] lonBndsName = ds[lonVarName].attrs['bounds'] zBndsName = ds[zVarName].attrs['bounds'] @@ -938,7 +938,7 @@ def run_task(self): dLon = numpy.deg2rad(lonBnds[:, 1] - lonBnds[:, 0]) lat = numpy.deg2rad(ds[latVarName]) dz = zBnds[:, 1] - zBnds[:, 0] - radius = 6378137.0 + radius = cime_constants['SHR_CONST_REARTH'] area = radius**2*numpy.cos(lat)*dLat*dLon volume = dz*area ds[volVarName] = volume diff --git a/mpas_analysis/polar_regions.cfg b/mpas_analysis/polar_regions.cfg index f0b05f03a..ce81fcb1e 100644 --- a/mpas_analysis/polar_regions.cfg +++ b/mpas_analysis/polar_regions.cfg @@ -22,7 +22,7 @@ depths = ['top', -50, -200, -400, -600, -800] [climatologyMapWoaTemperature] ## options related to plotting climatology maps of potential temperature ## at various levels, including the sea floor against control model results -## and WOA18 climatological data +## and WOA23 climatological data # A dictionary with keywords for the norm normArgsResult = {'vmin': -2., 'vmax': 2.} @@ -30,7 +30,7 @@ normArgsResult = {'vmin': -2., 'vmax': 2.} [climatologyMapWoaSalinity] ## options related to plotting climatology maps of salinity ## at various levels, including the sea floor against control model results -## and WOA18 climatological data +## and WOA23 climatological data # A dictionary with keywords for the norm normArgsResult = {'vmin': 33.8, 'vmax': 35.0} @@ -40,7 +40,7 @@ normArgsResult = {'vmin': 33.8, 'vmax': 35.0} ## options related to plotting T/S diagrams of ocean regions # the names of region groups to plot, each with its own section below -regionGroups = ['Ocean Basins', 'Arctic Ocean Regions', +regionGroups = ['Ocean Basins', 'Arctic Ocean Regions', 'Greenland Regions', 'Antarctic Regions'] [TSDiagramsForGreenlandRegions] diff --git a/mpas_analysis/shared/climatology/remap_observed_climatology_subtask.py b/mpas_analysis/shared/climatology/remap_observed_climatology_subtask.py index 162d1f9e0..3deebaab1 100644 --- a/mpas_analysis/shared/climatology/remap_observed_climatology_subtask.py +++ b/mpas_analysis/shared/climatology/remap_observed_climatology_subtask.py @@ -159,8 +159,7 @@ def run_task(self): stage='climatology', season=season, comparisonGridName=comparisonGridName) - if 'month' in ds.variables.keys() and \ - 'year' in ds.variables.keys(): + if 'month' in ds.coords and 'year' in ds.coords: # this data set is not yet a climatology, so compute # the climatology monthValues = constants.monthDictionary[season] @@ -171,7 +170,8 @@ def run_task(self): # climatology so assume this already is one seasonalClimatology = ds - write_netcdf_with_fill(seasonalClimatology, climatologyFileName) + write_netcdf_with_fill(seasonalClimatology, + climatologyFileName) remapper = self.remappers[comparisonGridName] diff --git a/preprocess_observations/preprocess_woa23.py b/preprocess_observations/preprocess_woa23.py new file mode 100755 index 000000000..6afccbdb3 --- /dev/null +++ b/preprocess_observations/preprocess_woa23.py @@ -0,0 +1,183 @@ +#!/usr/bin/env python +import os +from datetime import datetime + +import gsw +import numpy as np +import xarray as xr +from mpas_analysis.shared.climatology import compute_climatology +from mpas_analysis.shared.constants import constants +from mpas_analysis.shared.io.download import download_files +from mpas_analysis.shared.io import write_netcdf +from mpas_tools.cime.constants import constants as cime_constants + + +def download_woa(): + """ + Download 0.25 degree WOA23 annual and monthly files + """ + base_url = \ + 'https://www.ncei.noaa.gov/thredds-ocean/fileServer/woa23/DATA' + + for field in ['temperature', 'salinity']: + for index in range(13): + woa_filename = f'woa23_decav91C0_{field[0]}{index:02d}_04.nc' + woa_url = f'{base_url}/{field}/netcdf/decav91C0/0.25/' + + if os.path.exists(woa_filename): + continue + + download_files(fileList=[woa_filename], + urlBase=woa_url, + outDir='.') + + +def combine(): + """ + Run this step of the test case + """ + now = datetime.now() + datestring = now.strftime("%Y%m%d") + out_filename = f'woa23_decav_04_pt_s_mon_ann.{datestring}.nc' + if os.path.exists(out_filename): + return + + ds_ann = xr.open_dataset('woa23_decav91C0_t00_04.nc', decode_times=False) + + ds_out = xr.Dataset() + + for var in ['lon', 'lat', 'depth']: + ds_out[var] = ds_ann[var] + ds_out[f'{var}_bnds'] = ds_ann[f'{var}_bnds'] + + for field in ['temperature', 'salinity']: + var = f'{field[0]}_an' + ann_filename = f'woa23_decav91C0_{field[0]}00_04.nc' + ds_ann = xr.open_dataset( + ann_filename, decode_times=False).isel(time=0).drop_vars('time') + + time_slices = [] + for month in range(1, 13): + depth_slices = [] + ds_month = xr.open_dataset( + f'woa23_decav91C0_{field[0]}{month:02d}_04.nc', + decode_times=False).isel(time=0).drop_vars('time') + for depth_index in range(ds_ann.sizes['depth']): + if depth_index < ds_month.sizes['depth']: + ds = ds_month + else: + ds = ds_ann + depth_slices.append(ds[var].isel(depth=depth_index)) + + da_month = xr.concat(depth_slices, dim='depth') + print(month, da_month) + time_slices.append(da_month) + + ds_out[var] = xr.concat(time_slices, dim='time') + ds_out[var] = ds_out[var].transpose('time', 'lat', 'lon', 'depth') + ds_out[var].attrs = ds_ann[var].attrs + + ds_out['month'] = ('time', np.arange(1, 13)) + ds_out = _temp_to_pot_temp(ds_out) + ds_out.to_netcdf(out_filename) + + +def compute_obs_ts_climatology(season): + """ + season : str + The season over which to compute the climatology + """ + now = datetime.now() + datestring = now.strftime("%Y%m%d") + + in_filename = f'woa23_decav_04_pt_s_mon_ann.{datestring}.nc' + out_filename = f'woa23_{season}_decav_04_pt_s_z_vol.{datestring}.nc' + + if os.path.exists(out_filename): + return + + print("\n computing T S climatogy for WOA23...") + + print(f' Reading from {in_filename}...') + ds = xr.open_dataset(in_filename) + + # compute volume from lat, lon, depth bounds + print(' Computing volume...') + lat_bnds = ds['lat_bnds'] + lon_bnds = ds['lon_bnds'] + depth_bnds = ds['depth_bnds'] + dlat = np.deg2rad(lat_bnds[:, 1] - lat_bnds[:, 0]) + dlon = np.deg2rad(lon_bnds[:, 1] - lon_bnds[:, 0]) + lat = np.deg2rad(ds['lat']) + dz = depth_bnds[:, 1] - depth_bnds[:, 0] + radius = cime_constants['SHR_CONST_REARTH'] + area = radius**2*np.cos(lat)*dlat*dlon + ds['volume'] = dz*area + + ds = ds.rename({'time': 'Time'}) + ds = ds.set_coords('month') + ds.coords['year'] = np.ones(ds.sizes['Time'], int) + + month_values = constants.monthDictionary[season] + print(' Computing climatology...') + ds = compute_climatology(ds, month_values, maskVaries=True) + + print(' Broadcasting z coordinate...') + attrs = ds['depth'].attrs + ds['depth'] = -ds['depth'] + ds['depth'].attrs = attrs + ds['depth'].attrs['positive'] = 'up' + + _, _, z = xr.broadcast(ds['pt_an'], ds['s_an'], ds['depth']) + + ds['zBroadcast'] = z + + print(f' writing {out_filename}...') + write_netcdf.write_netcdf_with_fill(ds, out_filename) + print(' Done!') + + +def main(): + download_woa() + combine() + for season in ['ANN', 'JFM', 'JAS']: + compute_obs_ts_climatology(season) + + +def _temp_to_pot_temp(ds): + dims = ds.t_an.dims + + slices = list() + for depth_index in range(ds.sizes['depth']): + temp_slice = ds.t_an.isel(depth=depth_index) + in_situ_temp = temp_slice.values + salin = ds.s_an.isel(depth=depth_index).values + lat = ds.lat.broadcast_like(temp_slice).values + lon = ds.lon.broadcast_like(temp_slice).values + z = -ds.depth.isel(depth=depth_index).values + pressure = gsw.p_from_z(z, lat) + mask = np.isfinite(in_situ_temp) + SA = gsw.SA_from_SP(salin[mask], pressure[mask], lon[mask], + lat[mask]) + pot_temp = np.nan * np.ones(in_situ_temp.shape) + pot_temp[mask] = gsw.pt_from_t(SA, in_situ_temp[mask], + pressure[mask], p_ref=0.) + pot_temp_slice = xr.DataArray(data=pot_temp, dims=temp_slice.dims, + attrs=temp_slice.attrs) + + slices.append(pot_temp_slice) + + ds['pt_an'] = xr.concat(slices, dim='depth').transpose(*dims) + + ds.pt_an.attrs['standard_name'] = \ + 'sea_water_potential_temperature' + ds.pt_an.attrs['long_name'] = \ + 'Objectively analyzed mean fields for ' \ + 'sea_water_potential_temperature at standard depth levels.' + + ds = ds.drop_vars('t_an') + return ds + + +if __name__ == '__main__': + main()