Open blaylockbk opened 9 months ago
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)
Very cool. Would love to see a PR if you have the time.
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.