Sea-Ice Extent Verification Using PointStat and python embedding

[malloryprow]

In one of the METplus User Telecons (Oct 3, 2022), I brought up how as a part of EVS we need to do verification for sea-ice extent (area considered ice covered). My initial approach involved using MODE and using the AREA output from the line type to calculate sea-ice extent, but as we touched on in that telecon the AREA is just the number of grid points, and we'd still need to calculate the area and keep track of the grid spacing to do so. The more I think about it this approach seems cumbersome.

The current RTOFS verification code uses python to calculate sea-ice extent. My new approach of thinking is that this pre-existing sea-ice extent calculation code can be adapted for the python embedding in METplus. The python embedding code will take the gridded input file and calculate sea-ice extent, which would result in a single value. This will need to be done for both the forecast and the observations. The value will get assigned to a lat,lon (the same one for both) and the PointStat can be used for getting stat file(s).

I haven't done python embedding before, but I know python well enough to know I can figure it out. Is there a good script to reference for using python embedding with PointStat?

[JohnHalleyGotway]

@malloryprow, there are three types of python embedding supported in MET.

1. To pass a single field of gridded data to any of the MET tools that read gridded data.
2. To pass point observations to any of the MET tools that read point observations.
3. To pass matched fcst/obs pairs to Stat-Analysis (described in this section).

Based on the situation you describe, I suspect that the 3rd option would work best. Use python compute the forecast ice area coverage and the verifying observation ice area coverage. Then format and pass the result as an MPR (matched pair) line type to Stat-Analysis. So you'd be calling Stat-Analysis instead of Point-Stat.

Honestly though, I don't know what piece of work you'd actually want Stat-Analysis to perform. Usually, you run Stat-Analysis to aggregate data across multiple cases through time. Presumably though your use case would just generate a handful of MPR matched pair lines. I suppose you could have Stat-Analysis convert that to SL1L2 partial sums and/or apply threshold to generated CTC contingency table count output lines.

Just want to make sure that Stat-Analysis is actually performing a useful step.

[malloryprow]

@JohnHalleyGotway Thanks! StatAnalysis also makes sense. We just need the SL1L2 partial sums, so that would be great.

[JohnHalleyGotway]

@malloryprow, OK, so the work here is developing a python embedding script that serves up the paired data in the MPR line type format. Here's a very simple example of running a python script that literally just reads a "*_mpr.txt" file and passes the result to Stat-Analysis to compute partial sums. I ran this from the top-level MET installation directory, so you'd need to adjust paths accordingly.

Of course, if you already have "*_mpr.txt" files, you'd typically just pass that as input to Stat-Analysis directly. This is just a simple example to demonstrate how the python embedding of MPR data works.

bin/stat_analysis \
-lookin python share/met/python/read_ascii_mpr.py out/point_stat/*_mpr.txt \
-job aggregate_stat -line_type MPR -out_line_type SL1L2 \
-by FCST_VAR,FCST_LEV -out_stat agg_sl1l2.stat

The real work here is copying/modifying the read_ascii_mpr.py script to compute the fcst/obs sea ice area, define all the metadata (i.e. time stamps, message type, variable name/level, etc) for the 21 .stat header columns, and then formatting it.

Let me know if you get stuck along the way. As you work through the details, please let us know of any suggestions you have to make the process easier in the future.

For reference, here's the output I see when running the command above:

musial6:MET-develop johnhg$ bin/stat_analysis -lookin python share/met/python/read_ascii_mpr.py out/point_stat/*_mpr.txt -job aggregate_stat -line_type MPR -out_line_type SL1L2 -by FCST_VAR,FCST_LEV -out_stat agg_sl1l2.stat
DEBUG 1: Start stat_analysis by johnhg(6088) at 2022-12-20 17:17:28Z  cmd: bin/stat_analysis -lookin python share/met/python/read_ascii_mpr.py out/point_stat/point_stat_360000L_20070331_120000V_mpr.txt -job aggregate_stat -line_type MPR -out_line_type SL1L2 -by FCST_VAR,FCST_LEV -out_stat agg_sl1l2.stat 
DEBUG 1: Default Config File: /Volumes/d1/projects/MET/MET_development/MET-develop/share/met/config/STATAnalysisConfig_default
Python Script:	'share/met/python/read_ascii_mpr.py'
Input File:	'out/point_stat/point_stat_360000L_20070331_120000V_mpr.txt'
Data Length:	27984
Data Type:	<class 'list'>
DEBUG 2: 
DEBUG 2: Processing Job 1: -job aggregate_stat -line_type MPR -by FCST_VAR -by FCST_LEV -out_stat agg_sl1l2.stat -out_line_type SL1L2 
DEBUG 1: Creating STAT output file "agg_sl1l2.stat"
GSL_RNG_TYPE=mt19937
GSL_RNG_SEED=820945251
JOB_LIST:      -job aggregate_stat -line_type MPR -by FCST_VAR -by FCST_LEV -out_stat agg_sl1l2.stat -out_line_type SL1L2 
DEBUG 2: Computing output for 3 case(s).
DEBUG 2: For case "TMP:P900-750", found 2 unique OBTYPE values: ADPUPA,ADPSFC
DEBUG 2: For case "TMP:P900-750", found 3 unique VX_MASK values: DTC165,DTC166,LMV
WARNING: For case "TMP:P900-750", found 3 unique INTERP_MTHD values: NEAREST,MEDIAN_SQUARE,DW_MEAN_SQUARE
WARNING: For case "TMP:P900-750", found 2 unique INTERP_PNTS values: 1.00000,9.00000.
DEBUG 2: For case "UGRD:Z10", found 3 unique VX_MASK values: DTC165,DTC166,LMV
WARNING: For case "UGRD:Z10", found 3 unique INTERP_MTHD values: NEAREST,MEDIAN_SQUARE,DW_MEAN_SQUARE
WARNING: For case "UGRD:Z10", found 2 unique INTERP_PNTS values: 1.00000,9.00000.
DEBUG 2: For case "VGRD:Z10", found 3 unique VX_MASK values: DTC165,DTC166,LMV
WARNING: For case "VGRD:Z10", found 3 unique INTERP_MTHD values: NEAREST,MEDIAN_SQUARE,DW_MEAN_SQUARE
WARNING: For case "VGRD:Z10", found 2 unique INTERP_PNTS values: 1.00000,9.00000.
DEBUG 2: Job 1 used 27984 out of 27984 STAT lines.
DEBUG 1: Finish stat_analysis by johnhg(6088) at 2022-12-20 17:17:45Z

[malloryprow]

Thanks. I have been playing around with a script this morning, and am making some progress. I'll keep you all updated!

[JohnHalleyGotway]

@malloryprow wondering if we can mark this discussion as being answered and lock it to prevent future posts? Of course you can always create a new discussion if new issues or questions arise.

[malloryprow]

Yes, feel free to mark as answered!