Add an accessor to extract data nearest a particular latitude/longitude point

[blaylockbk]

This could be a useful feature for some people interested in timeseries of values (e.g., cloud top temperature) at a one or more locations.

Ideally, a user should be able to provide the desired lat/lon point, then the method would calculate the location in the geostationary projection, determine the nearest grid point in the GOES data, and return the value at that point. Probably could use cartopy to transform the lat/lon coordinate to geostationary projection coordinate.

[blaylockbk]

Possible starting point...

• cartopy.crs.CRS — cartopy 0.22.0 documentation (scitools.org.uk)
• Parse_cf in Xarray Integration (metpy.xarray) — MetPy 1.6 (unidata.github.io)
• Xarray Integration (metpy.xarray) — MetPy 1.6 (unidata.github.io)

[ibarlet]

I implemented this feature in an iPython notebook using the Xarray preprocessing function. I'll try to make some time to wrap it into a class and make a PR but in case I don't get around to it and anyone is looking for this quickly here's a simplistic starting point.

import xarray as xr
import numpy as np
from functools import partial

from goes2go import GOES, config

def lat_lon_to_scan_angles(latitude, longitude, goes_imager_projection, decimal_degrees=True):
    """
    Convert latitude and longitude to ABI scan angle coordinates.
    This converts geodetic coordinates into to the GRS80 elliopsoid as used by GOES-R.
    https://www.goes-r.gov/users/docs/PUG-L1b-vol3.pdf (section 5.1.2.8.2)
    """
    if decimal_degrees:
        latitude = np.radians(latitude)
        longitude = np.radians(longitude)
    
    r_pol = goes_imager_projection.semi_minor_axis
    r_eq = goes_imager_projection.semi_major_axis
    H = goes_imager_projection.perspective_point_height + r_eq
    e = 0.0818191910435  # eccentricity of GRD 1980 ellipsoid
    lambda_0 = np.radians(goes_imager_projection.longitude_of_projection_origin) # Always in degrees from the GOES data

    phi_c = np.arctan((r_pol**2/r_eq**2)*np.tan(latitude))  # Geocentric latitude
    r_c = r_pol/(np.sqrt(1-(e**2*np.cos(phi_c)**2)))  # Geocentric distance to the point on the ellipsoid

    s_x = H - r_c*np.cos(phi_c)*np.cos(longitude-lambda_0)
    s_y = -r_c*np.cos(phi_c)*np.sin(longitude-lambda_0)
    s_z = r_c*np.sin(phi_c)

    # Confirm location is visible from the satellite
    if (H*(H-s_x) < (s_y**2 + (r_eq**2/r_pol**2)*s_z**2)).any():
        raise ValueError("One or more points not visible by satellite")
    
    y = np.arctan(s_z/s_x)
    x = np.arcsin(-s_y / (np.sqrt(s_x**2+s_y**2+s_z**2)))

    return x, y  # Always returned in radians

def _preprocess(ds, target_lat, target_lon, decimal_degrees=True):
    x_target, y_target = lat_lon_to_scan_angles(target_lat, target_lon, ds["goes_imager_projection"], decimal_degrees)
    return ds.sel(x=x_target, y=y_target, method="nearest")

target_latitude = 38.897957
target_longitude = -77.036560

G = GOES(satellite=goes16,
         product="ABI-L2-TPW",
         domain="C").timerange(
             start='2024-04-08 17:00',
             end='2024-04-08 20:00',
             return_as='filelist'
         )

partial_func = partial(_preprocess, target_lat=target_latitude, target_lon=target_longitude, decimal_degrees=True)
preprocessed_ds = xr.open_mfdataset([str(config['timerange']['save_dir']) + "/" + f for f in G['file'].to_list()],
                  concat_dim='t',
                  combine='nested',
                  preprocess=partial_func)

print(preprocessed_ds.TPW)

[blaylockbk]

Very cool. Would love to see a PR if you have the time.