From 8e08808c852af65e08e9e10a167dcce069bae694 Mon Sep 17 00:00:00 2001
From: Nils Brueggemann <nils.brueggemann@mpimet.mpg.de>
Date: Wed, 5 Mar 2025 08:16:11 +0100
Subject: [PATCH] pyic_view: Allows different projections.

---
 scripts/pyic_view.py | 74 +++++++++++++++++++++++++++++++++++++++-----
 1 file changed, 66 insertions(+), 8 deletions(-)

diff --git a/scripts/pyic_view.py b/scripts/pyic_view.py
index 85ae4d3..f4a7fc0 100755
--- a/scripts/pyic_view.py
+++ b/scripts/pyic_view.py
@@ -120,10 +120,21 @@ class view(object):
             "cmo.thermal", "cmo.haline", "cmo.curl",
         ]
         self.res_all = [1., 0.3, 0.1, 0.02]
+        self.proj_all = [
+          "None", 
+          "+proj=latlong",
+          "+proj=stere +lat_0=90 +lon_0=0",
+          "+proj=stere +lat_0=-90 +lon_0=0",
+          "+proj=eqearth",
+          "+proj=moll",
+        ]
         self.font_size = 6*self.fig_size_fac
         self.res = 0.3
         self.it = 0
         self.iz = 0
+        self.proj = self.proj_all[0]
+        if self.proj!="None":
+            self.transformer = Proj.from_pipeline(self.proj)
 
         # Opean data set
         self.load_data()
@@ -134,6 +145,7 @@ class view(object):
         self.color_limits = tk.StringVar(value="auto")  # Default color limits
         self.lon_lat_reg_tk = tk.StringVar(value="-180,180,-90,90")
         self.selected_res = tk.StringVar(value="0.3")
+        self.selected_proj = tk.StringVar(value="None")
 
         # Create figure and axis
         self.fig, self.ax, self.cax = generate_axes(asp=0.5)
@@ -167,11 +179,11 @@ class view(object):
         print('Setup sliders')
         self.slider_t = tk.Scale(root, from_=0, to=len(self.ds.time)-1, 
             orient="horizontal", label="time", command=self.update_plot)
-        self.slider_t.grid(row=1, column=0, sticky="ew")
+        self.slider_t.grid(row=1, column=0, columnspan=2, sticky="ew")
         
         self.slider_d = tk.Scale(root, from_=0, to=len(self.ds.depth)-1, 
             orient="horizontal", label="depth", command=self.update_plot)
-        self.slider_d.grid(row=1, column=1, sticky="ew")
+        self.slider_d.grid(row=1, column=2, columnspan=2, sticky="ew")
         
         # Create dropdown menus
         print('Setup var dropdown')
@@ -331,9 +343,27 @@ class view(object):
             self.fpath_ckdtree, self.lon_reg, self.lat_reg
         )
         self.Lon, self.Lat = np.meshgrid(self.dai.lon.data, self.dai.lat.data)
+        if self.proj=="None":
+            self.X, self.Y = self.Lon, self.Lat
+        else:
+            self.X, self.Y = self.transformer.transform(self.Lon, self.Lat, direction='FORWARD')
+            #valid = np.isfinite(self.X) & np.isfinite(self.Y)
+            #self.dai[valid==False] = np.nan
         # make plot
-        self.hm = pyic.shade(self.dai.lon, self.dai.lat, self.dai, 
+        valid = np.isfinite(self.X) & np.isfinite(self.Y)
+        self.hm = pyic.shade(
+            #self.X[valid], self.Y[valid], self.dai.data[valid], 
+            self.X, self.Y, self.dai.data, 
             ax=self.ax, cax=self.cax)
+        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())
         #pyic.plot_settings(self.ax, xlim=self.lon_reg, ylim=self.lat_reg)
         self.ax.set_facecolor('0.7')
         self.update_cmap()
@@ -346,7 +376,12 @@ class view(object):
         self.update_title()
 
     def make_new_axis(self, *args):
+
         self.res = float(self.selected_res.get())
+        self.proj = self.selected_proj.get()
+
+        if self.proj!="None":
+            self.transformer = Proj.from_pipeline(self.proj)
 
         try:
           self.ax.remove()
@@ -355,8 +390,21 @@ class view(object):
           pass
 
         self.update_lon_lat_reg()
-        #proj = ccrs.PlateCarree()
-        proj = None
+
+        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])
+        self.fig, self.ax, self.cax = generate_axes(asp, generate_figure=False)
+
+        ##proj = ccrs.PlateCarree()
+        #proj = None
 
         self.ax = self.fig.add_subplot(
           position=self.pos_ax, projection=proj)
@@ -465,12 +513,22 @@ class view(object):
     def on_click(self, event):
         # Avoid clicking outside the axes
         if event.xdata is not None and event.ydata is not None:  
-            lon0, lat0 = event.xdata, event.ydata
+            if self.proj!="None":
+                lon_click, lat_click = self.transformer.transform(
+                    event.xdata, event.ydata, 
+                    direction='INVERSE',
+                ) 
+            else:
+                lon_click, lat_click = event.xdata, event.ydata
             ind = np.argmin(
-                (self.Lon.flatten()-lon0)**2+(self.Lat.flatten()-lat0)**2
+                (self.Lon.flatten()-lon_click)**2+(self.Lat.flatten()-lat_click)**2
             )
             data_click = self.dai.data.flatten()[ind]
-            print(f"lon:{event.xdata:.2f}, lat: {event.ydata:.2f}, data: {data_click:.4f}")
+                
+            txt = f"lon:{lon_click:.2f}, lat: {lat_click:.2f}, data: {data_click:.4f}"
+            self.ht_point.set_text(txt)
+            print(txt)
+            self.canvas.draw()
 
 def main():
     import argparse 
-- 
GitLab