Point_stat: Python Embedding for Forecast Fields, attrs dictionary question

[justintsu]

Hello,

I am trying to pass in a model field into point_stat using python embedding. I have a python routine that is communicating two main variables to point_stat: met_data (numpy array) and attrs (standard dictrionary).

• met_data is obtained by reading in a model binary field and reformatting it into an m x n matrix.
• attrs is obtained by:

1. gleaning information from the filename such as valid time, lead time, level, etc.
2. reading in a different binary meta data file that contains aspects about the grid such as projection. The data from this file will populate the grid 'nested' dictionary that will be placed in attrs. Here is how I am defining the grid, based on the projection (nproj) that is defined in the metadata file:

if nproj == 1:
    lat_str = 'latitu_sfc_000000_000000_{1}_{0}_00000000_fcstfld'.format(dtg_str,tokens[4])
    lon_str = 'longit_sfc_000000_000000_{1}_{0}_00000000_fcstfld'.format(dtg_str,tokens[4])
    lat_f = '/'.join(ff_str.split('/')[:-1]+[lat_str])
    lon_f = '/'.join(ff_str.split('/')[:-1]+[lon_str])
    lats = np.float64(np.fromfile(lat_f, '>f'))
    lons = np.float64(np.fromfile(lon_f, '>f'))
    grid['type']   = "Mercator"
    grid['lat_ll'] = np.min(lats)
    grid['lon_ll'] = np.min(lons)
    grid['lat_ur'] = np.max(lats)
    grid['lon_ur'] = np.max(lons)
elif nproj == 2:
    grid['type']        = "Lambert Conformal"
    grid['scale_lat_1'] = float(datahd[3])
    grid['scale_lat_2'] = float(datahd[4])
    grid['lon_orient']  = float(datahd[5])
    grid['lat_pin']     = float(datahd[6])
    grid['lon_pin']     = float(datahd[7])
    grid['x_pin']       = float(datahd[ind+3])
    grid['y_pin']       = float(datahd[ind+4])

    #checking for negative lon values, put on 0-360 scale
    if grid['lon_orient'] < 0. : grid['lon_orient'] = grid['lon_orient'] + 360.
    if grid['lon_pin'] < 0.: grid['lon_pin'] = grid['lon_pin'] + 360.
elif nproj == 3:
    grid['type']       = "Polar Stereographic"
    grid['hemisphere'] = "N"
    grid['scale_lat']  = float(datahd[3])
    grid['lat_pin']    = float(datahd[6])
    grid['lon_pin']    = float(datahd[7])
    grid['x_pin']      = float(datahd[ind+3])
    grid['y_pin']      = float(datahd[ind+4])
    grid['lon_orient'] = float(datahd[5])
    grid['d_km']       = float(datahd[ind+6])/1000.
    grid['r_km']       = 6361.2
    if grid['lon_orient'] < 0. : grid['lon_orient'] = grid['lon_orient'] + 360.
    if grid['lon_pin'] < 0.: grid['lon_pin'] = grid['lon_pin'] + 360.

How does point_stat use the grid sub dictionary to determine its forecast grid?
I have separate latitude and longitude files that I can read in that explicitly define the latitude and longitude grids in the same m x n dimensions as met_data. Is it possible to just explicitly define the grid and pass this into point_stat without having to do what I am doing above?

[jprestop]

Hi @justintsu. I see that this question was submitted on July 26, but has gone unanswered. I apologize that this question was missed. Unfortunately, the person best able to answer your question is out until later in the week. @JohnHalleyGotway could you please take a look at this question when you get back into the office?

[justintsu]

Thank you for tagging up. I am looking forward to a response when John gets back.

[DanielAdriaansen]

Hi @justintsu, thanks for reaching out, and for your patience.

If the separate latitude and longitude files you have defining the grid are in a supported format that MET can read, then I think you can provide the path to this file as the grid entry in the attrs dictionary. Here is an example test I just did:

