gen vx mask

[fearon10]

Hi There,

This question is similar to #1463 , but a little different. I am also using python embedding. I would like to create standalone masks for land, ocean, and sea ice based on the following mask.JPG, with land=1, ocean=0, and sea ice=2.

I able to create a land.nc mask using the following:

gen_vx_mask PYTHON_NUMPY PYTHON_NUMPY land.nc -type data -input_field 'name="readdata.py -f data.hdf5 -l 10 -lt model_level -v sea_land_ice_mask";' -mask_field 'name="readdata.py -f data.hdf5 -l 10 -lt model_level -v sea_land_ice_mask";' -thresh eq1 intersection

I am also able to create a sea ice mask with the same logic, except thresh eq 2 is used:

gen_vx_mask PYTHON_NUMPY PYTHON_NUMPY ice.nc -type data -input_field 'name="readdata.py -f data.hdf5 -l 10 -lt model_level -v sea_land_ice_mask";' -mask_field 'name="readdata.py -f data.hdf5 -l 10 -lt model_level -v sea_land_ice_mask";' -thresh eq2 -intersection

However, I am not able to create an ocean mask using thresh eq0:

gen_vx_mask PYTHON_NUMPY PYTHON_NUMPY sea.nc -type data -input_field 'name="readdata.py -f data.hdf5 -l 10 -lt model_level -v sea_land_ice_mask";' -mask_field 'name="readdata.py -f data.hdf5 -l 10 -lt model_level -v sea_land_ice_mask";' -thresh eq0 -intersection

I've tried other combinations too. For example, using my new land.nc as input and the mask field as follows:

gen_vx_mask land.nc land.nc sea.nc -type data -mask_field 'name="data_mask"; level="(,)";' -thresh eq1 -complement

This and the previous command gives a raster with all grid cells equal to 1 as follows:

I think I'm just misunderstanding some aspect of the gen_vx_mask command. Any assistance would be much appreciated.

Thanks,
Matt

[JohnHalleyGotway]

Hello Matt,

I see that you'd like to run Gen-Vx-Mask 3 times to create 3 different stand-alone masks. I think that the commands you're running are slightly more complicated than they need to be.

Gen-Vx-Mask has 4 required arguments:

Usage: gen_vx_mask
	input_grid
	mask_file
	out_file
	-type string

The input_grid must define the grid for which the mask should be defined. That can be a named grid, the path to a gridded data file, or an explicit grid specification string. You're currently using python embedded to serve up gridded data. Based on the image you sent, this looks like a global lat/lon grid. And that might already be defined in MET as one of NCEP's predefined grids. Check this list of grid definitions. For example, if this is a global 0.5 degree lat/lon grid, that'd correspond to G004 in MET. And that'd simplify your commands like this:

# For Land
gen_vx_mask G004 PYTHON_NUMPY land.nc \
-type data \
-mask_field 'name="readdata.py -f data.hdf5 -l 10 -lt model_level -v sea_land_ice_mask";' \
-thresh eq1

# For Ice
gen_vx_mask G004 PYTHON_NUMPY ice.nc \
-type data \
-mask_field 'name="readdata.py -f data.hdf5 -l 10 -lt model_level -v sea_land_ice_mask";' \
-thresh eq2

# For Water
gen_vx_mask G004 PYTHON_NUMPY water.nc \
-type data \
-mask_field 'name="readdata.py -f data.hdf5 -l 10 -lt model_level -v sea_land_ice_mask";' \
-thresh eq0

If G004 is NOT correct and you can't find a pre-defined one in that list, another option is just explicitly defining a grid specification string.

Or if you would like to purse the "double python embedding" approach you're using for other reasons, I'm happy to debug further. I'd just need the readdata.py and data.hdf5 files to replicate the behavior. You can post sample data to our anonymous ftp site following these instructions. Please just let me know if you decide to post data so I know to go grab it.

Thanks,
John

[fearon10]

Thanks John.
I took your advice and simplified my gen_vx_mask steps. To accomplish this, I converted my hdf5 file to netcdf outside of MET (filename=test.nc). I figured this would allow me to set aside the python embedding for now and perform a direct read of the file into MET. I also found the appropriate grid specification which is G193. With these changes, I found that the gen_vx_mask steps you suggested work, with the caveat that my netcdf lat/long order is exactly the same as the G193 grid.

