WIP: Add antarctic melt time series

[xylar]

This is a work-in-progress PR describing the addition of an analysis task for plotting time series of melt rates under individual ice shelves and groups of ice shelves.

[xylar]

@stephenprice, I rebased this branch onto add_sose_T_S after making some small fixes to the Edison config files. If you've already made your own fork of this branch, you might want to rebase or reset --hard to pick up my changes.

[stephenprice]

Will do @xylar. Thanks -

[stephenprice]

@xylar -- I did some spot checking here of places where either melt flux or average melt rate from our simulations roughly agrees with Rignot et al. One or the other is close in a few places, but never both, and this appears to be because the shelf areas based on the cellMasks are too large in all cases I've checked (by 30-50%). Overall, the total shelf area estimate for all of Antarctica is ~2.7x too large relative to Rignot et al. I can imagine how this could be if we are always treating a grid cell as floating if any part of it is partially floating. Is that the case? This problem will get better at higher resolution, but it may make comparing modeled and obs. melt rates difficult for our lower-res. simulations. It's possible there is some other bug here as well, but I thought I'd get your thoughts before digging too deep on this. Thanks.

[xylar]

Sorry, @stephenprice, I am missing a factor here. freshwaterFlux should be multiplied by landIceFraction, since many cells are only partly covered by land ice. I'll make a note of this here and fix it momentarily.

[xylar]

I noticed that the numbers seemed weird but I didn't realize how bad the problem was or what the cause was until your comment.

[xylar]

This should be changed to:

areaCell = dsRestart.landIceFraction.isel(Time=0)*dsRestart.areaCell

[stephenprice]

Ok, great! Well that was an easy fix (I didn't know landIceFraction even existed). I'll give this a try and report back.

[stephenprice]

At first glance this looks much better. I'll do some more careful checks on the reported areas later this morning.

[stephenprice]

With this change, the area of the indiv. shelves I checked are within ~5% and the total ice shelf area of the AIS is within ~7% (relative to the values reported in Rignot et al.). This seems reasonable to me considering the coarse res. of the 60to30 mesh. Also, all of the modeled shelf areas are biased too-small which also seems reasonable when you consider the fractal nature of coastlines.

[xylar]

Great! I don't recall the exact algorithm I used to get landIceFraction but there's definitely some smoothing involved. I think there's also something that says, if a cell has more than 50% grounded ice, it's grounded, so that would likely account for the bias toward less ice-shelf area. This should decrease with higher res but stay biased in the same direction.

[stephenprice]

@xylar -- FYI, I am making code updates and pushing them to my fork now, so you don't have to make these same changes on your fork anymore if you don't want to.

[stephenprice]

Ok, thanks. Yeah, I remember that Ross and FR did not actually look to bad in these runs. But there were still huge biases elsewhere.

[stephenprice]

@xylar -- I'd like your opinion on next steps in terms of comparing melt rate time series from indiv. shelves to observations. The easy thing to do would be to make lookup tables (dictionaries) of the Rignot et al. melt rates and uncertainties (ditto for Depoorter and other estimates) for each shelf and then use those for plotting against modeled values (either as a mean and unc. bounds on the melt time series plots, or as some kind of bar chart like you suggested from Timmermann and Hellmer).

The harder (and more costly) version would be to re-grid the obs. data onto the relevant ocean mesh and re-calculate the "observed" values that we are comparing with. The downside to that would be that 1) we're actually changing the reported obs. values that we are comparing with and 2) the reported unc. loose a bit of their meaning.

I think we'll want the the harder version at some point (e.g. for comparing spatial patterns of melt). But that is a different exercise with different goals than this one. Note also that for some obs. datesets (like Depoorter et al.) there isn't really an option to remap reported values onto our grids.

[xylar]

The easy thing to do would be to make lookup tables (dictionaries) of the Rignot et al. melt rates and uncertainties (ditto for Depoorter and other estimates) for each shelf and then use those for plotting against modeled values (either as a mean and unc. bounds on the melt time series plots, or as some kind of bar chart like you suggested from Timmermann and Hellmer).

Yes, this is the first step. I would put the Rignot and Depoorter values (best estimate and uncertainties) into a .csv file (a simple spreadsheet format) within the Melt observations directory. csv files are really easy to read into numpy arrays in python, and I'm sure it's easy to take one column as dictionary keys and another as values so you could make this a dictionary. I'm hesitant to hard-code a dictionary of observations into the code itself.

[xylar]

The harder (and more costly) version would be to re-grid the obs. data onto the relevant ocean mesh and re-calculate the "observed" values that we are comparing with. The downside to that would be that 1) we're actually changing the reported obs. values that we are comparing with and 2) the reported unc. loose a bit of their meaning.

