diff --git a/README.md b/README.md index e4bec4701..90dd17edf 100644 --- a/README.md +++ b/README.md @@ -5,7 +5,7 @@ Analysis for simulations produced with Model for Prediction Across Scales (MPAS) components and the Energy Exascale Earth System Model (E3SM), which used those components. -![sea surface temperature](docs/_static/sst_example.png) +![sea surface temperature](docs/users_guide/_static/sst_example.png) ## conda-forge @@ -127,7 +127,7 @@ for more details. 3. If you installed the `mpas-analysis` package, run: `mpas_analysis myrun.cfg`. This will read the configuration first from `mpas_analysis/default.cfg` and then replace that - configuraiton with any changes from from `myrun.cfg` + configuration with any changes from from `myrun.cfg` 4. If you want to run a subset of the analysis, you can either set the `generate` option under `[output]` in your config file or use the `--generate` flag on the command line. See the comments in @@ -178,7 +178,7 @@ Note: for older runs, mpas-seaice files will be named: * `mpascice.rst.0002-01-01_00000.nc` * `streams.cice` * `mpas-cice_in` - Also, for older runs mpaso-in will be named: + Also, for older runs `mpaso_in` will be named: * `mpas-o_in` @@ -221,13 +221,13 @@ If you are running from a git repo: 2. If using the `mpas-analysis` conda package, download the job script and/or sample config file from the [example configs directory](https://github.com/MPAS-Dev/MPAS-Analysis/tree/develop/configs). - 2. Modify the number of parallel tasks, the run name, the output directory + 3. Modify the number of parallel tasks, the run name, the output directory and the path to the config file for the run. - 3. Note: the number of parallel tasks can be anything between 1 and the + 4. Note: the number of parallel tasks can be anything between 1 and the number of analysis tasks to be performed. If there are more tasks than parallel tasks, later tasks will simply wait until earlier tasks have finished. - 4. Submit the job using the modified job script + 5. Submit the job using the modified job script diff --git a/docs/users_guide/tasks/climatologyMapAntarcticMelt.rst b/docs/users_guide/tasks/climatologyMapAntarcticMelt.rst index 90e1acc90..72d0f930c 100644 --- a/docs/users_guide/tasks/climatologyMapAntarcticMelt.rst +++ b/docs/users_guide/tasks/climatologyMapAntarcticMelt.rst @@ -4,7 +4,7 @@ climatologyMapAntarcticMelt =========================== An analysis task for comparison of Antarctic maps of melt rates against -observations from `Rignot et al. (2013)`_. +observations from `Adusumilli et al. (2020) `_. Component and Tags:: @@ -76,7 +76,7 @@ For more details, see: Observations ------------ -:ref:`rignot_melt` +:ref:`adusumilli_melt` Example Result -------------- @@ -84,5 +84,3 @@ Example Result .. image:: examples/ant_melt.png :width: 720 px :align: center - -.. _`Rignot et al. (2013)`: http://doi.org/10.1126/science.1235798 \ No newline at end of file diff --git a/docs/users_guide/tasks/timeSeriesAntarcticMelt.rst b/docs/users_guide/tasks/timeSeriesAntarcticMelt.rst index 263ee35d5..6dcef1795 100644 --- a/docs/users_guide/tasks/timeSeriesAntarcticMelt.rst +++ b/docs/users_guide/tasks/timeSeriesAntarcticMelt.rst @@ -4,7 +4,8 @@ timeSeriesAntarcticMelt ======================= An analysis task for plotting time series of mean melt rates per ice shelf or -Antarctic region along with observations from `Rignot et al. (2013)`_. +Antarctic region along with observations from `Rignot et al. (2013)`_ +and `Adusumilli et al. (2020) `_. Component and Tags:: @@ -102,7 +103,8 @@ Other Options Observations ------------ -:ref:`rignot_melt` +* :ref:`rignot_melt` +* :ref:`adusumilli_melt` Example Result -------------- diff --git a/mpas_analysis/docs/parse_quick_start.py b/mpas_analysis/docs/parse_quick_start.py index d388d3866..3134d0358 100644 --- a/mpas_analysis/docs/parse_quick_start.py +++ b/mpas_analysis/docs/parse_quick_start.py @@ -12,8 +12,8 @@ def build_quick_start(): replace = {'# MPAS-Analysis': '# Quick Start Guide\n', '[![Build Status]': '', '[![Documentation Status]': '', - '![sea surface temperature](docs/_static/sst_example.png)': - '![sea surface temperature](_static/sst_example.png)\n'} + '![sea surface temperature](docs/users_guide/_static/sst_example.png)': + '![sea surface temperature](users_guide/_static/sst_example.png)\n'} skip = [('## conda-forge', '## Installation')] outContent = '' diff --git a/mpas_analysis/obs/observational_datasets.xml b/mpas_analysis/obs/observational_datasets.xml index 481958b6c..24b9e4743 100755 --- a/mpas_analysis/obs/observational_datasets.xml +++ b/mpas_analysis/obs/observational_datasets.xml @@ -572,6 +572,70 @@ + + + Antarctic melt rates and fluxes + + + ocean + + + Melt rates and melt fluxes from Adusumilli et al. (2020) + + + [Data from: Interannual variations in meltwater input to the Southern Ocean from Antarctic ice shelves](https://doi.org/10.6075/J04Q7SHT) + + + Under copyright (US) + + Use: This work is available from the UC San Diego Library. This digital + copy of the work is intended to support research, teaching, and private + study. + + Constraint(s) on Use: This work is protected by the U.S. Copyright Law + (Title 17, U.S.C.). Use of this work beyond that allowed by "fair use" + or any license applied to this work requires written permission of the + copyright holder(s). Responsibility for obtaining permissions and any + use and distribution of this work rests exclusively with the user and + not the UC San Diego Library. Inquiries can be made to the UC San Diego + Library program having custody of the work. + + + [Adusumilli et al. (2020)](https://doi.org/10.1038/s41561-020-0616-z) + + + @article{Adusumilli2020, + title = {Interannual variations in meltwater input to the {Southern} {Ocean} from {Antarctic} ice shelves}, + copyright = {2020 The Author(s), under exclusive licence to Springer Nature Limited}, + issn = {1752-0908}, + url = {https://www.nature.com/articles/s41561-020-0616-z}, + doi = {10.1038/s41561-020-0616-z}, + journal = {Nature Geoscience}, + author = {Adusumilli, Susheel and Fricker, Helen Amanda and Medley,Brooke and Padman, Laurie and Siegfried, Matthew R.}, + month = aug, + year = {2020}, + note = {Publisher: Nature Publishing Group}, + pages = {1--5} + } + + + - http://library.ucsd.edu/dc/object/bb0448974g/_3_1.h5 + + + preprocess_observations/preprocess_adusumilli_melt.py + + + - climatologyMapAntarcticMelt + - timeSeriesAntarcticMelt + + + Ocean/Melt/Adusumilli + + + adusumilli_melt + + + HadISST Nino 3.4 Index @@ -1215,41 +1279,41 @@ ERA5 is the fifth generation ECMWF reanalysis for the global climate and weather for the past 4 to 7 decades. Currently data is available from 1950, with Climate Data Store entries for 1950-1978 (preliminary back extension) and from 1959 onwards (final release plus timely updates, this page). ERA5 replaces the ERA-Interim reanalysis. - - Reanalysis combines model data with observations from across the world into a globally complete and - consistent dataset using the laws of physics. This principle, called data assimilation, is based on + + Reanalysis combines model data with observations from across the world into a globally complete and + consistent dataset using the laws of physics. This principle, called data assimilation, is based on the method used by numerical weather prediction centres, where every so many hours (12 hours at ECMWF) - a previous forecast is combined with newly available observations in an optimal way to produce a new + a previous forecast is combined with newly available observations in an optimal way to produce a new best estimate of the state of the atmosphere, called analysis, from which an updated, improved forecast is issued. - Reanalysis works in the same way, but at reduced resolution to allow for the provision of a dataset + Reanalysis works in the same way, but at reduced resolution to allow for the provision of a dataset spanning back several decades. Reanalysis does not have the constraint of issuing timely forecasts, - so there is more time to collect observations, and when going further back in time, to allow for the + so there is more time to collect observations, and when going further back in time, to allow for the ingestion of improved versions of the original observations, which all benefit the quality of the reanalysis product. - + ERA5 provides hourly estimates for a large number of atmospheric, ocean-wave and land-surface quantities. An uncertainty estimate is sampled by an underlying 10-member ensemble at three-hourly intervals. Ensemble mean and spread have been pre-computed for convenience. - Such uncertainty estimates are closely related to the information content of the available + Such uncertainty estimates are closely related to the information content of the available observing system which has evolved considerably over time. They also indicate flow-dependent sensitive areas. To facilitate many climate applications, monthly-mean averages have been pre-calculated too, though monthly means are not available for the ensemble mean and spread. - + ERA5 is updated daily with a latency of about 5 days (monthly means are available around the 6th of each month). In case that serious flaws are detected in this early release (called ERA5T), this data could be different from the final release 2 to 3 months later. In case that this occurs users are notified. - + The data set presented here is a regridded subset of the full ERA5 data set on native resolution. It is online on spinning disk, which should ensure fast and easy access. It should satisfy the requirements for most common applications. - + An overview of all ERA5 datasets can be found in this [article](https://confluence.ecmwf.int/display/CKB/The+family+of+ERA5+datasets). Information on access to ERA5 data on native resolution is provided in these [guidelines](https://confluence.ecmwf.int/display/CKB/How+to+download+ERA5). - - Data has been regridded to a regular lat-lon grid of 0.25 degrees for the reanalysis and + + Data has been regridded to a regular lat-lon grid of 0.25 degrees for the reanalysis and 0.5 degrees for the uncertainty estimate (0.5 and 1 degree respectively for ocean waves). There are four main sub sets: hourly and monthly products, both on pressure levels (upper air fields) and single levels (atmospheric, ocean-wave and land surface quantities). - + The present entry is "ERA5 monthly mean data on single levels from 1959 to present". @@ -1280,7 +1344,7 @@ } - - https://cds.climate.copernicus.eu/cdsapp#!/dataset/reanalysis-era5-single-levels-monthly-means?tab=form + - https://cds.climate.copernicus.eu/cdsapp#!/dataset/reanalysis-era5-single-levels-monthly-means?tab=form @@ -1303,21 +1367,21 @@ ocean - The ESA Sea State Climate Change Initiative (CCI) project has produced global merged multi-sensor + The ESA Sea State Climate Change Initiative (CCI) project has produced global merged multi-sensor time-series of monthly gridded satellite altimeter significant wave height (referred to as Level 4 (L4) data) with a particular focus for use in climate studies. - + This dataset contains the Version 1.1 Remote Sensing Sea Surface Height product, gridded over a global - regular cylindrical projection (1°x1° resolution), averaging valid and good measurements from all + regular cylindrical projection (1°x1° resolution), averaging valid and good measurements from all available altimeters on a monthly basis (using the L2P products also available). These L4 products are meant for statistics and visualization. - + This first version of the Sea State CCI products is inherited from the GlobWave project, - building on experience and existing outputs. It extends and improves the GlobWave products, + building on experience and existing outputs. It extends and improves the GlobWave products, which were a post-processing over existing L2 altimeter agency products with additional filtering, - corrections and variables. A major improvement consists in a new denoised sea surface height + corrections and variables. A major improvement consists in a new denoised sea surface height variable using Empirical Mode Decomposition, which was used as input to these monthly statistical fields. - + The altimeter data used in the Sea State CCI dataset v1.1 come from multiple satellite missions spanning from 1991 to 2018 (ERS-1, ERS-2, Topex, Envisat, GFO, CryoSat-2, Jason-1, Jason-2, Jason-3, SARAL). Many altimeters are bi-frequency (Ku-C or Ku-S) and only measurements in Ku band were used, @@ -1350,7 +1414,7 @@ } - - https://data.ceda.ac.uk/neodc/esacci/sea_state/data/v1.1_release/l4/v1.1 + - https://data.ceda.ac.uk/neodc/esacci/sea_state/data/v1.1_release/l4/v1.1 diff --git a/mpas_analysis/ocean/climatology_map_antarctic_melt.py b/mpas_analysis/ocean/climatology_map_antarctic_melt.py index 309dcc0aa..d7e9c8aad 100644 --- a/mpas_analysis/ocean/climatology_map_antarctic_melt.py +++ b/mpas_analysis/ocean/climatology_map_antarctic_melt.py @@ -10,6 +10,8 @@ # https://raw.githubusercontent.com/MPAS-Dev/MPAS-Analysis/main/LICENSE import os import csv + +import numpy as np import xarray as xr import dask from multiprocessing.pool import ThreadPool @@ -57,7 +59,7 @@ def __init__(self, config, mpasClimatologyTask, regionMasksTask, regionMasksTask : ``ComputeRegionMasks`` A task for computing region masks - controlconfig : mpas_tools.config.MpasConfigParser + controlConfig : mpas_tools.config.MpasConfigParser Configuration options for a control run """ # Authors @@ -81,8 +83,8 @@ def __init__(self, config, mpasClimatologyTask, regionMasksTask, seasons = config.getexpression(sectionName, 'seasons') if len(seasons) == 0: - raise ValueError('config section {} does not contain valid list ' - 'of seasons'.format(sectionName)) + raise ValueError(f'config section {sectionName} does not contain ' + f'valid list of seasons') comparisonGridNames = config.getexpression(sectionName, 'comparisonGrids') @@ -98,8 +100,8 @@ def __init__(self, config, mpasClimatologyTask, regionMasksTask, self.add_subtask(tableSubtask) if len(comparisonGridNames) == 0: - raise ValueError('config section {} does not contain valid list ' - 'of comparison grids'.format(sectionName)) + raise ValueError(f'config section {sectionName} does not contain ' + f'valid list of comparison grids') # the variable 'timeMonthly_avg_landIceFreshwaterFlux' will be added to # mpasClimatologyTask along with the seasons. @@ -115,17 +117,30 @@ def __init__(self, config, mpasClimatologyTask, regionMasksTask, if controlConfig is None: refTitleLabel = \ - 'Observations (Rignot et al, 2013)' + 'Observations (Adusumilli et al, 2020)' observationsDirectory = build_obs_path( config, 'ocean', 'meltSubdirectory') + comparison_res = config.getfloat( + 'climatology', 'comparisonAntarcticStereoResolution') + + # the maximum available resolution that is not coarser than the + # comparison + avail_res = np.array([10., 4., 1.]) + valid = avail_res >= comparison_res + if np.count_nonzero(valid) == 0: + res = np.amin(avail_res) + else: + res = np.amax(avail_res[valid]) + obsFileName = \ - '{}/Rignot_2013_melt_rates_6000.0x6000.0km_10.0km_' \ - 'Antarctic_stereo.nc'.format(observationsDirectory) + f'{observationsDirectory}/Adusumilli/Adusumilli_2020_' \ + f'iceshelf_melt_rates_2010-2018_v0_6000x6000km_{res:g}km_' \ + f'Antarctic_stereo.20230504.nc' refFieldName = 'meltRate' - outFileLabel = 'meltRignot' - galleryName = 'Observations: Rignot et al. (2013)' + outFileLabel = 'meltAdusumilli' + galleryName = 'Observations: Adusumilli et al. (2020)' remapObservationsSubtask = RemapObservedAntarcticMeltClimatology( parentTask=self, seasons=seasons, fileName=obsFileName, @@ -138,7 +153,7 @@ def __init__(self, config, mpasClimatologyTask, regionMasksTask, remapObservationsSubtask = None controlRunName = controlConfig.get('runs', 'mainRunName') galleryName = None - refTitleLabel = 'Control: {}'.format(controlRunName) + refTitleLabel = f'Control: {controlRunName}' refFieldName = mpasFieldName outFileLabel = 'melt' @@ -329,7 +344,7 @@ def __init__(self, parentTask, mpasClimatologyTask, controlConfig, mpasClimatologyTask : ``MpasClimatologyTask`` The task that produced the climatology to be remapped and plotted - controlconfig : mpas_tools.config.MpasConfigParser + controlConfig : mpas_tools.config.MpasConfigParser Configuration options for a control run (if any) regionMasksTask : ``ComputeRegionMasks`` @@ -347,7 +362,7 @@ def __init__(self, parentTask, mpasClimatologyTask, controlConfig, tags = ['climatology', 'table'] if subtaskName is None: - subtaskName = 'table{}'.format(season) + subtaskName = f'table{season}' # call the constructor from the base class (AnalysisTask) super(AntarcticMeltTableSubtask, self).__init__( @@ -357,7 +372,6 @@ def __init__(self, parentTask, mpasClimatologyTask, controlConfig, componentName=parentTask.componentName, tags=tags) - config = parentTask.config self.season = season self.mpasClimatologyTask = mpasClimatologyTask self.controlConfig = controlConfig @@ -484,16 +498,17 @@ def run_task(self): regionNames = decode_strings(ds.regionNames) - outDirectory = '{}/antarcticMelt/'.format( - build_config_full_path(config, 'output', 'tablesSubdirectory')) + tableBase = build_config_full_path(config, 'output', + 'tablesSubdirectory') + outDirectory = f'{tableBase}/antarcticMelt/' try: os.makedirs(outDirectory) except OSError: pass - tableFileName = '{}/antarcticMeltRateTable_{}.csv'.format(outDirectory, - self.season) + tableFileName = \ + f'{outDirectory}/antarcticMeltRateTable_{self.season}.csv' with open(tableFileName, 'w', newline='') as csvfile: writer = csv.DictWriter(csvfile, fieldnames=fieldNames) @@ -501,15 +516,15 @@ def run_task(self): writer.writeheader() for index, regionName in enumerate(regionNames): row = {'Region': regionName, - 'Area': '{}'.format(ds.area[index].values), - mainRunName: '{}'.format(ds.meltRates[index].values)} + 'Area': f'{ds.area[index].values}', + mainRunName: f'{ds.meltRates[index].values}'} if dsControl is not None: row[controlRunName] = \ - '{}'.format(dsControl.meltRates[index].values) + f'{dsControl.meltRates[index].values}' writer.writerow(row) - tableFileName = '{}/antarcticMeltFluxTable_{}.csv'.format(outDirectory, - self.season) + tableFileName = \ + f'{outDirectory}/antarcticMeltFluxTable_{self.season}.csv' with open(tableFileName, 'w', newline='') as csvfile: writer = csv.DictWriter(csvfile, fieldnames=fieldNames) @@ -517,9 +532,9 @@ def run_task(self): writer.writeheader() for index, regionName in enumerate(regionNames): row = {'Region': regionName, - 'Area': '{}'.format(ds.area[index].values), - mainRunName: '{}'.format(ds.totalMeltFlux[index].values)} + 'Area': f'{ds.area[index].values}', + mainRunName: f'{ds.totalMeltFlux[index].values}'} if dsControl is not None: row[controlRunName] = \ - '{}'.format(dsControl.totalMeltFlux[index].values) + f'{dsControl.totalMeltFlux[index].values}' writer.writerow(row) diff --git a/mpas_analysis/ocean/time_series_antarctic_melt.py b/mpas_analysis/ocean/time_series_antarctic_melt.py index 1935c198d..8c172f7f5 100644 --- a/mpas_analysis/ocean/time_series_antarctic_melt.py +++ b/mpas_analysis/ocean/time_series_antarctic_melt.py @@ -15,6 +15,7 @@ import matplotlib.pyplot as plt from geometric_features import FeatureCollection, read_feature_collection +from geometric_features.aggregation import get_aggregator_by_name from mpas_analysis.shared.analysis_task import AnalysisTask @@ -57,7 +58,7 @@ def __init__(self, config, mpasTimeSeriesTask, regionMasksTask, regionMasksTask : ``ComputeRegionMasks`` A task for computing region masks - controlconfig : mpas_tools.config.MpasConfigParser, optional + controlConfig : mpas_tools.config.MpasConfigParser, optional Configuration options for a control run (if any) """ # Authors @@ -169,8 +170,7 @@ def __init__(self, parentTask, startYear, endYear, mpasTimeSeriesTask, taskName=parentTask.taskName, componentName=parentTask.componentName, tags=parentTask.tags, - subtaskName='computeMeltRates_{:04d}-{:04d}'.format(startYear, - endYear)) + subtaskName=f'computeMeltRates_{startYear:04d}-{endYear:04d}') self.mpasTimeSeriesTask = mpasTimeSeriesTask self.run_after(mpasTimeSeriesTask) @@ -182,8 +182,8 @@ def __init__(self, parentTask, startYear, endYear, mpasTimeSeriesTask, self.restartFileName = None self.startYear = startYear self.endYear = endYear - self.startDate = '{:04d}-01-01_00:00:00'.format(self.startYear) - self.endDate = '{:04d}-12-31_23:59:59'.format(self.endYear) + self.startDate = f'{self.startYear:04d}-01-01_00:00:00' + self.endDate = f'{self.endYear:04d}-12-31_23:59:59' self.variableList = \ ['timeMonthly_avg_landIceFreshwaterFlux'] @@ -246,16 +246,17 @@ def run_task(self): mpasTimeSeriesTask = self.mpasTimeSeriesTask config = self.config - outputDirectory = '{}/iceShelfFluxes/'.format( - build_config_full_path(config, 'output', 'timeseriesSubdirectory')) + timeSeriesBase = build_config_full_path(config, 'output', + 'timeseriesSubdirectory') + outputDirectory = f'{timeSeriesBase}/iceShelfFluxes/' try: os.makedirs(outputDirectory) except OSError: pass - outFileName = '{}/iceShelfFluxes_{:04d}-{:04d}.nc'.format( - outputDirectory, self.startYear, self.endYear) + outFileName = f'{outputDirectory}/iceShelfFluxes_' \ + f'{self.startYear:04d}-{self.endYear:04d}.nc' # Load data: inputFile = mpasTimeSeriesTask.outputFile @@ -271,14 +272,14 @@ def run_task(self): if numpy.all(dsOut.Time.values == dsIn.Time.values): return else: - self.logger.warning('File {} is incomplete. Deleting ' - 'it.'.format(outFileName)) + self.logger.warning(f'File {outFileName} is incomplete. ' + f'Deleting it.') os.remove(outFileName) except OSError: # something is potentially wrong with the file, so let's delete # it and try again - self.logger.warning('Problems reading file {}. Deleting ' - 'it.'.format(outFileName)) + self.logger.warning(f'Problems reading file {outFileName}. ' + f'Deleting it.') os.remove(outFileName) restartFileName = \ @@ -310,7 +311,7 @@ def run_task(self): datasets = [] nTime = dsIn.sizes['Time'] for tIndex in range(nTime): - self.logger.info(' {}/{}'.format(tIndex+1, nTime)) + self.logger.info(f' {tIndex + 1}/{nTime}') freshwaterFlux = \ dsIn.timeMonthly_avg_landIceFreshwaterFlux.isel(Time=tIndex) @@ -320,7 +321,7 @@ def run_task(self): totalMeltFluxes = numpy.zeros((nRegions,)) for regionIndex in range(nRegions): - self.logger.info(' {}'.format(regionNames[regionIndex])) + self.logger.info(f' {regionNames[regionIndex]}') cellMask = \ dsRegionMask.regionCellMasks.isel(nRegions=regionIndex) @@ -402,18 +403,18 @@ def run_task(self): # ------- # Xylar Asay-Davis - outputDirectory = '{}/iceShelfFluxes/'.format( - build_config_full_path(self.config, 'output', - 'timeseriesSubdirectory')) + timeSeriesBase = build_config_full_path(self.config, 'output', + 'timeseriesSubdirectory') + outputDirectory = f'{timeSeriesBase}/iceShelfFluxes/' - outFileName = '{}/iceShelfFluxes_{:04d}-{:04d}.nc'.format( - outputDirectory, self.startYears[0], self.endYears[-1]) + outFileName = f'{outputDirectory}/iceShelfFluxes_' \ + f'{self.startYears[0]:04d}-{self.endYears[-1]:04d}.nc' if not os.path.exists(outFileName): inFileNames = [] for startYear, endYear in zip(self.startYears, self.endYears): - inFileName = '{}/iceShelfFluxes_{:04d}-{:04d}.nc'.format( - outputDirectory, startYear, endYear) + inFileName = f'{outputDirectory}/iceShelfFluxes_' \ + f'{startYear:04d}-{endYear:04d}.nc' inFileNames.append(inFileName) ds = xarray.open_mfdataset(inFileNames, combine='nested', @@ -461,7 +462,7 @@ def __init__(self, parentTask, iceShelf, regionIndex, controlConfig): regionIndex : int The index into the dimension ``nRegions`` of the ice shelf to plot - controlconfig : mpas_tools.config.MpasConfigParser, optional + controlConfig : mpas_tools.config.MpasConfigParser, optional Configuration options for a control run (if any) """ # Authors @@ -474,7 +475,7 @@ def __init__(self, parentTask, iceShelf, regionIndex, controlConfig): taskName=parentTask.taskName, componentName=parentTask.componentName, tags=parentTask.tags, - subtaskName='plotMeltRates_{}'.format(iceShelf.replace(' ', '_'))) + subtaskName=f'plotMeltRates_{iceShelf.replace(" ", "_")}') self.iceShelfMasksFile = parentTask.iceShelfMasksFile self.iceShelf = iceShelf @@ -503,9 +504,9 @@ def setup_and_check(self): self.xmlFileNames = [] for prefix in ['melt_flux', 'melt_rate']: + iceShelfSuffix = self.iceShelf.replace(" ", "_") self.xmlFileNames.append( - '{}/{}_{}.xml'.format(self.plotsDirectory, prefix, - self.iceShelf.replace(' ', '_'))) + f'{self.plotsDirectory}/{prefix}_{iceShelfSuffix}.xml') return def run_task(self): @@ -516,8 +517,8 @@ def run_task(self): # ------- # Xylar Asay-Davis, Stephen Price - self.logger.info("\nPlotting Antarctic melt rate time series for " - "{}...".format(self.iceShelf)) + self.logger.info(f'\nPlotting Antarctic melt rate time series for ' + f'{self.iceShelf}...') self.logger.info(' Load melt rate data...') @@ -542,6 +543,10 @@ def run_task(self): refTotalMeltFlux, refMeltRates = \ self._load_ice_shelf_fluxes(self.controlConfig) + else: + controlRunName = None + refTotalMeltFlux = None + refMeltRates = None # Load observations from multiple files and put in dictionary based # on shelf key name @@ -554,8 +559,7 @@ def run_task(self): obsDict = {} # dict for storing dict of obs data for obsName in obsFileNameDict: - obsFileName = '{}/{}'.format(observationsDirectory, - obsFileNameDict[obsName]) + obsFileName = f'{observationsDirectory}/{obsFileNameDict[obsName]}' obsDict[obsName] = {} obsFile = csv.reader(open(obsFileName, 'r')) next(obsFile, None) # skip the header line @@ -579,8 +583,21 @@ def run_task(self): 'meltRate': meltRate, 'meltRateUncertainty': meltRateUncertainty} break - - # If areas from obs file used need to be converted from sq km to sq m + regionGroup = 'Ice Shelves' + _, prefix, date = get_aggregator_by_name(regionGroup) + + obsFileName = f'{observationsDirectory}/Adusumilli/Adusumilli_2020_' \ + f'iceshelf_melt_rates_2010-2018_v0.20230504.' \ + f'{prefix}{date}.nc' + with xarray.open_dataset(obsFileName) as ds_adusumilli: + region_names = [name.values for name in ds_adusumilli.regionNames] + index = region_names.index(self.iceShelf) + ds_shelf = ds_adusumilli.isel(nRegions=index) + obsDict['Adusumilli et al. (2020)'] = { + 'meltFlux': ds_shelf.totalMeltFlux.values, + 'meltFluxUncertainty': ds_shelf.meltFluxUncertainty.values, + 'meltRate': ds_shelf.meanMeltRate.values, + 'meltRateUncertainty': ds_shelf.meltRateUncertainty.values} mainRunName = config.get('runs', 'mainRunName') movingAveragePoints = config.getint('timeSeriesAntarcticMelt', @@ -610,22 +627,23 @@ def run_task(self): obsDict[obsName]['meltRateUncertainty']) else: # append NaN so this particular obs won't plot - self.logger.warning('{} observations not available for ' - '{}'.format(obsName, self.iceShelf)) + self.logger.warning(f'{obsName} observations not available ' + f'for {self.iceShelf}') obsMeltFlux.append(None) obsMeltFluxUnc.append(None) obsMeltRate.append(None) obsMeltRateUnc.append(None) title = self.iceShelf.replace('_', ' ') + suffix = self.iceShelf.replace(' ', '_') xLabel = 'Time (yr)' yLabel = 'Melt Flux (GT/yr)' timeSeries = totalMeltFlux.isel(nRegions=self.regionIndex) - filePrefix = 'melt_flux_{}'.format(self.iceShelf.replace(' ', '_')) - outFileName = '{}/{}.png'.format(self.plotsDirectory, filePrefix) + filePrefix = f'melt_flux_{suffix}' + outFileName = f'{self.plotsDirectory}/{filePrefix}.png' fields = [timeSeries] lineColors = [config.get('timeSeries', 'mainColor')] @@ -685,8 +703,8 @@ def run_task(self): savefig(outFileName, config) - caption = 'Running Mean of Total Melt Flux under Ice ' \ - 'Shelves in the {} Region'.format(title) + caption = f'Running Mean of Total Melt Flux under Ice Shelves in ' \ + f'the {title} Region' write_image_xml( config=config, filePrefix=filePrefix, @@ -704,8 +722,8 @@ def run_task(self): timeSeries = meltRates.isel(nRegions=self.regionIndex) - filePrefix = 'melt_rate_{}'.format(self.iceShelf.replace(' ', '_')) - outFileName = '{}/{}.png'.format(self.plotsDirectory, filePrefix) + filePrefix = f'melt_rate_{suffix}' + outFileName = f'{self.plotsDirectory}/{filePrefix}.png' fields = [timeSeries] lineColors = [config.get('timeSeries', 'mainColor')] @@ -750,8 +768,8 @@ def run_task(self): savefig(outFileName, config) - caption = 'Running Mean of Area-averaged Melt Rate under Ice ' \ - 'Shelves in the {} Region'.format(title) + caption = f'Running Mean of Area-averaged Melt Rate under Ice ' \ + f'Shelves in the {title} Region' write_image_xml( config=config, filePrefix=filePrefix, @@ -764,7 +782,8 @@ def run_task(self): imageDescription=caption, imageCaption=caption) - def _load_ice_shelf_fluxes(self, config): + @staticmethod + def _load_ice_shelf_fluxes(config): """ Reads melt flux time series and computes regional total melt flux and mean melt rate. @@ -773,14 +792,16 @@ def _load_ice_shelf_fluxes(self, config): # ------- # Xylar Asay-Davis - outputDirectory = '{}/iceShelfFluxes/'.format( - build_config_full_path(config, 'output', 'timeseriesSubdirectory')) + timeSeriesBase = build_config_full_path(config, 'output', + 'timeseriesSubdirectory') + + outputDirectory = f'{timeSeriesBase}/iceShelfFluxes/' startYear = config.getint('timeSeries', 'startYear') endYear = config.getint('timeSeries', 'endYear') - outFileName = '{}/iceShelfFluxes_{:04d}-{:04d}.nc'.format( - outputDirectory, startYear, endYear) + outFileName = f'{outputDirectory}/iceShelfFluxes_' \ + f'{startYear:04d}-{endYear:04d}.nc' dsOut = xarray.open_dataset(outFileName) return dsOut.totalMeltFlux, dsOut.meltRates diff --git a/preprocess_observations/preprocess_adusumilli_melt.py b/preprocess_observations/preprocess_adusumilli_melt.py new file mode 100755 index 000000000..a0bb7db05 --- /dev/null +++ b/preprocess_observations/preprocess_adusumilli_melt.py @@ -0,0 +1,328 @@ +#!/usr/bin/env python +import os +import shutil +import subprocess + +import h5py +import numpy as np +import pyproj +import xarray as xr +from geometric_features import FeatureCollection, GeometricFeatures +from geometric_features.aggregation import get_aggregator_by_name +from mpas_tools.cime.constants import constants as cime_constants +from pyremap import Remapper, ProjectionGridDescriptor +from pyremap.polar import get_antarctic_stereographic_projection + +from mpas_analysis.shared.io.download import download_files +from mpas_analysis.shared.constants import constants as mpas_constants + + +def download_adusumilli(out_filename): + """ + Remap the Adusumilli et al. (2020) melt rates at 1 km resolution to an MPAS + mesh + + Parameters + ---------- + out_filename : str + The original Adusumilli et al. (2020) melt rates + """ + + if os.path.exists(out_filename): + return + + download_files(fileList=['_3_1.h5'], + urlBase='http://library.ucsd.edu/dc/object/bb0448974g/', + outDir='.') + shutil.move('_3_1.h5', out_filename) + + +def adusumilli_hdf5_to_netcdf(in_filename, out_filename): + """ + Convert Adusumilli et al. (2020) melt rates to NetCDF format + + Parameters + ---------- + in_filename : str + The original Adusumilli et al. (2020) melt rates + + out_filename : str + The Adusumilli et al. (2020) melt rates in NetCDF format + """ + if os.path.exists(out_filename): + return + + print(f'Reading {in_filename}...') + h5_data = h5py.File(in_filename, 'r') + + x = np.array(h5_data['/x'])[:, 0] + y = np.array(h5_data['/y'])[:, 0] + melt_rate = np.array(h5_data['/w_b']) + melt_rate_uncertainty = np.array(h5_data['/w_b_uncert']) + print('done.') + + print('Creating an xarray dataset...') + ds = xr.Dataset() + + projection = pyproj.Proj('+proj=stere +lat_ts=-71.0 +lat_0=-90 +lon_0=0.0 ' + '+k_0=1.0 +x_0=0.0 +y_0=0.0 +ellps=WGS84') + latlon_projection = pyproj.Proj(proj='latlong', datum='WGS84') + + print('Computing lat/lon...') + x_2d, y_2d = np.meshgrid(x, y) + transformer = pyproj.Transformer.from_proj(projection, + latlon_projection) + lon, lat = transformer.transform(x_2d, y_2d) + print('done.') + + # Adusumilli et al. (2020) ice density (caption of Fig. 1 and Methods + # section) + rho_ice = 917. + rho_fw = cime_constants['SHR_CONST_RHOFW'] + ice_to_fw_equiv = rho_ice / rho_fw + + ds['x'] = (('x',), x) + ds['y'] = (('y',), y) + ds['lon'] = (('y', 'x'), lon) + ds.lon.attrs['units'] = 'degrees' + ds['lat'] = (('y', 'x'), lat) + ds.lat.attrs['units'] = 'degrees' + ds['meltRate'] = (('y', 'x'), ice_to_fw_equiv * melt_rate) + ds.meltRate.attrs['units'] = 'm/yr of freshwater' + ds['meltRateUncertainty'] = (('y', 'x'), + ice_to_fw_equiv * melt_rate_uncertainty) + ds.meltRateUncertainty.attrs['units'] = 'm/yr of freshwater' + print('Writing the dataset...') + ds.to_netcdf(out_filename) + print('done.') + + +def remap_adusumilli(in_filename, out_prefix, date, task_count=512): + """ + Remap Adusumilli et al. (2020) melt rates to 10 km comparison grid + + Parameters + ---------- + in_filename : str + The Adusumilli et al. (2020) melt rates in NetCDF format + + out_prefix : str + A prefix for the file to contain the Adusumilli et al. (2020) melt + rates and melt fluxes remapped to the comparison grid + + date : str + A date string to append to the file name. + + task_count : int + The number of MPI tasks to use to create the mapping file + """ + ds = xr.open_dataset(in_filename) + + melt_attrs = ds.meltRate.attrs + uncert_attrs = ds.meltRateUncertainty.attrs + + mask = ds.meltRate.notnull() + ds['meltRate'] = ds.meltRate.where(mask, 0.) + ds['meltMask'] = mask.astype(float) + mask = ds.meltRateUncertainty.notnull() + ds['meltRateUncertSqr'] = (ds.meltRateUncertainty**2).where(mask, 0.) + ds['uncertMask'] = mask.astype(float) + ds = ds.drop_vars(['lat', 'lon', 'meltRateUncertainty']) + + in_x = ds.x.values + in_y = ds.y.values + lx = np.abs(1e-3 * (in_x[-1] - in_x[0])) + ly = np.abs(1e-3 * (in_y[-1] - in_y[0])) + dx = np.abs(1e-3 * (in_x[1] - in_x[0])) + + in_grid_name = f'{lx:g}x{ly:g}km_{dx:g}km_Antarctic_stereo' + + in_projection = pyproj.Proj('+proj=stere +lat_ts=-71.0 +lat_0=-90 ' + '+lon_0=0.0 +k_0=1.0 +x_0=0.0 +y_0=0.0 ' + '+ellps=WGS84') + + in_descriptor = ProjectionGridDescriptor.create( + in_projection, in_x, in_y, in_grid_name) + + width = 6000. + reses = [1., 4., 10.] + + for res in reses: + x_max = 0.5 * width * 1e3 + nx = int(width / res) + 1 + out_x = np.linspace(-x_max, x_max, nx) + + out_grid_name = f'{width:g}x{width:g}km_{res:g}km_Antarctic_stereo' + + out_projection = get_antarctic_stereographic_projection() + + out_descriptor = ProjectionGridDescriptor.create( + out_projection, out_x, out_x, out_grid_name) + + method = 'conserve' + + map_filename = f'map_{in_grid_name}_to_{out_grid_name}_{method}.nc' + + remapper = Remapper(in_descriptor, out_descriptor, map_filename) + + if not os.path.exists(map_filename): + remapper.build_mapping_file(method=method, mpiTasks=task_count, + esmf_parallel_exec='srun') + + ds_out = remapper.remap(ds) + mask = ds_out.meltMask > 0. + ds_out['meltRate'] = ds_out.meltRate.where(mask) + ds_out.meltRate.attrs = melt_attrs + mask = ds_out.uncertMask > 0. + ds_out['meltRateUncertainty'] = \ + (np.sqrt(ds_out.meltRateUncertSqr)).where(mask) + ds_out.meltRateUncertainty.attrs = uncert_attrs + ds_out = ds_out.drop_vars(['meltRateUncertSqr']) + ds_out.to_netcdf(f'{out_prefix}_{out_grid_name}.{date}.nc') + + +def compute_regional_means(in_filename, out_prefix, date, chunk_size=50000, + core_count=128, multiprocessing_method='spawn'): + """ + Remap the Adusumilli et al. (2020) melt rates at 1 km resolution to an MPAS + mesh + + Parameters + ---------- + in_filename : str + The Adusumilli et al. (2020) melt rates in NetCDF format + + out_prefix : str + A prefix for the file to contain the Adusumilli et al. (2020) mean melt + rates and melt fluxes aggregated over ice-shelf regions + + date : str + A date string to append to ``out_prefix``. + + chunk_size : int, optional + The number of grid points that are processed in one operation when + creating region masks + + core_count : int, optional + The number of processes to use to compute masks. + + multiprocessing_method : str, optional + The multiprocessing method use for python mask creation ('fork', + 'spawn' or 'forkserver') + """ + region_group = 'Ice Shelves' + aggregation_function, prefix, region_date = \ + get_aggregator_by_name(region_group) + out_filename = f'{out_prefix}.{date}.{prefix}{region_date}.nc' + mask_files = _compute_masks(in_filename, aggregation_function, chunk_size, + core_count, multiprocessing_method) + + ds = xr.open_dataset(in_filename) + cell_area = ((ds.x.values[1] - ds.x.values[0]) * + (ds.y.values[1] - ds.y.values[0])) + + rho_fw = cime_constants['SHR_CONST_RHOFW'] + kg_per_gt = mpas_constants.kg_per_GT + gt_per_m3 = rho_fw / kg_per_gt + + ds_out = xr.Dataset() + mean_melt_rate = [] + total_melt_flux = [] + area_melt = [] + melt_rate_uncert = [] + melt_flux_uncert = [] + area_uncert = [] + for region_index, region_name in enumerate(mask_files.keys()): + mask_filename = mask_files[region_name] + ds_mask = xr.open_dataset(mask_filename) + region_mask = ds_mask.regionMasks.isel(nRegions=0) == 1 + + melt_mean, melt_area = _compute_mean_and_area( + ds.meltRate.values, region_mask, cell_area) + melt_total = gt_per_m3 * melt_area * melt_mean + mean_melt_rate.append(melt_mean) + area_melt.append(melt_area) + total_melt_flux.append(melt_total) + uncert_sqr = ds.meltRateUncertainty.values**2 + uncert_mean_sqr, uncert_area = _compute_mean_and_area( + uncert_sqr, region_mask, cell_area) + uncert_rms = np.sqrt(uncert_mean_sqr) + area_uncert.append(uncert_area) + uncert_total = gt_per_m3 * uncert_area * uncert_rms + melt_rate_uncert.append(uncert_rms) + melt_flux_uncert.append(uncert_total) + + ds_out['regionNames'] = ('nRegions', list(mask_files.keys())) + + ds_out['meanMeltRate'] = (('nRegions',), mean_melt_rate) + ds_out.meanMeltRate.attrs['units'] = 'm/yr of freshwater' + ds_out['meltRateUncertainty'] = (('nRegions',), melt_rate_uncert) + ds_out.meltRateUncertainty.attrs['units'] = 'm/yr of freshwater' + ds_out['totalMeltFlux'] = (('nRegions',), total_melt_flux) + ds_out.totalMeltFlux.attrs['units'] = 'GT/yr' + ds_out['meltFluxUncertainty'] = (('nRegions',), melt_flux_uncert) + ds_out.meltFluxUncertainty.attrs['units'] = 'GT/yr' + ds_out['meltArea'] = (('nRegions',), area_melt) + ds_out.meltArea.attrs['units'] = 'm^2' + ds_out['uncertaintyArea'] = (('nRegions',), area_uncert) + ds_out.uncertaintyArea.attrs['units'] = 'm^2' + ds_out.to_netcdf(out_filename) + + +def main(): + prefix = 'Adusumilli_2020_iceshelf_melt_rates_2010-2018_v0' + date = '20230504' + hdf5_filename = f'{prefix}.h5' + netcdf_filename = f'{prefix}.{date}.nc' + + download_adusumilli(hdf5_filename) + adusumilli_hdf5_to_netcdf(hdf5_filename, netcdf_filename) + remap_adusumilli(netcdf_filename, prefix, date) + compute_regional_means(netcdf_filename, prefix, date) + + +def _compute_masks(mesh_filename, aggregation_function, chunk_size, core_count, + multiprocessing_method): + gf = GeometricFeatures() + fc = aggregation_function(gf) + try: + os.makedirs('ice_shelf_masks') + except FileExistsError: + pass + files = {} + for feature in fc.features: + region_name = feature['properties']['name'] + region_prefix = region_name.replace(' ', '_') + region_mask_filename = f'ice_shelf_masks/{region_prefix}.nc' + files[region_name] = region_mask_filename + if os.path.exists(region_mask_filename): + continue + fc_region = FeatureCollection(features=[feature]) + geojson_filename = f'ice_shelf_masks/{region_prefix}.geojson' + fc_region.to_geojson(geojson_filename) + + args = ['compute_projection_region_masks', + '-i', mesh_filename, + '--lon', 'lon', + '--lat', 'lat', + '-g', geojson_filename, + '-o', region_mask_filename, + '--chunk_size', f'{chunk_size}', + '--process_count', f'{core_count}', + '--multiprocessing_method', f'{multiprocessing_method}', + '--show_progress'] + subprocess.run(args=args, check=True) + return files + + +def _compute_mean_and_area(field, region_mask, cell_area): + valid_melt = np.isfinite(field) + mask = np.logical_and(region_mask, valid_melt) + valid_count = np.count_nonzero(mask) + area = cell_area * valid_count + field_mean = np.sum(field[mask] * cell_area) / area + return field_mean, area + + +if __name__ == '__main__': + main()