From e79449ce7b76be2e23cb2240d956bba9e963ec28 Mon Sep 17 00:00:00 2001
From: Lilli Freischem <lilli@freischem.biz>
Date: Mon, 27 Nov 2023 14:19:32 +0000
Subject: [PATCH] function to attach latlon coordinates to goes xarray dataset

---
 notebooks/goes-scripts/xy_to_latlon.py | 73 ++++++++++++++++++++++++++
 1 file changed, 73 insertions(+)
 create mode 100644 notebooks/goes-scripts/xy_to_latlon.py

diff --git a/notebooks/goes-scripts/xy_to_latlon.py b/notebooks/goes-scripts/xy_to_latlon.py
new file mode 100644
index 0000000..806d155
--- /dev/null
+++ b/notebooks/goes-scripts/xy_to_latlon.py
@@ -0,0 +1,73 @@
+import numpy as np
+
+'''
+script adapted from https://lsterzinger.medium.com/add-lat-lon-coordinates-to-goes-16-goes-17-l2-data-and-plot-with-cartopy-27f07879157f
+'''
+
+
+def calc_latlon(ds):
+    '''
+    Takes GOES dataset (one image) and computes latitude and 
+    longitude for each pixel using horizontal scan angles x
+    and vertical scan angles y.
+
+    Input:
+        ds  xarray.Dataset
+    Output:
+        ds  xarray.Dataset with lat and lon values added for each datapoint
+            and used as indeces.
+    '''
+    # The math for this function was taken from 
+    # https://makersportal.com/blog/2018/11/25/goes-r-satellite-latitude-and-longitude-grid-projection-algorithm 
+
+    x = ds.x
+    y = ds.y
+    goes_imager_projection = ds.goes_imager_projection
+    
+    x,y = np.meshgrid(x,y)
+    
+    r_eq = goes_imager_projection.attrs["semi_major_axis"] # earth radius at equator
+    r_pol = goes_imager_projection.attrs["semi_minor_axis"] # earth radius at pole
+    l_0 = goes_imager_projection.attrs["longitude_of_projection_origin"] * (np.pi/180) # lambda0
+    h_sat = goes_imager_projection.attrs["perspective_point_height"] # distance satellite to nearest equator surface point
+    H = r_eq + h_sat # distance satellite to earth centre
+    
+    a = np.sin(x)**2 + (np.cos(x)**2 * (np.cos(y)**2 + (r_eq**2 / r_pol**2) * np.sin(y)**2))
+    b = -2 * H * np.cos(x) * np.cos(y)
+    c = H**2 - r_eq**2
+    
+    r_s = (-b - np.sqrt(b**2 - 4*a*c))/(2*a)
+    
+    s_x = r_s * np.cos(x) * np.cos(y)
+    s_y = -r_s * np.sin(x)
+    s_z = r_s * np.cos(x) * np.sin(y)
+
+    # latitude and longitude
+    lat = np.arctan((r_eq**2 / r_pol**2) * (s_z / np.sqrt((H-s_x)**2 +s_y**2))) * (180/np.pi)
+    lon = (l_0 - np.arctan(s_y / (H-s_x))) * (180/np.pi)
+    
+    ds = ds.assign_coords({
+        "lat":(["y","x"],lat),
+        "lon":(["y","x"],lon)
+    })
+    ds.lat.attrs["units"] = "degrees_north"
+    ds.lon.attrs["units"] = "degrees_east"
+
+    return ds
+
+def get_xy_from_latlon(ds, lats, lons):
+    lat1, lat2 = lats
+    lon1, lon2 = lons
+
+    lat = ds.lat.data
+    lon = ds.lon.data
+    
+    x = ds.x.data
+    y = ds.y.data
+    
+    x,y = np.meshgrid(x,y)
+    
+    x = x[(lat >= lat1) & (lat <= lat2) & (lon >= lon1) & (lon <= lon2)]
+    y = y[(lat >= lat1) & (lat <= lat2) & (lon >= lon1) & (lon <= lon2)] 
+    
+    return ((min(x), max(x)), (min(y), max(y)))