Polar Data to Read

[rgrumbine]

For sea ice and polar verification I've been using data sets it seems METplus doesn't yet deal with

1. IABP/IPAB Arctic and Antarctic drifting buoy observations https://iabp.apl.uw.edu/data.html

• ASCII
• Prime use is to extract the buoy locations and use that to evaluate drift guidance
• Secondary (for me) is that some of these include surface met. packages to support weather model evaluation.

2. OSI-SAF sea ice concentrations https://osi-saf.eumetsat.int/products/osi-401-d (e.g.)

• Netcdf
• Polar Stereographic 10 km for N and S hemisphere
• Lat-lon values embedded in each day's file

3. NSIDC v4 Climate Data Record ftp://sidads.colorado.edu/pub/DATASETS/G02202_V4

• Polar Stereographic 25 km for N and S hemisphere
• Netcdf
• Lat-lon are separate auxiliary files to the daily fields

4. National Ice Center (NIC) high resolution ice edge analysis (global, daily)

• ASCII (sample below)

OSI-SAF and NSIDC are for evaluating model concentration field (threat scores, Heidke Skill score, integrated ice extent error, ...)
IABP/IPAB are for evaluating sea ice motion field
NIC is for evaluating the ice edge location

I have codes (C++ or python) to read all of these.

NIC sample:
NATICE ICE EDGE DATE 20230815
WGS84 DECIMAL DEGREES x 1000
LINE 70180N146420W 70299N146558W 70452N146990W
LINE 70452N146990W 70664N147622W 70725N147699W 70782N147940W
70844N148403W 70897N149115W 70986N150069W 71050N150244W
71113N150690W 71164N151598W 71193N152594W 71207N152969W
71192N153142W 71108N153143W 70987N152529W 70939N151894W
70856N151355W 70722N151128W 70662N151004W 70578N150418W
70565N150253W 70463N150000W
LINE 68824N179118W 68909N178487W 68911N178103W 68853N177783W
68780N177629W 68713N177504W 68652N177408W 68565N177250W
68508N177130W 68399N176847W 68330N176560W 68342N176302W
68493N176012W 68660N176071W 68850N176244W 68974N176402W
69087N176593W 69169N176766W 69366N176914W 69524N176756W
69652N176583W 69869N176414W 69979N176500W 70061N176538W
70159N176525W 70281N176544W 70406N176702W 70493N176831W
70602N176809W 70640N176551W 70831N175930W 70956N175838W
71026N175799W 71128N175598W 71170N175286W 71292N175014W
... continuing through all the line segments in the hemisphere

Thanks,
Bob Grumbine

[JohnHalleyGotway]

Hi @rgrumbine, thanks for starting up this discussion. I wanted to follow up on item 2. OSI-SAF sea ice concentrations.

I pulled a sample data file from their ftp site (ftp://osisaf.met.no/reprocessed/ice/conc-cont-reproc/v1p2). I see that this is a CF-compliant netCDF file:

ncdump -h ice_conc_nh_ease-125_cont-reproc_201908071200.nc
...
	:Conventions = "CF-1.4" ;

So MET should be able to read it directly, especially after support for Lambert Azimuthal Equal Area grids was added in MET version 11.1.0 via dtcenter/MET#2565. While testing with the plot_data_plane tool does produce output, the map data is noticeably absent. That indicates that MET is NOT parsing the grid information properly, which is confirmed by running at verbosity level 4.

MET-main_v11.1/bin/plot_data_plane \
ice_conc_nh_ease-125_cont-reproc_201908071200.nc plot.ps \
'name="ice_conc"; level="(0,*,*)";' -v 4
...
DEBUG 4: Grid information:
DEBUG 4:    Projection: Lambert Azimuthal Equal Area
DEBUG 4:    Nx: 849
DEBUG 4:    Ny: 849
DEBUG 4:    SpheroidName: Undefined
DEBUG 4:    Lat_LL: nan
DEBUG 4:    Lon_LL: nan
DEBUG 4:    Lat_Pole: 90.000
DEBUG 4:    Lon_Pole: -0.000

Here's the image it creates but its not correctly placed on the earth. The gray data we see is just missing data values over land. We should see black map outlines.

I wrote up a bugfix GitHub issue in the MET repository: dtcenter/MET#2651

[rgrumbine]

Sorry I didn't specify better. It is the polar stereographic grids that we want (if only to be using only one projection that isn't lat-lon) -- ice_conc_nh_polstere-100_multi_202301131200.nc for instance.

The ease grids are ... well, fundamentally the same info, but slightly lower resolution and included by OSISAF more for a technical reason that doesn't apply to us.

[JohnHalleyGotway]

@rgrumbine good news on that front. MET seems to handle those files just fine, as tested through the plot_data_plane tool:

MET-main_v11.1/bin/plot_data_plane \
ice_conc_nh_polstere-100_cont-reproc_201908071200.nc plot.ps \
'name="ice_conc"; level="(0,*,*)";'

It plots this data with the map outlines matching the contours of the data well.

The following log message is printed:

DEBUG 2: NcCfFile::get_grid_mapping_polar_stereographic() --> This is an ellipsoidal earth.

