Add antarctic melt time series

[stephenprice]

This pull request adds the ability to plot time series of submarine melt rates for individual ice shelves.

[stephenprice]

@xylar - FYI, my first few commit messages here were done incorrectly, so I repaired them and did a force push to correct.

[xylar]

This is replaces #228, which contains some discussion related to this PR.

[xylar]

I would do the conversion to float here, rather than down below, because you never use the string values.

[stephenprice]

I tried this previously in my little test script. For some reason, it doesn't like trying to convert to float at that point:

ValueError: could not convert string to float: B

In [269]: float()?
File "", line 1
float()?
^
SyntaxError: invalid syntax

[xylar]

That's very odd. Can you give some more details of what the exact line is that causes that error? There should be no difference between doing the cast to float here or later on. It's just cleaner to do it here.

[stephenprice]

If you grab the csv file from:
/global/project/projectdirs/acme/observations/Ocean/Melt/
(or set a path to it)

... and then just try w/ this little script in ipython you'll likely see the same complaint I'm seeing:

import csv
fn = open('Rignot_2013_melt_rates.csv', 'rU')
data = csv.reader(fn)
obsDict = {}
for line in data:
shelfName = line[0]
meltFlux = float( line[6] )
meltRate = float( line[8] )
obsDict[shelfName] = { 'meltFlux': meltFlux, 'meltRate': meltRate }

Ahhh wait ... I bet I know what is happening. It is trying to convert the first line of the data file to floats, and that first line is column headings. So I probably need to skip that first line, and then it will work. Let me try that. Thanks -

[xylar]

Ah, I'm pretty sure you're right and you need to skip the first row, say with:

First = True
for line in data:
    if First:
        First = False
        continue
    ...

[xylar]

Actually, you can probably just do:

for line in data[1:]:
   ...

[stephenprice]

I found another fix for this on stack overflow. It's working now w/ the float conversion at the data-read in point. Thanks -

[xylar]

A suggested function for plotting these observations (has not been debugged!):

def plot(self, config, modelValues, obsMean, obsUncertainty, 
         title, xlabel, ylabel, fileout):

    """
    Plots the list of time series data sets and stores the result in an image
    file.
    Parameters
    ----------
    config : instance of ConfigParser
        the configuration, containing a [plot] section with options that
        control plotting

    modelValues : xarray DataSet
        the data set to be plotted

    obsMean : float
        a single observed mean value to plot as a constant line

    obsUncertainty : float
        The observed uncertainty, to plot as a shaded rectangle around the mean 

    title : str
        the title of the plot

    xlabel, ylabel : str
        axis labels

    fileout : str
        the file name to be written

    Authors
    -------
    Xylar Asay-Davis, Stephen Price
    """

    # play with these as you need
    figsize = (15, 6)
    dpi = 300
    modelLineStyle = 'b-'
    modelLineWidth = 1.2
    obsColor = 'gray'
    obsLineWidth = 2

    maxXTicks = 20
    calendar = self.calendar

    axis_font = {'size': config.get('plot', 'axisFontSize')}
    title_font = {'size': config.get('plot', 'titleFontSize'),
                  'color': config.get('plot', 'titleFontColor'),
                  'weight': config.get('plot', 'titleFontWeight')}

    plt.figure(figsize=figsize, dpi=dpi)

    minDays = modelValues.Time.min()
    maxDays = modelValues.Time.max()
    plt.plot(modelValues.Time.values, modelValues.values, modelLineStyle,
             linewidth=modelLineWidth)

    ax = plt.gca()

    # this makes a "patch" with a single rectangular polygon, where the pairs are time and melt 
    # rate for the 4 corners, and "True" means it is a closed polygon:
    patches = [Polygon([[minDays, obsMean-obsUncertainty], [maxDays, obsMean-obsUncertainty],
                        [maxDays, obsMean+obsUncertainty], [minDays, obsMean+obsUncertainty]], 
                       closed=True)]

    # make the polygon gray and mostly transparent
    p = PatchCollection(patches, color=obsColor, alpha=0.4)
    # add it to the plot on the current axis
    ax.add_collection(p)

    # also plot a line
    plt.plot([minDays, maxDays], [obsMean, obsMean], color=obsColor, linewidth=obsLineWidth)

    # this will need to be imported from shared.plot.plotting
    plot_xtick_format(plt, calendar, minDays, maxDays, maxXTicks)

    if title is not None:
        plt.title(title, **title_font)
    if xlabel is not None:
        plt.xlabel(xlabel, **axis_font)
    if ylabel is not None:
        plt.ylabel(ylabel, **axis_font)
    if fileout is not None:
        plt.savefig(fileout, dpi=dpi, bbox_inches='tight', pad_inches=0.1)

    if not config.getboolean('plot', 'displayToScreen'):
        plt.close()

[milenaveneziani]

sorry the lack of imagination: what does WIP stand for?

