Looking for Examples of Python Embedding with Xarray

[johnlwagner]

Greetings
We are working with python embedding for the first time. Does anyone have a working example of a script using python embedding with xarray that they could share with us? We have been following the examples in Appendix F, but have not gotten a successful test yet. We have looked at some other discussions but haven't found anything that we got to work for us.
We are working on WCOSS and are extracting a single QPF06 grid from a netCDF file that has multiple elements and lead times.

In addition to having issues with python embedding, we may also have an issue with a MET variable. When trying to check a netcdf file with plot_data_plane, we get the following error:
can't open file '/apps/ops/para/libs/intel/19.1.3.304/met/11.1.0-rcl/share/met/wrappers/write_tmp_dataplane.py': [Errno 2] No such file or directory
That file does exist in /apps/ops/para/libs/intel/19.1.3.304/met/11.1.0-rc1/share/met/python/pyembed and in share/met/wrappers for other builds. Testing other versions of MET did not fix this for us. Is there a variable that can be set to change the location of the wrappers directory?

[DanielAdriaansen]

Hi @johnlwagner, thanks for your question. I developed this tutorial for Python Embedding which included writing a Python embedding script to use with plot_data_plane using both PYTHON_NUMPY and PYTHON_XARRAY to plot GOES-16 data:
https://dtcenter.org/metplus-practical-session-guide-version-5-0/session-9-python-embedding/writing-python-script-python-embedding

Can you look at the script there and compare with yours and see if anything stands out? Also please feel free to share your script and command with us, and we will be able to better assist you, and also log files. Without seeing those it's hard to know exactly what might be causing you trouble when trying to use Python embedding with Xarray.

Before I reply about the second part of your question, can you please share your plot_data_plane command and any Python embedding scripts that are being used? There shouldn't be a need to force MET to look for a particular file when using Python embedding.

[johnlwagner]

Thanks Dan. Your link did help and I reworked my python script to match the example more. I believe I have the attributes correct, but there is still an issue with the grid_info. I have a LatLon grid and here is what I have for grid_info:

grid_info= {'type': 'LatLon', 'name': 'accesse', 'lat_ll': -89850000.0, 'lon_ll': 225000.0, 'delta_lat': 450000.0, 'delta_lon': 300000.0, 'Nlat': 800, 'Nlon': 600}

The error I get is

Traceback (most recent call last):
  File "/apps/ops/para/libs/intel/19.1.3.304/met/11.0.2/share/met/wrappers/write_tmp_dataplane.py", line 55, in <module>
    nx, ny = met_in.met_data.shape
ValueError: too many values to unpack (expected 2)
ERROR  : 
ERROR  : tmp_nc_dataplane() -> command "/apps/spack/python/3.8.6/intel/19.1.3.304/pjn2nzkjvqgmjw4hmyz43v5x4jbxjzpk/bin/python3 /apps/ops/para/libs/intel/19.1.3.304/met/11.0.2/share/met/wrappers/write_tmp_dataplane.py /tmp/tmp_met_nc_2859290_0 $SCRverif/practice_gridded_pyembed $ncFile $testDir $proj" failed ... status = 256
ERROR  : 

I have Nlat and Nlon set, so I'm not sure how nx,ny have too many values. My scripts are below. I wasn't able to attach them.
If there is a better way to send you scripts, please let me know. test_regrid.sh calls practice_gridded_pyembed.py. Please let me know if you see anything that I'm doing wrong.

test_regrid.sh:

#!/bin/sh
#
export proj=12
export init_date=2025011500
export valid_date=$($NDATE +$proj $init_date)
export initYYYYMMDD=`echo $init_date | cut -c1-8`
export initHH=`echo $init_date | cut -c9-10 | sed "s/^0//g"`
export validYYYYMMDD=`echo $valid_date | cut -c1-8`
export validHH=`echo $valid_date | cut -c9-10 | sed "s/^0//g"`
#
. /lfs/h3/mdl/verif/noscrub/usr/john.l.wagner/mdl-verification/verif_v2.1.0/fix/defined_paths_h3.env
export ncFile=/lfs/h1/ops/prod/com/blend/v4.2/blend.${initYYYYMMDD}/$(printf %02d $initHH)/modeldata/blend.t$(printf %02d $initHH)z.accesse_01_apcp.nc
if [[ ! -s $ncFile ]]; then echo "ncFile $ncFile not found"; exit; fi
export testDir=/lfs/h3/mdl/ptmp/john.l.wagner/component/accesse_01/co/dbnc/${validYYYYMMDD}/${valid_date}/q6_component
mkdir -p $testDir
export MET_PYTHON_LD="-L/apps/spack/python/3.8.6/intel/19.1.3.304/pjn2nzkjvqgmjw4hmyz43v5x4jbxjzpk/lib  -lcrypt -lintl -ldl  -lutil -lm -lm"
export MET_PYTHON_CC="-I/apps/spack/python/3.8.6/intel/19.1.3.304/pjn2nzkjvqgmjw4hmyz43v5x4jbxjzpk/include/python3.8 -I/apps/spack/python/3.8.6/intel/19.1.3.304/pjn2nzkjvqgmjw4hmyz43v5x4jbxjzpk/include/python3.8  -fno-strict-aliasing -Wsign-compare -Wunreachable-code -fPIC -DNDEBUG -g  -O3 -Wall"
export MET_PYTHON_EXE=/apps/spack/python/3.8.6/intel/19.1.3.304/pjn2nzkjvqgmjw4hmyz43v5x4jbxjzpk/bin/python3
echo $SCRverif
echo "$SCRverif/test_nc.py $ncFile $testDir $proj"
#
plot_data_plane PYTHON_XARRAY $ncFile \
  'name="$SCRverif/practice_gridded_pyembed.py $ncFile $testDir $proj";' \
  -title "Python enabled plot_data_plane"

