From 2ad230e73a5f473f1b4a06b65db4053c580e42b9 Mon Sep 17 00:00:00 2001
From: Nils Brueggemann <nils.brueggemann@mpimet.mpg.de>
Date: Mon, 10 Mar 2025 15:36:01 +0100
Subject: [PATCH] pyic_view: Improve handling of projection and data ranges.

---
 scripts/pyic_view.py | 111 ++++++++++++++++++++++++++++++-------------
 1 file changed, 77 insertions(+), 34 deletions(-)

diff --git a/scripts/pyic_view.py b/scripts/pyic_view.py
index 00718e1..7fb0e3e 100755
--- a/scripts/pyic_view.py
+++ b/scripts/pyic_view.py
@@ -53,6 +53,27 @@ def generate_axes(asp, generate_figure=True):
     cax.set_position([x0+axw0+daxcax, y0, caxw, axh0])
     return fig, ax, cax
 
+def get_xlim_ylim(lon_reg, lat_reg, proj, transformer):
+    if proj!='None':
+        #x_reg, y_reg = transformer.transform(
+        #    lon_reg, lat_reg,
+        #    direction='FORWARD',
+        #)
+        #xtest = np.linspace(x_reg[0], x_reg[1], 9)
+        #ytest = np.linspace(y_reg[0], y_reg[1], 9)
+        lon_test = np.linspace(lon_reg[0], lon_reg[1], 9)
+        lat_test = np.linspace(lat_reg[0], lat_reg[1], 9)
+        Lon_test, Lat_test = np.meshgrid(lon_test, lat_test)
+        Xtest, Ytest = transformer.transform(
+            Lon_test, Lat_test,
+            direction='FORWARD',
+        )
+        xlim = np.array([Xtest.min(), Xtest.max()])
+        ylim = np.array([Ytest.min(), Ytest.max()])
+    else:
+        xlim, ylim = lon_reg, lat_reg
+    return xlim, ylim
+
 def str_to_array(string):
   string = string.replace(' ', '')
   array = np.array(string.split(','), dtype=float)
@@ -138,11 +159,17 @@ class view(object):
         }
         self.font_size = 6*self.fig_size_fac
         self.res = 0.3
+        self.lon_reg = np.array([-180., 180.])
+        self.lat_reg = np.array([-90., 90.])
+        self.xlim = self.lon_reg
+        self.ylim = self.lat_reg
         self.it = 0
         self.iz = 0
         self.proj = self.proj_dict[list(self.proj_dict)[0]]
         if self.proj!="None":
             self.transformer = Proj.from_pipeline(self.proj)
+        else:
+            self.transformer = "None"
 
         # Opean data set
         self.load_data()
@@ -281,11 +308,11 @@ class view(object):
         reset_zoom_button.pack(side="left", padx=5)
 
         # proj entry
-        print('Setup proj dropdown')
+        self.message('Setup proj dropdown')
         res_menu = ttk.Combobox(frame3, textvariable=self.selected_proj, 
             values=list(self.proj_dict), state="readonly")
         res_menu.pack(side="left", padx=5)
-        res_menu.bind("<<ComboboxSelected>>", self.make_new_axis)
+        res_menu.bind("<<ComboboxSelected>>", self.update_projection)
 
         # --- end tk objects
 
@@ -301,7 +328,7 @@ class view(object):
     # reset zoom
     def reset_zoom(self):
         self.selected_res.set("0.3")
-        self.lon_lat_reg_tk.set("-180,180,-90,90")
+        self.set_default_lon_lat_reg()
         self.make_new_axis()
         self.update_data()
 
@@ -440,15 +467,8 @@ class view(object):
             adjust_axlims=False,
             )
         # set ax limits
-        if self.proj=="+proj=stere +lat_0=90 +lon_0=0":
-            self.ax.set_xlim([-4660515.349812048,  4660515.349812048])
-            self.ax.set_ylim([-4658959.2511977535, 4658959.2511977535])
-        elif self.proj=="+proj=stere +lat_0=-90 +lon_0=0":
-            self.ax.set_xlim([-5965970.154575175, 5965970.154575175])
-            self.ax.set_ylim([-5963978.177895851, 5963978.177895851])
-        else:
-            self.ax.set_xlim(self.X.min(), self.X.max())
-            self.ax.set_ylim(self.Y.min(), self.Y.max())
+        self.ax.set_xlim(self.xlim)
+        self.ax.set_ylim(self.ylim)
         # set ax labels
         if self.proj=="None" or self.proj=="+proj=latlong":
             #self.ax.set_xticks(np.arange(-180.,180.,45.))