[xylar]

WIP stands for work-in-progress.

[stephenprice]

@xylar -- with a few tweaks, I've got your local plotting script mostly working, e.g. ...

You are only seeing the single 'mean' value line there because I'm having trouble with the patch object for displaying the range as a shaded patch. This fig. was made by commenting out the section that adds that patch to the plot. The offending line (351) seems to be this one:

p = PatchCollection(patches, color=obsColor, alpha=0.4)

, giving the following error:

AttributeError: 'Polygon' object has no attribute 'get_transform'
ERROR: analysis task timeSeriesAntarcticMelt failed during run
Traceback (most recent call last):
File "./run_analysis.py", line 409, in
run_analysis(config, analyses, configFiles, args.subtask)
File "./run_analysis.py", line 249, in run_analysis
raise lastException
AttributeError: 'Polygon' object has no attribute 'get_transform'

For now, I am going to commit and push the sort-of working version to my fork w/ this section commented out. Let me know if you have any thoughts on the error.

[xylar]

@stephenprice, so this was a good guess but it's not the Polygon we need. Should be:

from matplotlib.patches import Polygon

I think that's what's causing your error.

[xylar]

If you take a look at the example I emailed you from POPSICLES, you can see which modules I import there.

[stephenprice]

Ok, that works now. Thanks! I will try to tidy up the csv file so that we can re-test applying to all of the ice shelves we are modeling.

[stephenprice]

@xylar -- getting back to this after a long hiatus. For some reason, it doesn't look like I added a new plot showing that the melt rate unc. can now be plotted as well (but this was done / working ~2 weeks ago). I cleaned up the names in the csv table, so now I should be able to do a test w/ many more ice shelves included in the analysis (I will try this tomorrow).

One potential issue is that any ice shelf region for which there is something like a "_1", "_2" etc. or "_East", "_West" will need to first be combined (e.g. "Ronne_1" and "Ronne_2" would just have to become "Ronne") for comparison to Rignot. Maybe this has already been done.

I'd also like to talk with you at some point about what needs to be done to extend this to the 30to10 mesh. For the advis. board meeting that Todd and I have in early Oct., it would be nice to show a few comparative figures from a few diff. simulations (e.g., highlighting that biases appear to be better in the higher-res. simulation, at least based on what we know so far).

[stephenprice]

Shoot ... hit "comment and close" by accident rather than "comment" (reopened now)

[milenaveneziani]

hit "comment and close" by accident rather than "comment"

no worries, this was a rite of passage mistake for me too on github..

[stephenprice]

It looks like this is largely working now. Some example figs. below. The good ...

and the bad ...

But whatever. It will be great to have this capability now. I will try to apply to @JeremyFyke's longer B-case run tomorrow.

One interesting thing is that the melt flux vs. melt rate plots don't quite line up the same w.r.t. the observations. E.g. ...

I suspect this is because our shelf areas aren't always close to the areas that Rignot et al. used. I'm not sure if there's a way to fix this, or if we should just plan to mainly look at melt flux rather than melt rate if we want a really accurate comparison.

I should note these are from the old AGU 2016 runs. The same analysis applied to the more recent runs looks a bit better. I'll post some figs for tomorrows meeting on this.

[xylar]

One potential issue is that any ice shelf region for which there is something like a "_1", "_2" etc. or "_East", "_West" will need to first be combined (e.g. "Ronne_1" and "Ronne_2" would just have to become "Ronne") for comparison to Rignot. Maybe this has already been done.

I thought I was combining these already. If not, I agree that there's work to do.

[xylar]

@stephenprice, so if you're using the correct mask files, you should not see any of the _1 and _2 regions:

cd /project/projectdirs/acme/mpas_analysis/region_masks/
ncdump -v regionNames oEC60to30v3wLI_iceShelfMasks.nc | tail -n 120

I see only "complete" ice shelves when I do this. Note, these were generated using the driver script driver_scripts/setup_ice_shelves.py. Are you using this mask file or one generated with this driver script? If so, you shouldn't see any of the numbered subregions.

[xylar]

I'd also like to talk with you at some point about what needs to be done to extend this to the 30to10 mesh. For the advis. board meeting that Todd and I have in early Oct., it would be nice to show a few comparative figures from a few diff. simulations (e.g., highlighting that biases appear to be better in the higher-res. simulation, at least based on what we know so far).

Nothing at all needs to be done. This is already supported (at least on Edison) because I generated the needed mask file: oRSS30to10wLI_iceShelfMasks.nc.

[xylar]

One interesting thing is that the melt flux vs. melt rate plots don't quite line up the same w.r.t. the observations.

That's pretty unavoidable.

I suspect this is because our shelf areas aren't always close to the areas that Rignot et al. used. I'm not sure if there's a way to fix this, or if we should just plan to mainly look at melt flux rather than melt rate if we want a really accurate comparison.

