Mask Pressure Levels Below Ground Level with Grid-Stat

[ShawnMurdzek-NOAA]

I'm trying to verify model output using Grid-Stat with my "truth" dataset being output from another model. Output from both models include standard fields (T, humidity, winds) on several pressure levels. Some of these pressure levels, however, lie beneath the surface (e.g., 925 mb output in the intermountain west). I would like to mask those grid points that lie beneath the surface so they are not included in the verification. For what it's worth, both models produce non-masked values of T, humidity, winds, etc. for all pressure levels, even if they lie beneath the surface (i.e., the model output does not mask any of the values below the surface).

My first thought was to try to create a mask similar to the output from gen_vx_mask. I would need such a mask to have dimensions of (lat, lon, pressure) because there are several pressure levels I am performing verification for. I have thus far been unable to create a 3D mask using gen_vx_mask, however. I tried the following, but still ended up with a 2D mask with the dimensions of (lat, lon):

gen_vx_mask \
  ${in_file} \
  ${in_file} \
  ${out_file} \
  -type data \
  -mask_field 'name="HGT"; level="P900, P875";' \
  -thresh 'gt1000' \

Are 3D masks not supported by MET? If not, is there a better way for me to go about masking pressure levels that lie beneath the surface?

Thanks in advance for any guidance and suggestions!

[JohnHalleyGotway]

Hello @ShawnMurdzek-NOAA I see you have a question about excluding data from your evaluation that falls below the model surface.

First, you're correct to note that the gen_vx_mask tool always generated simple 2D mask fields with lat/lon dimensions. It does not generate 3-dimensional data... for varying output times or pressure levels.

Next, this does not answer your question, but is somewhat related. The topo_mask configuration option provides related functionality, but for the Point-Stat tool. When verifying at the surface, it compares the model elevation to the height of the observing location. And thresholding that difference (based on config settings), allows for the point observations to be excluded from the statistics and/or determine which forecast grid points should be used when interpolating to that point observation location. Unfortunately, no analogous configuration option is provided for Grid-Stat. So no easy/obvious fix for this.

I tested with some sample GFS data which includes the heights (HGT) of the surface and many pressure levels. For example:

164:30852225:d=2012040900:HGT:975 mb:12 hour fcst:
178:33861945:d=2012040900:HGT:1000 mb:12 hour fcst:
180:34581404:d=2012040900:HGT:surface:12 hour fcst:

It seems to me that you'd only want to verify at grid points where surface height is >= the height of that pressure level. Is that correct?

One option is writing a Python script to apply this logic, setting the below-ground points to MET's internal bad data value of -9999. And then you could use Python embedding in MET to serve up this data to Grid-Stat. But if you're not already familiar and comfortable with Python embedding, then this would take time.

Another option is running PCP-Combine multiple times to compute the difference described above. For example:

pcp_combine -subtract gfs_2012040900_F012.grib2 'name="HGT"; level="P1000";' gfs_2012040900_F012.grib2 'name="HGT"; level="L0";' -name HGT_P1000_minus_surface HGT_P1000_minus_surface.nc
pcp_combine -subtract gfs_2012040900_F012.grib2 'name="HGT"; level="P850";' gfs_2012040900_F012.grib2 'name="HGT"; level="L0";' -name HGT_P850_minus_surface HGT_P850_minus_surface.nc

And then in Grid-Stat, you could use "data masking" to read data and threshold it:

obs = {

  field = [
    {
      name = "TMP";
      level = "P1000";
      mask = {
        poly = [ "HGT_P1000_minus_surface.nc { name=\"HGT_P1000_minus_surface\"; level=\"(*,*)\"; } ge0.0" ];
      }
    },
    {
      name = "TMP";
      level = "P850";
      mask = {
        poly = [ "HGT_P850_minus_surface.nc { name=\"HGT_P850_minus_surface\"; level=\"(*,*)\"; } ge0.0" ];
      }
    }
   ];
}

So this demonstrates that you can define masks separately for each verification task and the mask.poly option enables you to specify an input gridded data file, followed by the data to be read from that file, and followed by a threshold to apply to define the region of interest.

I do think it sure would be nice to enhance the MET statistics tools to address this situation more directly without needing to run pcp_combine and/or gen_vx_mask multiple times upstream.

[ShawnMurdzek-NOAA]

Hi @JohnHalleyGotway,

Thanks for the response! To answer your question, yes, I only want to verify grid points if the height at that pressure level >= the surface height.

Both avenues you suggested seem like good paths forward. I am going to try to attempt the first path, but without the Python embedding (i.e., modify the model fields before running MET). I'm running verification on several model runs with the same "truth" run being used for verification in all cases. Based on your answer, it seems like I should be able to simply change any pressure levels beneath the surface in that "truth" model run to be -9999, and then MET will ignore those grid points during verification.

Can I keep this question open until I confirm that this path works? Hopefully it will only take a day or two to test.