@@ -476,9 +496,6 @@ class view(object):
         self.res = float(self.selected_res.get())
         self.proj = self.proj_dict[self.selected_proj.get()]
 
-        if self.proj!="None":
-            self.transformer = Proj.from_pipeline(self.proj)
-
         try:
           self.ax.remove()
           self.cax.remove()
@@ -487,17 +504,10 @@ class view(object):
 
         self.update_lon_lat_reg()
 
-        if self.proj=="+proj=stere +lat_0=90 +lon_0=0":
-            asp = 1.0
-        elif self.proj=="+proj=stere +lat_0=-90 +lon_0=0":
-            asp = 1.0
-        elif self.proj=="+proj=eqearth":
-            asp = 0.4867169753874043
-        elif self.proj=="+proj=moll":
-            asp = 0.5
-        else:
-            asp = (self.lat_reg[1]-self.lat_reg[0])/(self.lon_reg[1]-self.lon_reg[0])
+        asp = (self.ylim[1]-self.ylim[0])/(self.xlim[1]-self.xlim[0])
         self.fig, self.ax, self.cax = generate_axes(asp, generate_figure=False)
+        self.ax.set_xlim(self.xlim)
+        self.ax.set_ylim(self.ylim)
 
         #print('------')
         #print(self.fig.get_size_inches())
@@ -506,19 +516,52 @@ class view(object):
         self.plot_data()
         self.canvas.draw()
 
+    def set_default_lon_lat_reg(self):
+        if self.proj=="+proj=stere +lat_0=90 +lon_0=0":
+            #self.xlim = [-4660515.349812048,  4660515.349812048]
+            #self.ylim = [-4658959.2511977535, 4658959.2511977535]
+            self.lon_reg = [-180, 180]
+            self.lat_reg = [50, 90]
+        elif self.proj=="+proj=stere +lat_0=-90 +lon_0=0":
+            #self.xlim = [-5965970.154575175, 5965970.154575175]
+            #self.ylim = [-5963978.177895851, 5963978.177895851]
+            self.lon_reg = [-180, 180]
+            self.lat_reg = [-90, -40]
+        else:
+            self.lon_reg = [-180, 180]
+            self.lat_reg = [-90, 90]
+        self.lon_lat_reg_tk = tk.StringVar(value=f"{self.lon_reg[0]},{self.lon_reg[1]},{self.lat_reg[0]},{self.lat_reg[1]}")
+        #self.lon_lat_reg_tk.set(lon_lat_reg_str)
+        self.xlim, self.ylim = get_xlim_ylim(
+            self.lon_reg, self.lat_reg, self.proj, self.transformer)
+        print(f'set_default_lon_lat_reg:')
+        print(self.lon_reg, self.lat_reg, self.xlim, self.ylim)
+        
+    def update_projection(self, *args):
+        self.proj = self.proj_dict[self.selected_proj.get()]
+        if self.proj!="None":
+            self.transformer = Proj.from_pipeline(self.proj)
+        self.set_default_lon_lat_reg()
+        self.make_new_axis()
+
     def update_data_range(self):
         xlim = self.ax.get_xlim()
         ylim = self.ax.get_ylim()
-        if self.proj!='None':
-            lon_lim, lat_lim = self.transformer.transform(
-                xlim, ylim,
-                direction='INVERSE',
-            ) 
+        if self.proj=="None":
+            self.lon_reg, self.lat_lim = xlim, ylim
+            self.xlim, self.ylim = xlim, ylim
         else:
-            lon_lim = xlim
-            lat_lim = ylim
-        lon_lat_reg_str = f'{lon_lim[0]:3g},{lon_lim[1]:3g},{lat_lim[0]:3g},{lat_lim[1]:3g}'
-        self.lon_lat_reg_tk.set(lon_lat_reg_str)
+            # FIXME: This is likely not working well since data
+            # interp_to_rectgrid_xr (which is currently workin in the back)
+            # can only operate with simple regional limitis but not as complicated
+            # as some projections require it to be
+            lon_lim, lat_lim = self.transformer.transform(
+                xlim, ylim, direction='INVERSE')
+            self.lon_reg, self.lat_reg = lon_lim, lat_lim
+            self.xlim, self.ylim = xlim, ylim
+            self.message('::: Warning: "Updating data range" with more complex projections can lead to strange results. :::')
+
+        self.lon_lat_reg_tk = tk.StringVar(value=f"{self.lon_reg[0]},{self.lon_reg[1]},{self.lat_reg[0]},{self.lat_reg[1]}")
         self.make_new_axis()
 
     def increase_slider(self, slider):
-- 
GitLab