But we are in a transition from an internal custom grid library to the more robust and widely used Proj library.

I would like to keep dtcenter/MET#2651 as an open bug however since MET should be able to read that data on a lambert azimuthal equal area grid.

[rgrumbine]

Good news John!

I agree with keeping the bug open. Lambert is a pretty common projection even if not one of my favorites.

I assume that users can change color bars? And it looks like the plot is for raw ice concentrations, rather than the qc'd. We'll be able to specify which fields are used from the file, right?

[JohnHalleyGotway]

Yes, of course, you can select which 2D fields are extracted from the NetCDF files through the name and level settings.

Looking in ice_conc_nh_polstere-100_cont-reproc_201908071200.nc I see variables for:

	short ice_conc(time, yc, xc) ;
	short standard_error(time, yc, xc) ;
	short smearing_standard_error(time, yc, xc) ;
	short algorithm_standard_error(time, yc, xc) ;
	byte status_flag(time, yc, xc) ;

You can easily retrieve data from one of these by setting name = NetCDF variable name. But, if you want to apply some additional logic to use the ice_conc value, but only when the status_flag meets some criteria (for example), that's an additional step. Options are...

• write a python embedding script to read this data, apply whatever logic you'd like, and serve up the resulting field of data in memory
• run the gen_vx_mask tool instead to implement the data filtering logic you'd like

Yes, you certainly can specify alternative color tables for plot_data_plane, but it really isn't intended as a full featured plotting utility. Instead, it's meant as a diagnostic tool to confirm that MET is reading your input data properly and placing it on the correct spot on the globe. That's why I always start with it for new datasets.

Here's an example of plotting the status_flag variable where I request the level using the time value and use a different color table.

plot_data_plane ice_conc_nh_polstere-100_cont-reproc_201908071200.nc status.ps \
'name="status_flag"; level="(@20190808,*,*)";' \
-color_table share/met/colortables/NCL_colortables/BlueWhiteOrangeRed.ctable

I'll note that the timestamp in the data is 12 hours > timestamp in the filename. Hopefully that's as it should be. Here's the plot_data_plane usage statement if you'd like to see more.

[rgrumbine]

The ice concentration grids are for plus or minus 12 hours of 00z on the given date, so the timestamp vs. filename checks. I don't remember why I started doing that a decade earlier than OSI-SAF, and they might have been copying me. Or we both copied NSIDC or someone.

Glad to know that METplus works as I'd expected on color bars and variable reading.

[JohnHalleyGotway]

Next, following up on item 3. NSIDC v4 Climate Data Record.

I pulled this sample data file:
https://polarwatch.noaa.gov/erddap/files/nsidcG02202v4shmday/seaice_conc_monthly_sh_202208_f17_v04r00.nc

And see that it also follows the NetCDF CF-convention:

ncdump -h seaice_conc_monthly_sh_202208_f17_v04r00.nc
...
:Conventions = "CF-1.6, ACDD-1.3" ;

Testing using plot_data_plane:

MET-main_v11.1/bin/plot_data_plane \
seaice_conc_monthly_sh_202208_f17_v04r00.nc plot.ps \
'name="cdr_seaice_conc_monthly"; level="(0,*,*)";'

And that produces this runtime error:

ERROR  : NcCfFile::get_grid_mapping_polar_stereographic() --> Didn't find X dimension (projection_x_coordinate) in netCDF file.

However, I do clearly see projection_x_coordinate identified as the standard name of the xgrid variable:

	float xgrid(x) ;
		xgrid:valid_range = -3950000.f, 3950000.f ;
		xgrid:units = "meters" ;
		xgrid:long_name = "projection_grid_x_centers" ;
		xgrid:standard_name = "projection_x_coordinate" ;
		xgrid:axis = "X" ;

If I rename the xgrid and ygrid variables to x and y, respectively, that turns them into coordinate variables (i.e. the variable and dimension names match).

ncrename -v xgrid,x -v ygrid,y \
seaice_conc_monthly_sh_202208_f17_v04r00.nc \
seaice_conc_monthly_sh_202208_f17_v04r00_MOD.nc

Doing so DOES enable plot_data_plane to run but we have a similar problem parsing the grid specification correctly:

MET-main_v11.1/bin/plot_data_plane \
seaice_conc_monthly_sh_202208_f17_v04r00_MOD.nc plot.ps \
'name="cdr_seaice_conc_monthly"; level="(0,*,*)";' -v 4
...
DEBUG 4: Grid information:
DEBUG 4:    Projection: Stereographic
DEBUG 4:    Nx: 316
DEBUG 4:    Ny: 332
DEBUG 4:    IsNorthHemisphere: false
DEBUG 4:    Lon_orient: -180.000
DEBUG 4:    Bx: nan
DEBUG 4:    By: nan
DEBUG 4:    Alpha: inf

So this warrants another bugfix issue in MET to...

• Determine whether coordinate variables are REQUIRED by the CF-convention or only recommended. If that latter, generalize MET to handle non-coordinate variables defining the X and Y coordinates.
• Investigate and fix the parsing of this polar stereographic grid spec.

Please see issue dtcenter/MET#2652.