practice_gridded_pyembed.py:

import xarray as xr
import os
import sys
import cftime
#
# Test script for python embedding in MET
#
# Get the command line arguments
input_file = sys.argv[1]
outdir = sys.argv[2]
proj = sys.argv[3]

# Open the input file as an xarray dataset
qpf06 = xr.open_dataset(input_file,decode_times=False)

# Initialize the time variables based on the names
# of the input file and output directory
init_date=input_file.split("/")[8].split(".")[1]
init_cycle=input_file.split("/")[9]
valid_date=outdir.split("/")[10]
valid_cycle=outdir.split("/")[11][-2:]
init=init_date+"_"+init_cycle+"0000"
valid=valid_date+"_"+valid_cycle+"0000"
lead=proj+"0000"
name=outdir.split("/")[7]


# Loop over the data_vars in the input file to find the 
# variable and lead time that we need
for var_name in qpf06.data_vars:
	var_data = qpf06[var_name]
	attrs2=var_data.attrs
	if var_name == 'apcp_06h':
		var_data.attrs=qpf06.attrs
		qpf06_data=var_data.sel(lead_time_06h=int(proj)) 
		#print(attrs2)
		#print(qpf06_data.attrs)
		#print(qpf06_data.attrs['model_family'])
		#print("coords=",qpf06_data.coords)
		#print("dims=",qpf06_data.dims)
		#print("values=",qpf06_data.values)
# Construct a python dictionary names attrs that contains
# the required attributes needed for python embedding
		attrs = {}
		attrs['valid'] = valid_date 
		attrs['init'] = init_date
		attrs['lead'] = lead
		attrs['accum'] = '000000'
		attrs['name'] = qpf06_data.attrs['model_family']
		attrs['long_name'] = attrs2['long_name']
		attrs['level'] = attrs2['level']
		attrs['units'] = attrs2['units']
# Set met_data and its attributes
		met_data = qpf06_data
		met_data.attrs = attrs
	elif var_name == 'grib2_gds':
		#print(var_name)
		#print(var_name,var_data.attrs)
		#print("grib2 coords=",var_data.coords)
		#print("grib2 values=",var_data.values)
		#print("grib2 dims=",var_data.dims)
# Construct a python dictionary names grid_info that contains
# the required grid information needed for python embedding
		grid_info = {}
		grid_info['type'] = 'LatLon'
		grid_info['name'] = attrs['name']
		grid_info['lat_ll'] = float(var_data.values[11])
		grid_info['lon_ll'] = float(var_data.values[12])
		grid_info['delta_lat'] = float(var_data.values[16])
		grid_info['delta_lon'] = float(var_data.values[17])
		grid_info['Nlat'] = var_data.values[7]
		grid_info['Nlon'] = var_data.values[8]
		print("grid_info=",grid_info)
		attrs['grid'] = grid_info

[DanielAdriaansen]

Hi @johnlwagner glad to hear the example helped! I welcome any feedback on how that could be made better if you found any of it confusing.

Back to your question, I have some questions about this part:

		grid_info['lat_ll'] = float(var_data.values[11])
		grid_info['lon_ll'] = float(var_data.values[12])
		grid_info['delta_lat'] = float(var_data.values[16])
		grid_info['delta_lon'] = float(var_data.values[17])
		grid_info['Nlat'] = var_data.values[7]
		grid_info['Nlon'] = var_data.values[8]

Earlier in your script, it appears var_data for the grid_info is set like this:

	var_data = qpf06['grib2_gds']

My guess here, without having your data file and being able to inspect what var_data looks like for the variable grib2_gds, is that this:

var_data.values[7]

is returning an object that is not just a single integer value like 800. Instead, it's returning a tuple or some other type of Python object that is not a single value. I would recommend inspecting var_data['grib2_gds'] and all of the elements of it you're using to make sure you are providing a single integer value for grid_info['Nlat'] and grid_info['Nlon'].

Let me know what you find out those items look like, which may help.

[johnlwagner]

It looks like the lat/lon values are indeed integers and not tuples. I set them both to ints just to be sure. grid_info['Nlat'] = int(var_data.values[7]) grid_info['Nlon'] = int(var_data.values[8]) I'm still getting the same error: Traceback (most recent call last): File "/apps/ops/para/libs/intel/19.1.3.304/met/11.0.2/share/met/wrappers/write_tmp_dataplane.py", line 55, in <module> nx, ny = met_in.met_data.shape ValueError: too many values to unpack (expected 2) ERROR : ERROR : tmp_nc_dataplane() -> command "/apps/spack/python/3.8.6/intel/19.1.3.304/pjn2nzkjvqgmjw4hmyz43v5x4jbxjzpk/bin/python3 /apps/ops/para/libs/intel/19.1.3.304/met/11.0.2/share/met/wrappers/write_tmp_dataplane.py /tmp/tmp_met_nc_529407_0 $SCRverif/practice_gridded_pyembed $ncFile $testDir $proj" failed ... status = 256 ERROR : Here's what all of my grid_info entries look like: grid_info= {'type': 'LatLon', 'name': 'accesse', 'lat_ll': -89850000.0, 'lon_ll': 225000.0, 'delta_lat': 450000.0, 'delta_lon': 300000.0, 'Nlat': 800, 'Nlon': 600}

…

On Wed, Jan 22, 2025 at 4:19 PM Dan Adriaansen ***@***.***> wrote: Hi @johnlwagner <https://github.com/johnlwagner> glad to hear the example helped! I welcome any feedback on how that could be made better if you found any of it confusing. Back to your question, I have some questions about this part: grid_info['lat_ll'] = float(var_data.values[11]) grid_info['lon_ll'] = float(var_data.values[12]) grid_info['delta_lat'] = float(var_data.values[16]) grid_info['delta_lon'] = float(var_data.values[17]) grid_info['Nlat'] = var_data.values[7] grid_info['Nlon'] = var_data.values[8] Earlier in your script, it appears var_data for the grid_info is set like this: var_data = qpf06['grib2_gds'] My guess here, without having your data file and being able to inspect what var_data looks like for the variable grib2_gds, is that this: var_data.values[7] is returning an object that is not just a single integer value like 800. Instead, it's returning a tuple or some other type of Python object that is not a single value. I would recommend inspecting var_data['grib2_gds'] and all of the elements of it you're using to make sure you are providing a single integer value for grid_info['Nlat'] and grid_info['Nlon']. Let me know what you find out those items look like, which may help. — Reply to this email directly, view it on GitHub <#2864 (reply in thread)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AGZ7TYIOXF6W4XVW7JBMO4D2MADM3AVCNFSM6AAAAABVH6Q2YKVHI2DSMVQWIX3LMV43URDJONRXK43TNFXW4Q3PNVWWK3TUHMYTCOJSGE4TCOI> . You are receiving this because you were mentioned.Message ID: ***@***.***>

-- John Wagner Verification Task Lead NOAA/National Weather Service Meteorological Development Laboratory Digital Forecast Services Division SSMC2 Room 10106 Silver Spring, MD 20910 (240) 622-9138 <https://mail.google.com/mail/u/0/#m_226873644979562856_> (Google Voice) (301) 427-9471 (office) (908) 902-4155 (cell/text)

[DanielAdriaansen]

OK thanks for checking on that. Another idea this morning after reading your error again, is that the actual data you are sending to MET is not 2D.

Can you verify that this line:

		met_data = qpf06_data

that the variable met_data actually has two dimensions and not three or more? For example, in the past I've encountered NetCDF files with a time dimension, such that the variable qpf06 has dimensions of (time,nlat,nlon). If that's the case, it could be that met_data has the shape of (1,800,600), which MET does not support. You'll need to reduce the shape of met_data to 2D prior to calling MET.

Let me know what you find. If this turns out to be a red herring, could you send your data file so I could try and run it myself? I am happy to help by reading error messages and guessing but having the data file would probably allow it to go faster. You can read how to share data here:
#954

[johnlwagner]

Yes, you were right, I did have a 3d array. Now that I have a 2d array, I am seeing different errors. I'm seeing a conflict with python v3.8 and numpy. When testing with python v3.10, there seem to be other python conflicts. I will see if I can fix the python issues on my end and follow up with you.