attrs['grid'] = '/d1/projects/MET/MET_test_data/unit_test/model_data/met_nc/nam/nam_2012040900_F036_APCP12.nc'

And MET reports this:

DEBUG 3: Use the grid defined by file "/d1/projects/MET/MET_test_data/unit_test/model_data/met_nc/nam/nam_2012040900_F036_APCP12.nc".

You can also use a named grid supported by MET:
https://met.readthedocs.io/en/develop/Users_Guide/appendixB.html#grids

attrs['grid'] = 'G218'

Or a grid specification string:
https://met.readthedocs.io/en/develop/Users_Guide/appendixB.html#grid-specification-strings

attrs['grid'] = 'lambert 614 428 12.190 -133.459 -95.0 12.19058 6367.47 25.0 N'

All of these options are outlined in the MET user guide Appendix F for Python embedding:
https://met.readthedocs.io/en/develop/Users_Guide/appendixF.html#pyembed-grid-attrs

I will respond separately in your comment from yesterday. Let me know if this information about providing grid metadata is sufficient, or whether there may be an additional approach that would be helpful for users that we could consider. Thanks!

[justintsu]

Still haven't heard back from the MET team yet but I dug up some more information on this in my debugging. I set point_stat debug verbosity to 4 so that it would print out the rejected observations while running the point observation PYTHON embedding routine in interactive mode. For some reason, I believe that some part of point_stat is setting my longitude to negative.

My python routine reads in raobs from a IODA file, organizes it in a dataframe and then converts it in a MET compatible dictionary called met_point_obs.

The point_stat debug statement prints this out:

**DEBUG 4: For rh/pre_000500 versus rh/P500, skipping observation off the grid where (x, y) = (140, 57) and grid (nx, ny) = (121, 121):
DEBUG 4:    raob DN717 20230330_000000 51.18000 _-16.59000_ -9999.00000 rh 300.00000 -9999.00000 NAVGEM 56.20312**

It makes sense that this observation was rejected because the lvl (300) does not match up with the 500 mb lvl of the forecast. However, the longitudes (italicized) of the observation point_stat registers and the observation I pass in do not match up (see below)

From Python:
['raob', 'DN717', '20230330_000000', 51.18000030517578, **_16.59000015258789_**, 'n/a', 'rh', 300.0, 'n/a', 'NAVGEM', 56.203121185302734]

Note that the 4th index in the list is a positive 16.59 but point_stat thinks this is a negative 16.59.

Another separate error I am noticing is in my forecast grid being passed into point_stat via attrs. when I run the python embedding function on its own, the grid I pass into point_stat is:

In [4]: grid
Out[4]:
{'name': 'COAMPS Test Grid',
 'nx': 121,
 'ny': 121,
 'type': 'Polar Stereographic',
 'hemisphere': 'N',
 'scale_lat': 60.0,
 'lat_pin': 71.0492,
 'lon_pin': 318.2404,
 'x_pin': 61.0,
 'y_pin': 61.0,
 'lon_orient': 317.789,
 'd_km': 45.0,
 'r_km': 6361.2}

The grid that point_stat tells me it interprets from python embedding is:

DEBUG 4: Grid information:
DEBUG 4:    Projection: Stereographic
DEBUG 4:    Nx: 121
DEBUG 4:    Ny: 121
DEBUG 4:    IsNorthHemisphere: true
DEBUG 4:    Lon_orient: 42.211
DEBUG 4:    Bx: 60.653
DEBUG 4:    By: 105.024
DEBUG 4:    Alpha: 263.7814
DEBUG 4:

Notice that it takes the complement of my lon_orient, which is not the same (it should be -42.211, if point_stat conforms to a [-180 180]. Everything else seems to pass correctly (although, i am not sure what Bx, By and Alpha are).

[DanielAdriaansen]

Hi @justintsu. MET expects longitudes to be provided in degrees east, so a hypothetical longitude somewhere over the United States would be provided as negative degrees east, whereas a hypothetical longitude over India would be provided as positive degrees east. This is mentioned in the MET user guide here:

https://met.readthedocs.io/en/develop/Users_Guide/appendixA.html#q-what-are-met-s-conventions-for-latitude-longitude-azimuth-and-bearing-angles

"MET considers north latitude and east longitude positive."

However, due to the original implementation design of the MET Earth geometry code, internally MET actually uses the opposite reference system (west longitude positive, east longitude negative). We recognize this is confusing and is something we'd like to change to provide more clarity to the user. Thus, you may see longitudes reported with the opposite sign in various debug statements in the code as you provided it because of how MET uses it internally.

You can always perform tests using the MET plot_point_obs tool to ensure MET plots your radiosonde site on the Earth where you expect it, or the MET plot_data_plane to plot your COAMPS data using Python embedding to ensure MET is placing your grid on the Earth where you expect using your Python embedding script and grid definition. Let me know if you'd like a sample command (or commands) for that, and let us know if you have any other questions!

Thanks again for your patience.

[justintsu]

Thanks for the response, can you provide me with an example using plot_point_obs and plot_data_plane if I need to use python embedding to pass in the COAMPS field/observation file? I've updated my python routines to ensure my lon is on a -180 to 180 span but I'd like to plot my obs and my forecast field and print this out to an image file to check and compare the latitude/longitudes. Can you offer more insight into how MET converts latitude and longitude inputs into grid points? For example, in the debug statement below:

DEBUG 4: raob 10304 20230405_060000 52.17000 -7.10000 -9999.00000 rh 150.00000 -9999.00000 NAVGEM 1.23108
DEBUG 4: For rh/pre_000700 versus rh/P700, skipping observation off the grid where (x, y) = (129, 46) and grid (nx, ny) = (121, 121):

How does MET convert 52.17N, -7E into (129,46)? On paper, this latitude longitude pair should be within my grid but it seems to fall outside the bounds (121, 121)

[DanielAdriaansen]

Hi Justin,

Just a note of clarification- this observation:

DEBUG 4: raob 10304 20230405_060000 52.17000 -7.10000 -9999.00000 rh 150.00000 -9999.00000 NAVGEM 1.23108

Is actually over Germany, because longitude is +7.10E, and thus it is off your grid. MET reports it as -7.10 because it uses degrees west internally, but requires degrees east as input. Therefore, the DEBUG statement there is printing the longitude in the units used internally not the units we request users to provide their longitude in. Indeed, an observation of -7.10E would be on your grid, but it would be printed in a DEBUG statement as +7.10 (degrees west). So behind the scenes in the code MET is:

1. Reading longitude in degrees east
2. Converting to degrees west (multiply by -1)
3. Convert back to degrees east before writing output files (different than DEBUG statements) by multiplying by -1 again

Let's use plot_point_obs to help demonstrate this. In units of degrees east, the site on your grid is over Ireland (-7.10E). Here is a Python embedding script for a single observation to plot this observation with plot_point_obs:

import pandas as pd
obs_string = 'raob 10304 20230405_060000 52.17000 -7.10000 -9999.00000 rh 0 0 NAVGEM 100.0'
obs_vals = obs_string.split()
df = pd.DataFrame([obs_vals])
df.columns=['typ','sid','vld','lat','lon','elv','var','lvl','hgt','qc','obs']
df.astype({'typ':'string','sid':'string','vld':'string','lat':'float','lon':'float','elv':'float','var':'string','lvl':'float','hgt':'float','qc':'string','obs':'float'})
point_data = df.values.tolist()

Using this command:

plot_point_obs PYTHON_NUMPY=test_python_obs.py test.ps -v 4

Produces an observation over Ireland:

But what about point_stat? Here is the DEBUG output from using that same Python embedding script with point_stat:

DEBUG 2: Searching 1 observations from 1 messages.
DEBUG 9: For TEST/L0 versus rh/L0, for observation type raob, over region FULL, for interpolation method NEAREST(1), using observation:
DEBUG 9:    raob 10304 20230405_060000 52.17000 7.10000 -9999.00000 rh 0.00000 0.00000 NAVGEM 100.00000
DEBUG 2: 
DEBUG 2: --------------------------------------------------------------------------------
DEBUG 2: 
DEBUG 2: Processing TEST/L0 versus rh/L0, for observation type raob, over region FULL, for interpolation method NEAREST(1), using 1 matched pairs.

Notice the sign of the longitude is different? I was providing the same data via Python but the debug statement reports a different sign, even though plot_point_obs places it over Ireland where it should be, and point_stat reports correctly that it's on your grid (it finds a matched pair).

I realize this is very confusing, so please ask for more clarification if this example doesn't help.

Regarding using plot_data_plane to visualize your gridded data, here is a sample Python embedding script using your grid definition, but populated with a constant field of 1 at each grid cell:

import numpy as np

grid_attrs = {'name': 'COAMPS Test Grid',
 'nx': 121,
 'ny': 121,
 'type': 'Polar Stereographic',
 'hemisphere': 'N',
 'scale_lat': 60.0,
 'lat_pin': 71.0492,
 'lon_pin': 318.2404,
 'x_pin': 61.0,
 'y_pin': 61.0,
 'lon_orient': 317.789,
 'd_km': 45.0,
 'r_km': 6361.2}

attrs = {}
attrs['valid'] = '20230405_060000'
attrs['init'] = attrs['valid']
attrs['lead'] = '000000'
attrs['accum'] = '000000'
attrs['name'] = 'TEST'
attrs['long_name'] = 'COAMPS Test'
attrs['level'] = 'L0'
attrs['units'] = 'None'
attrs['grid'] = grid_attrs

met_data = np.ones((121,121))

And here is the command to call plot_data_plane:

plot_data_plane PYTHON_NUMPY coamps_grid.ps 'name="test_python_grid.py";' -v 4

Which produces this grid:

Let us know if this helps answer your question, and if not please feel free to continue the discussion. Thanks!

[justintsu]

Thank you for the response. I've noted that the sign convention in the debug statement is in W convention even though MET requires longitude input to be in E. Thank you for providing me with examples on how to use plot_data_plane and plot_point_obs. Is there a setting with these functions to add latitude and longitude tick lines and labels on the plots?

[DanielAdriaansen]

No problem at all! Here are all the details about plot_data_plane:
https://met.readthedocs.io/en/develop/Users_Guide/plotting.html#plot-data-plane-usage

And plot_point_obs:
https://met.readthedocs.io/en/develop/Users_Guide/plotting.html#plot-point-obs-usage

I am not aware that there is an option for lat/lon lines, but there is limited support for colors and titles. The intent of these tools is mostly confirmation of placement of data spatially, and quick looks for confirming fields and obs look as a user would expect and less for more detailed visualization.

[justintsu]

great thanks. I ran this for both my forecast and obs for valid time 2023040212 for relative humidity

FORECAST:

OBS:

I only put half of the global obs. MET looks like it is rotating the forecast area (Arctic region above North America). Are these the exact grids that MET is using to find ob-fcst pairs and if so, does this imply that any matched pairs are not actually at the same lat/lon locations because the grids are transposed?

[DanielAdriaansen]

plot_data_plane is centering the view using the lon_orient value in your grid definition, which is 317.789, or ~ -42, or roughly the center of Greenland.

plot_point_obs uses a global GFS grid by default for placing observation points on the Earth/map, so maybe this is causing some confusion.

If you wish to plot the observation points on your forecast grid, you can point to one of your COAMPS test grid files, but MET must be able to read the file. If MET cannot read this file, you can also use a grid specification string (https://met.readthedocs.io/en/develop/Users_Guide/appendixB.html#grid-specification-strings), but you'll need to determine some additional metadata like the lat and lon of the lower left corner of the grid.

Does that help?