0

I have an xarray dataset that contains a data variable 'tmpc_obs' that was created using the interpolate_to_grid function within MetPy. I would like to be able to interpolate specific point values from lat lon pairs using xarray.interp, but I need assign a crs and then assign x and y values. MetPy has options to do this, but I am not having much luck getting it to work.

I created map_crs using cartopy and is what I will transform my lat lon points to, and then create the interpolated field. latlon_crs is the crs for lat lon points.


map_crs = ccrs.AlbersEqualArea(central_longitude=-79.0, 
                               central_latitude=35.5, 
                               standard_parallels=(34.33,36.16),
                               )

latlon_crs = ccrs.PlateCarree()

Transform the points to the map_crs and then interpolate the point values to a grid.

xp, yp, _ = map_crs.transform_points(latlon_crs, obs['lon'], obs['lat']).T
obsx, obsy, tmpc = interpolate_to_grid(xp, yp, obs['tmpc'].values,
                                      interp_type='rbf',
                                      rbf_func='thin_plate',
                                      hres=5000)

This turns the 3 arrays generated into an xarray dataset. x_p and y_p are x and y coordinates that are transformed to the map crs.

obs_ds = xr.Dataset(
    data_vars=dict(
        tmpc_obs=(["x", "y"], tmpc),
    ),
    coords=dict(
        x_p=(["x", "y"], obsx),
        y_p=(["x", "y"], obsy),
    ),
)

It's not possible however to interpolate point values at this time using the following code block as it will only give you nan's.

vals = []
for x, y in zip(stat['lon'], stat['lat']):
    x_t, y_t = map_crs.transform_point(x, y, src_crs=latlon_crs)
    tmp_c = obs_ds["tmpc_obs"].interp(x=x_t, y=y_t, method='cubic').values.item()
    tmp_c = float("{:.2f}".format(tmp_c))
    vals.append(tmp_c)

I believe you need to add y and x to the coordinates of the dataset that reference the dimensions. I am having issues doing this. I think I should use the same crs info from the mapcrs, but I don't know what to use for semi major and semi minor axis. Whether I use the axes or not, I get an error, and it doesnt attach the x and y to the coordinates, it only attaches the metpy_crs

# Assign the grid mapping to be able to index specific lat/lon pairs
obs_ds = obs_ds.metpy.assign_crs({
    'semi_major_axis': ???,
    'semi_minor_axis': ???,
    'grid_mapping_name': 'albers_conical_equal_area',
    "standard_parallel": [34.33,36.16],
    "latitude_of_projection_origin": 35.5,
    "longitude_of_central_meridian": -79.0,
}).metpy.assign_y_x()

Error:

UserWarning: No y and x coordinates assigned since horizontal coordinates were not found warnings.warn('No y and x coordinates assigned since horizontal coordinates '

<xarray.Dataset>
Dimensions:    (x: 57, y: 75)
Coordinates:
    x_p        (x, y) float64 -1.963e+05 -1.913e+05 ... 1.692e+05 1.742e+05
    y_p        (x, y) float64 -1.663e+05 -1.663e+05 ... 1.137e+05 1.137e+05
    metpy_crs  object Projection: albers_conical_equal_area
Dimensions without coordinates: x, y
Data variables:
    tmpc       (x, y) float64 0.0618 0.07065 0.07945 ... 0.08191 0.09538 0.1064
Cane89
  • 3
  • 3

0 Answers0