Getting point values from TEMPO L2 ozone profile V04 in R

[kodiologist]

I've been looking at the recently released TEMPO_O3PROF_L2_V04 product. Given a NetCDF file (such as TEMPO_O3PROF_L2_V04_20250921T210803Z_S012G02.nc) and some points, I'd like to get the ozone values (`/product/ozone_profile`) from the bottom layer (layer 24) for the given points. Do you have any suggestions on how to do this, particularly in R? If you don't know a method for R, a method for Python would be helpful.

My first attempt used the `stars` package, but there may be a bug preventing this from working correctly: https://github.com/r-spatial/stars/issues/755

[dkaufas]

Hi there,

While I don't have a direct R solution to offer, here's a Python approach using `xarray` that might be helpful:

Code: Select all

import xarray as xr

# Open the NetCDF file
filepath = "TEMPO_O3PROF_L2_V04_20250921T210803Z_S012G02.nc"
dt = xr.open_datatree(filepath)

# Extract ozone profile values from layer 24 (index 23, since Python uses 0-based indexing)
single_layer_values = dt["product/ozone_profile"].isel(layer=23)

For R, you might also try the `ncdf4` or `RNetCDF` packages while waiting for the `stars` issue to be resolved.

Best,
Danny

[kodiologist]

That appears to read in an array well enough (which in R can be done with the `ncdf4` package), but without any geometry. I'm looking for a way to get values corresponding to points on land, such as longitude -73.953626, latitude 40.789693.

[dkaufas]

Hi,

Ahh yes, extracting data at specific spatial points requires additional steps since TEMPO Level-2 data use a 2D irregular grid. Here's a Python approach that should work:

Code: Select all

import numpy as np
import xarray as xr

# Define the filepath
filepath = "TEMPO_O3PROF_L2_V04_20250921T210803Z_S012G02.nc"

# Define the Point Of Interest (lon, lat).
POI = np.array((-73.953626, 40.789693))

# Open the NetCDF file and geolocation arrays.
dt = xr.open_datatree(filepath)
lon_array = dt["geolocation/longitude"].values
lat_array = dt["geolocation/latitude"].values

# Transform longitude and latitude arrays to an N x 2 array of coordinate pairs.
coordinate_pairs = np.column_stack((lon_array.flatten(), lat_array.flatten()))

# Find the data point closest to the point of interest.
idx_min_flat = np.sum( 
    (coordinate_pairs - POI) ** 2,  axis=1  # Calculate squared Euclidean distances
).argmin(axis=0)  # Find the index of the minimum squared distance

# Convert the flattened index to row and column indices.
row_index, col_index = np.unravel_index(idx_min_flat, lon_array.shape)

print(f"flat index: {idx_min_flat}\n"
      f"mirror_step and xtrack indices: {row_index}, {col_index}\n"
      f"coordinate closest to POI: {coordinate_pairs[idx_min_flat]}")

# Extract ozone profile values from layer 24 (index 23, since Python uses 0-based indexing)
value_at_POI = dt["product/ozone_profile"].isel(
    {"mirror_step": row_index, "xtrack": col_index, "layer": 23}
).values

print(f"value_at_POI: {value_at_POI}")

Note that this approach finds the nearest data point using Euclidean distance. If your point of interest falls outside the area of the data granule, I think you'd still get a result, but it would be from the edge of the granule. You could consider adding a distance check to ensure the point identified is within a reasonable distance threshold.

Another note: For even greater precision with geolocation, you could use the `geolocation/longitude_bounds` and `geolocation/latitude_bounds` variables to create polygons and do point-in-polygon checks, though I don't currently have an example for those computations.

Best,
Danny

[kodiologist]

Yeah, I think that the cell with the closest centroid is not always the cell that contains the point of interest. Previously, in a project in which I was using TROPOMI observations, I was constructing a polygon mesh for each file and using point-in-polygon, with only minor optimizations. It was really slow when I had a lot of points and a lot of files. You can't just do the computation only once per point, either, because the quadrilateral mesh changes slightly between overpasses of the same areas. Perhaps there are no shortcuts here. Thanks.