add inverse transforms

modules/util/geos_nav.py

@@ -34,7 +34,11 @@ r_eq = 6378.1690  # km
 # LOFF = 1.53644345e-01
 
 # GOES ------------------------------
-h = 42164.16 # barycentric height, km
+h = 42164.16  # barycentric height, km
+invf = 298.2572221  # inverse flattening
+f = 1.0/invf
+fp = 1.0/((1.0-f)*(1.0-f))
+d = h*h - r_eq*r_eq
 scan_geom = 'GOES'
 sub_lon = -75.0
 sub_lon *= DEG_TO_RAD
@@ -60,6 +64,13 @@ LOFF = 0.151844
 # 65536 = 2^16
 
 
+def goes_to_geos(lamda_goes, theta_goes):
+    theta_geos = np.asin(np.sin(theta_goes) * np.cos(lamda_goes))
+    lamda_geos = np.atan(np.tan(lamda_goes) / np.cos(theta_goes))
+
+    return lamda_geos, theta_geos
+
+
 class GEOSNavigation:
     def __init__(self, sub_lon=-75.0, barycentric_height=42164.16, scan_geom='GOES', CFAC=5.6E-05, LFAC=-5.6E-05, COFF=-0.151844, LOFF=0.151844, num_elems=5424, num_lines=5424):
         self.sub_lon = sub_lon * DEG_TO_RAD
@@ -153,6 +164,14 @@ class GEOSNavigation:
         ll = np.where(np.isnan(theta), -1, ll)
         return cc, ll
 
+    def lc_to_earth_s(self, cc, ll):
+        x = cc * self.CFAC + self.COFF
+        y = ll * self.LFAC + self.LOFF
+
+        lon, lat = self.sat_to_earth_s(x, y)
+
+        return lon, lat
+
     def earth_to_indexs(self, lons, lats, len_x):
         num = lons.shape[0]
         idxs = []
@@ -167,6 +186,30 @@ class GEOSNavigation:
 
         return idxs
 
+    def sat_to_earth_s(self, x, y):
+        if self.scan_geom == 'GOES':
+            x, y = goes_to_geos(x, y)
+
+        c1 = (h * np.cos(x) * np.cos(y)) * (h * np.cos(x) * np.cos(y))
+        c2 = (np.cos(y) * np.cos(y) + fp * np.sin(y) * np.sin(y)) * d
+
+        c1 = np.where(c1 < c2, np.nan, c1)
+        c2 = np.where(c1 < c2, np.nan, c2)
+
+        s_d = np.sqrt(c1 - c2)
+        s_n = (h * np.cos(x) * np.cos(y) - s_d) / (np.cos(y) * np.cos(y) + fp * np.sin(y) * np.sin(y))
+        s_1 = h - s_n * np.cos(x) * np.cos(y)
+        s_2 = s_n * np.sin(x) * np.cos(y)
+        s_3 = s_n * np.sin(y)
+        s_xy = np.sqrt(s_1 * s_1 + s_2 * s_2)
+
+        geographic_lon = np.atan(s_2 / s_1) + self.sub_lon
+        geographic_lat = np.atan(fp * (s_3 / s_xy))
+
+        geographic_lon *= RAD_TO_DEG
+        geographic_lat *= RAD_TO_DEG
+
+        return geographic_lon, geographic_lat
 
 # def compute_scale_offset(lon_a, lat_a, col_a, line_a, lon_b, lat_b, col_b, line_b):
 #     lamda_a, theta_a = earth_to_sat(lon_a, lat_a)