for water:
gen_vx_mask G193 test.nc water.nc -type data -mask_field 'name="sea_land_ice_mask"; level="(0,,)"; file_type=NETCDF_NCCF;' -thresh eq0
for land:
gen_vx_mask G193 test.nc land.nc -type data -mask_field 'name="sea_land_ice_mask"; level="(0,,)"; file_type=NETCDF_NCCF;' -thresh eq1
for ice:
for water:
gen_vx_mask G193 test.nc ice.nc -type data -mask_field 'name="sea_land_ice_mask"; level="(0,,)"; file_type=NETCDF_NCCF;' -thresh eq2

While the former steps work, the following do not, where the G193 grid is replaced with my test.nc file

for water:
gen_vx_mask test.nc test.nc land.nc -type data -mask_field 'name="sea_land_ice_mask"; level="(0,,)"; file_type=NETCDF_NCCF;' -thresh eq0

Nor,
gen_vx_mask test.nc test.nc land.nc -type data -input_field 'name="sea_land_ice_mask"; level="(0,,)"; file_type=NETCDF_NCCF;' -mask_field 'name="sea_land_ice_mask"; level="(0,,)"; file_type=NETCDF_NCCF;' -thresh eq0

Is there something very specific that needs to be in the input_grid metadata that gen_vx_mask uses? Strange how using test.nc as the input_grid does not work. Here is a header dump of test.nc.

netcdf test {
dimensions:
latitude = 721 ;
longitude = 1440 ;
level = 15 ;
time = UNLIMITED ; // (1 currently)
forecast_reference_time = 1 ;
variables:
double time(time) ;
time:standard_name = "forecast_valid_time" ;
time:long_name = "forecast_valid_time" ;
time:units = "hours since 1900-01-01 00:00:00.0" ;
time:calendar = "gregorian" ;
time:axis = "T" ;
double forecast_reference_time(forecast_reference_time) ;
forecast_reference_time:standard_name = "forecast_initialization_time" ;
forecast_reference_time:long_name = "forecast_initialization_time" ;
forecast_reference_time:units = "hours since 1900-01-01 00:00:00.0" ;
forecast_reference_time:calendar = "gregorian" ;
float latitude(latitude, longitude) ;
latitude:standard_name = "latitude" ;
latitude:long_name = "latitude" ;
latitude:units = "degrees_north" ;
latitude:axis = "Y" ;
float longitude(latitude, longitude) ;
longitude:standard_name = "longitude" ;
longitude:long_name = "longitude" ;
longitude:units = "degrees_east" ;
longitude:axis = "X" ;
float sea_land_ice_mask(time, latitude, longitude) ;
sea_land_ice_mask:_FillValue = -32767.f ;
sea_land_ice_mask:standard_name = "sea_land_ice_mask" ;
sea_land_ice_mask:long_name = "sea_land_ice_mask" ;
sea_land_ice_mask:units = "unitless" ;
sea_land_ice_mask:missing_value = -32767.f ;
sea_land_ice_mask:coordinates = "latitude longitude" ;

// global attributes:
:source = "test NETCDF4" ;
:history = "Created Wed Apr 19 01:00:22 2023" ;
:description = "Conventions = CF-1.6" ;
}

Thanks for the help,
Matt

[JohnHalleyGotway]

Matt, please try specifying the global conventions attribute in a CF-compliant way.
Replace :description = "Conventions = CF-1.6" ;
With :Conventions = "CF-1.6" ;

You can read about that in this section. Does that produce a better result?

[fearon10]

Hi John,
Specifying the global attribute in a CF-complaint way, i.e., ":Conventions = "CF-1.6" ;" worked. I can generate an independent mask for land, ocean, and ice now with netcdf.
In addition, I resolved my python embedding issue. The input_grid and mask_grid are the same, but input_field and mask_field must be different. For example, the following now works for water.

gen_vx_mask PYTHON_NUMPY PYTHON_NUMPY ocean.nc -type data -input_field 'name="readdata.py -f data.hdf5 -l 10 -lt model_level -v zonal_wind_at_10m";' -mask_field 'name="readdata.py -f data.hdf5 -l 10 -lt model_level -v sea_land_ice_mask";' -thresh eq0