It would be useful to confirm this. Maybe it's worth thinking about how to plot our areas (times landIceFraction) vs. theirs to get an idea of which shelves are particularly bad and what we want to do about it (if anything). I tend to agree that melt fluxes are likely more important when areas disagree.

I should note these are from the old AGU 2016 runs. The same analysis applied to the more recent runs looks a bit better. I'll post some figs for tomorrows meeting on this.

I wonder if the too-low melt rates in the Amundsen region and under Filchner might be better in more recent runs. I'm really surprised, by the way, that Filchner is low. Melt climatology maps indicated that it would be way too high, if I recall. Fimbul and the other Dronning Maud Land ice shelves will probably still have large biases in the newer runs. Perhaps the coastal current is too strong there?

[stephenprice]

Hi @xylar. Responding in order to your comments above:

1. Yes, at NERSC, you've already got these regions combined. I posted this before finding where the ice shelf masks .nc files were. So this is not a problem right now. Mainly, I just wanted to make it clear this was a requirement for using the Rignot data (and presumably others as well). I also noted it as a requirement in a README file that I put together for the Rignot .csv file.

2. I see the 30to10 regions masks there as well. I'm going to try to copy those over to Theta today and see if I can apply this whole process to Mark's recent 30to10 run there.

3. Regarding shelf area differences, I included columns of the Rignot shelf areas in the .csv file for exactly this reason (e.g., to provide some fractional statistic for how our areas align with theirs). I'm not reading them in right now, but they are there if we want to read them in and use them for something like this.

4. W.r.t. melt rates in the ASE ... this IS a recent run, although only ~10 yrs out. This is the same set of runs for which we plotted / discussed these figures a few weeks ago:

https://acme-climate.atlassian.net/wiki/spaces/OCNICE/pages/146309188/2017-09-07+AT+NOON+Ocean+Ice+Special+Topics+ALCC+status+update

What we saw here -- very roughly -- was that melt rates were LOW on the ASE and too high in say, Queen Maud Land, and that as we went from 60to30 to 30to10, both of these biases got better (ASE melt rates went up, QML melt rates went down). This is why I want to redo this same set of plots for the 30to10 run sitting on theta.

[stephenprice]

@xylar -- I've been thinking that when we add a few more obs. datasets in here (e.g., Depoorter and Moholdt), it might be better to include the obs. and unc. as simple error bars. Here's an example for what that might look like (I just plotted the same set of errors 3 times, w/ the 2nd two instances divided by a scalar - just to get the idea across). I still need to look at the logic for iterating over multiple obs. data sets, but I don't think that should be too hard. Let me know what you think. Flagging @JeremyFyke for his thoughts as well.

[xylar]

@stephenprice, yes, I agree, that's a nice idea. Do you have a good algorithm for spacing out the "times" where the observations are plotted for variable numbers of observations? Probably not too hard but also not trivial.

[stephenprice]

@xylar -- no, haven't gone there yet. The spacing is not that hard to fiddle with. For the plot here, I just divided up the horiz. distance based on fractions of "maxDays - minDays". I think that would work fine for starters. Even if there's only one obs., you can plot it in the same place (say, 1/3 of the way across). If there are more, you add them in a bit farther along at some constant increment (I used 1/10th of maxDays - minDays here). The non-trivial logic would be deciding when you do / don't have another obs. to plot on there. I'll push what I've got so far so you can take a look.

[xylar]

@stephenprice, @JeremyFyke would like to be able to work with with this data in NetCDF form. It doesn't look like you're currently writing out the time series anywhere, correct? We could add that as either part of this PR or in a follow-up PR. What do you think?

[xylar]

I have now cleaned up the plotting and moved it to the shared module. I also cleaned up the code a little bit and added support for plotting ice shelves even when observations are not available. I believe this is ready to merge

[xylar]

@stephenprice, I'm going to merge this as soon as my test on Theta finishes (assuming everything goes okay).

One think that would be nice is to add Antarctica (and maybe also Filchner-Ronne and Ross) to the observations. With Antarctica, the question is presumably whether to use the "total surveyed" or "total Antarctica" values (there's not an easy way to do both unless we treat "Antarctica" as a special region). With FRIS and RIS, the question is how to combine the uncertainties. My suggestion would be just to add them for the fluxes but I don't know what to do for the melt rates. Anyway, this issue isn't standing in the way of this PR because we can update the observations at any time.

[stephenprice]

@xylar -- I've looked over your last 3 commits and this all looks great to me. Thanks for the clean-up on the obs. plotting side. I'm fine with you merging this if there are no other problems w/ your testing. Or, let me know and I can also merge later today.

[xylar]

Some last test results: everything seems to work on Theta. Since NERSC is down for maintenance today, I'm just going to upload a few example images here:

We may want to rethink the placement of the legend, but I don't think that's a high priority issue for this PR.