[rgrumbine]

Thank you.

[JohnHalleyGotway]

@rgrumbine, regarding item 1. IABP/IPAB Arctic and Antarctic drifting buoy observations

I took a look at the raw data provided at:
https://iabp.apl.uw.edu/WebData/

And I see that data is stored in a separate file for each buoy. So I understand that the changing lat/lon location through time provides relevant information. But it isn't immediately obvious to me how we'd include this in existing MET tools and use it to compare to model output.

I do see that 39967.dat, for example, contains pressure (BP) and temperature (Ts) data and can see how we'd use that in Point-Stat. Would you envision that we take a simple approach here?

• Enhance MET's ascii2nc tool to handle this type of data.
• Process as many input files as are passed in on the command line and just store the data in NetCDF format for Point-Stat to use.
• There wouldn't be any timing info in the output file names since the output file would contain data for all times available in the input.
• When it's processed by Point-Stat, most of the obs would be ignored since they won't occur very near in time to the forecast valid time... and that's a little inefficient.

Do you have any advice or requests for how the timing info should be handled? For example, should ascii2nc support a time window option... so we could store buoy data by day... or month... or year?

And do you have any thoughts about how to incorporate the changing lat/lon drift information into an algorithm in MET?

[rgrumbine]

For the point usage (Ta, P, ...) a lot of this is a matter of Met internals, which I don't know about. For me as a user,

• METplus transcribing from ASCII to Netcdf is fine.
• Time information is important to me, not least so as to be able to tell whether the model was 10 degrees wrong because it was just that wrong, versus being that wrong because it was 30 minutes off in timing a frontal passage.
• I'll suggest that current model output timeliness (1/3/6/24 hours, depending) shouldn't constrain the capability in Met. If observations are every 10 minutes, as some are here, well, some day models will be doing that.
• Storing buoy information by chunks (I'm agnostic as to frequency) makes sense.
• The total data volume of the IABP archive is small -- about 400 Mb. I agree about the potential inefficiency of storing 10 minute data when used only, say, daily. But in this case, the volume is only a few days' sea ice grids, where I'll be looking at 30+ years of grids from multiple centers.

n.b. I say 'I', but that isn't because I'm the only person who will be doing this, just that I am one of the people who will and I know my own thinking the best.

[rgrumbine]

The track usage is a different sort of creature than typical for METplus (afaik). But one that's not unique to sea ice. In particular, it would be useful for hurricane modelers to be able to do similar evaluations of forecast distance and direction of motion versus the observed.

The fundamental is that the sea ice drift model/output gives the starting time and location of a point and then the forecast locations it will be through forecast lead time, alternately, the drift distance and direction. If the point was exactly a buoy, life is simple. But in general, it won't be because sea ice forecasters need to know about the motion of the entire ice pack, not just the central portion that has buoys. So it is necessary to match up model starting points to buoy locations, for those points close enough to buoys to be evaluated.

The SIDFEx -- Sea ice dynamics forecast experiment -- has a fair amount on this kind of forecast and its evaluation (including code, in R) https://sidfex.polarprediction.net/

I have a document discussing which metrics, of 13, are most informative for sea ice drift -- https://polar.ncep.noaa.gov/mmab/papers/tn315/

I don't think this really answers your questions, but hope it puts us on a convergent iteration path.

[JohnHalleyGotway]

@rgrumbine regarding item 4. National Ice Center (NIC) high resolution ice edge analysis, I see that this is a collection of lat/lon locations describing the ice edge. As such, it's similar to an ESRI shapefile or a polyline mask file used by MET. Both ESRI shapefiles and ASCII polyline files can be read MET's gen_vx_mask tool and used to define which grid points are inside/outside of the enclosed region.

Any chance that this data is already available in shapefile format?

It strikes me that this ASCII format likely doesn't follow a widely adopted and generally used standard. But please correct me if I'm wrong. If not, I'd propose ingesting this data with Python. One possible workflow is:

1. Write a python script to read the file and reformat it as an ESRI shapefile or following MET's polyline format. A shapefile would be preferable since its more widely used.
2. We could run the custom user script wrapper in METplus to automate calls to this python script.
3. Next, run gen_vx_mask using this shapefile or polyline as input and define a 0/1 mask on your target grid - presumably you'd use your model output domain.
4. Next, run grid_stat (or other grid-to-grid tools in MET) to compare to model output and compute categorical stats, like the false alarms and misses we discussed in the meeting today.

Does this workflow sound right? Or are you hoping to use this data in a different way?

[rgrumbine]

Translating from this ASCII to a shapefile sounds reasonable. I've started learning about shapefiles and it should be relatively straightforward to do.

These line segments, though, probably can't be used as you describe as they don't necessarily close on land.

More importantly, the model grid and categorical stats is not how these data are intended to be used. For the categorical stats, we'll use the concentration grids.

Rather, the plan is to find the model's ice edge and compute the distance between the model's edge and the observed edge(s). This is already a figure of merit for some of our users.

This capability, of computing distance between two curves, should be of general interest -- model's front vs. analyzed, model's 540 dm line vs. observed, model freezing line vs. observed, ....