Previously, my logic with embedding was the same except my input_field and mask_field were the same, i.e., sea_land_ice_mask for both, and that must cause an issue. Specifying another variable in the hdf5 file (with the same dimensions as the mask variable) for input_field seems to do the trick.

Thank you,
Matt

[JohnHalleyGotway]

Matt,

Excellent! I'm glad you got it all sorted out.

Let me clarify one detail about how -input_field and -mask_field options work. Hopefully this will shed some light on the confusing behavior you witnessed. When doing -type data masking, the -mask_field is required. It defines the gridded data values to be thresholded to define the current masking operation.

However, the -input_field is always optional. It is generally intended as a way of combining masks across multiple runs to define complex regions. By default, with data masking (and all Gen-Vx-Mask masking operations), the initial value at each grid point is set to 0, as seen on this line of code. When you specify the -input_field command line option, it initializes the value at each grid point to whatever field you read in rather than 0. Any -input_field value = 0 is interpreted as the input mask being "off" and any value != 0 is interpreted as the mask being "on".

So the logic proceeds like this:

1. We have the field of initial values - either all 0's or whatever was specified with -input_field. Remember, =0 means mask "off" and !=0 means mask "on".
2. We apply the user-specified threshold to the -mask_field values to create a 0/1 binary mask.
3. We apply the user-configurable logic (-union, -intersection, -symdiff) to combine the input value with the mask value to determine the output value. That combination logic can be seen starting on this line of code.
4. The resulting values from the combination logic are written to the output file.

You note that when you used the same field for -input_field and -mask_field with -thresh eq0, you got a field of all 1's in the output. Looking at the default behavior from the combination logic linked above, it makes some sense:

            // Default behavior is to apply the mask value or pass through
            // the data value
            default:
               if(v_mask) v = mask_val;
               else       v = dp_data(x, y);
               break;

For each grid point...

• if it equals 0 (i.e. is water), then v_mask is true and the output value is set to mask_val (which is 1 by default).
• if it does not equal 0 (i.e. is not water), then we just pass whatever the -input_field value is to the output.

If you want to play around with this more, you could use -value 100 (or something) to write a mask value of 100 instead of the default mask value of 1. That'd enable you to differentiate in the output what values are being "passed through" versus what are being "masked out" in this step.

But here's the real punchline, this is all WAY TOO CONFUSING for the simple operation you'd like to perform. And that's why I recommended avoiding using -input_field at all because it introduces complexity to the logic. It is very powerful and flexible when you run Gen-Vx-Mask multiple times, but it's definitely overkill for what you're doing.

If you find it convenient to continue using -input_field, you can simplify by resetting all the input data values to 0. That'll revert back to the default logic and the tool will behave more intuitively. Use the censor_thresh and censor_val options like this to do so:

-mask_field 'name="readdata.py -f data.hdf5 -l 10 -lt model_level -v sea_land_ice_mask"; censor_thresh=ne-9999; censor_val=0.0;'

That'll replace any non-bad data value (i.e. -9999) with 0.

Hopefully that helps clarify.

[fearon10]

Thanks John.
This additional information you have provided is very helpful. One follow-up question I have is: if my desired grid is in a format not directly readable by MET (e.g., hdf5), isn't the input_field required? In my case, the input_grid will be PYTHON_NUMPY.

[JohnHalleyGotway]

@fearon10 just want to clarify a bit here. Gen-Vx-Mask has 3 required arguments, and we are talking about the first one of those, named input_grid:

Usage: gen_vx_mask
	input_grid
	mask_file
	out_file

You have 3 choices for how to specify that input_grid:

1. The name of a pre-defined grid, such as G004 for NCEP Grid number 4.
2. An explicit grid specification string, as described in Appendix B.
3. As the path a gridded data file.

And you're talking about option (3)... but using python embedding to serve up that data rather than a gridded data file that MET can read natively.

The answer to your question is yes. If you're using (3) to supply gridded data via python to define the input grid, then you'll need to supply the -input_field command line option to provide the python script to be run.

It's up to you to decide if option (1), (2), or (3) is easiest for you. If you stick with (3), I'd recommend serving up data values of all 0's so as to avoid the potentially confusing logic of combining the input field value with the current mask value to define the output value.