I think we could do this if we want to cross-correlate melt patterns. It would be easy to add this to the analysis that performs plots of melt fields, particularly if a climatology instead of a time series is what's needed.

I think we'll want the the harder version at some point (e.g. for comparing spatial patterns of melt). But that is a different exercise with different goals than this one. Note also that for some obs. datesets (like Depoorter et al.) there isn't really an option to remap reported values onto our grids.

Yes, I think this is a different analysis than we're currently working on.

[stephenprice]

Thanks for the feedback on this @xylar. I'm going to start w/ a short csv file w/ only one or two shelves in it. Once that is working I can update to include more. If there are any examples of were you currently use .csv files (or other tabular values) as inputs for analysis, let me know. Otherwise, I'll see what I can find on my own.

[xylar]

Here's stackoverflow post that might help you get started with csv:
https://stackoverflow.com/questions/26903304/reading-data-from-a-csv-file-in-python
To construct the file name for the cvs observations, do something like the following:

from ..shared.io.utility import build_config_full_path
...
observationsDirectory = build_config_full_path(config, 'oceanObservations', 'meltSubdirectory')
obsFileName = '{}/Rignot_2013_melt_rates.csv'.format(observationsDirectory)

We haven't had any need for csv file so far, so there aren't any examples in MPAS-Analysis yet.

[xylar]

By the way, there's an Excel file as part of the supplementary material for Rignot 2013 that you can convert to csv. So no need to do that work yourself, just download the example. You might have to modify a few ice-shelf names but that shouldn't be hard.

[xylar]

Following up on the stackoverflow example, the code might look like:

observationsDirectory = build_config_full_path(config, 'oceanObservations', 'meltSubdirectory')
obsFileName = '{}/Rignot_2013_melt_rates.csv'.format(observationsDirectory)
obsFile = csv.reader(obsFileName)
rignotDict = {}
for line in obsFile:
    shelfName = line[1]
    meltFlux = line[2]
    melFluxUncertainty = line[3]
    meltRate = line[4]
    meltRateUncertainty = line[5]
    rignotDict[shelfName] = {'meltFlux': meltFlux,
                             'meltFluxUncertainty': meltFluxUncertainty,
                             'meltRate': meltRate,
                             'meltRateUncertainty': meltRateUncertainty}

Then, you could access the values later as:

print 'The melt flux from Thwaites is', rignotDict['Thwaites']['meltFlux']

Again, I hope that helps.

[xylar]

By the way, a few recommendations for actually writing python code:

• install and use spyder as your editor (conda install spyder)
• add another useful library conda install pyflakes
• Go to Tools > Preferences > Editor > Advanced Settings and change "Indentation Characters" to "4 Spaces"
• Go to Tools > Preferences > Editor > Code Introspection/Analysis and turn on "Real-time code analysis".

This will give you warnings (orange bullets on the left) when you're not following the so-called PEP8 conventions for code formatting. We try to follow that standard pretty strictly.

If this is a lot of bother, no worries. We can always do a "reformatting" commit at the end.

[stephenprice]

Thanks for all these tips. I'll let you know how it goes.

[stephenprice]

I've got the full Rignot csv file reading into a dictionary now in a little test script. And while the original xls file is going to require some massaging first for this to work in practice, it is at least working in principle. Thanks for the suggestions.

[xylar]

Sounds great!

[stephenprice]

The i/o portion of this is working now (at NERSC). I'll spare you the further play-by-play, but just so you know, further devel. is happening on my fork at: https://github.com/stephenprice/MPAS-Analysis/commits/add_antarctic_melt_time_series

[stephenprice]

@xylar - I've gotten to a point where I can pull the obs values for melt flux / melt rate (or other values from the csv file) into the figures, but I'm not sure what the right syntax is for plotting them. Right now, I'm just pushing them onto the axis label. Can you point me to an example for how this is done for some other time series? E.g., I know sea ice plots against obs. values in their timeseries plots.

I've pushed the code to my fork if you want to take a look. Right now, it only works for select ice shelves because of name conflicts between our regions and the Rignot file. That will take a bit of work to fix, as some names / areas / values will need to be combined in the Rignot file first. But for now, if you just specify a few select ice shelves ('Totten', 'Getz', 'Abbot', 'Amery') in your config file it will work.

Lastly, should we be moving this discussion away from this page and over to something associated w/ my fork?

[xylar]

@stephenprice, I'm going to close this PR. Why don't you do a work-in-progress PR like this one but from your branch and we'll continue the discussion there. We can make a note that discussion started here.

[stephenprice]

Ok, will do. I suppose this is just a normal PR but w/ the specification in the description that it is a work-in-progress (and add don't merge tag,etc)?

[xylar]

Yes, exactly. Start the title with WIP and add a "don't merge" tag.