diff --git a/pyicon/pyicon_plotting.py b/pyicon/pyicon_plotting.py index 7d8839f32798eef300c416727c8ab1f967f0ac0e..0dd734d0060dcf9ac842e74efa4db5b0b9cbe75c 100644 --- a/pyicon/pyicon_plotting.py +++ b/pyicon/pyicon_plotting.py @@ -96,6 +96,10 @@ def hplot_base(IcD, IaV, clim='auto', cmap='viridis', cincr=-1., if projection=='none': ccrs_proj = None ccrs_transform = None + # Allow users to use cartopy projections directly + elif isinstance(projection, ccrs.Projection): + ccrs_proj = projection + ccrs_transform = ccrs.PlateCarree() elif projection=='RotatedPole': print ("generating cartopy projection with a rotated pole lon/lat = " + str(IcD.pol_lon) + "/" + str(IcD.pol_lat)) @@ -2048,10 +2052,11 @@ def plot(data, shade_proj = ccrs.PlateCarree() ccrs_proj = getattr(ccrs, projection) ccrs_proj = ccrs_proj(central_longitude=central_longitude) - else: + elif isinstance(projection, ccrs.Projection): shade_proj = ccrs.PlateCarree() ccrs_proj = projection - + else: + raise ValueError(f"`projection` must be a ccrs.Projection object or a string, not {type(projection)}") # --- rename dimensions if 'cell' in data.dims: data = data.rename(cell='ncells') diff --git a/pyicon/quickplots/pyicon_quickplots.py b/pyicon/quickplots/pyicon_quickplots.py index 3c24ad33813748979c9a001c8bcaa827296b1afb..de66471c90975b55ed784b17baecac515fe0370f 100644 --- a/pyicon/quickplots/pyicon_quickplots.py +++ b/pyicon/quickplots/pyicon_quickplots.py @@ -7,15 +7,6 @@ import numpy as np # --- reading data from netCDF4 import Dataset, num2date import datetime -# --- plotting -import matplotlib.pyplot as plt -import matplotlib -from matplotlib import ticker -#import my_toolbox as my -import cartopy -import cartopy.crs as ccrs -from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter -import cmocean # --- debugging from ipdb import set_trace as mybreak import pyicon as pyic @@ -124,12 +115,6 @@ def qp_hplot(fpath, var, IcD='none', depth=-1e33, iz=0, it=0, if units!='': IaV.units = units - # --- cartopy projection - if projection=='none': - ccrs_proj = None - else: - ccrs_proj = getattr(ccrs, projection)() - # --- do plotting (ax, cax, hm, diff --git a/scripts/pyic_view.py b/scripts/pyic_view.py index 00718e191199f94bbf6823493eaf7012036ca9d9..37e71a4f7342a6be97735edc522fa40790f0ca3c 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) @@ -106,7 +127,7 @@ def get_data(ds, var_name, it, iz, res, lon_reg, lat_reg): class view(object): def __init__(self, flist, path_grid, fig_size_fac=1.0): # Initialize Tkinter - print('setup TKinter') + self.message('setup TKinter') root = tk.Tk() root.title("pyicon view") root.geometry("1200x800") @@ -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() @@ -151,7 +178,7 @@ class view(object): self.selected_var = tk.StringVar(value=self.var_names[0]) self.selected_cmap = tk.StringVar(value=self.colormaps[0]) self.color_limits = tk.StringVar(value="auto") # Default color limits - self.lon_lat_reg_tk = tk.StringVar(value="-180,180,-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.save_fig_tk = tk.StringVar(value="./fig.pdf") self.selected_res = tk.StringVar(value="0.3") self.selected_proj = tk.StringVar(value="None") @@ -177,7 +204,7 @@ class view(object): self.canvas.mpl_connect("button_press_event", self.on_click) # --- 1st row - print('Setup sliders') + self.message('Setup sliders') pady_num = 8 frame1 = tk.Frame(root) frame1.pack(fill="x", pady=pady_num) @@ -213,7 +240,7 @@ class view(object): frame2.pack(fill="x", pady=pady_num) # Create dropdown menus - print('Setup var dropdown') + self.message('Setup var dropdown') var_menu = ttk.Combobox( frame2, textvariable=self.selected_var, values=list(self.ds.data_vars.keys()), state="readonly" @@ -222,14 +249,14 @@ class view(object): var_menu.bind("<<ComboboxSelected>>", self.update_data) # Cmap dropdown - print('Setup cmap dropdown') + self.message('Setup cmap dropdown') cmap_menu = ttk.Combobox(frame2, textvariable=self.selected_cmap, values=self.colormaps, state="readonly") cmap_menu.pack(side="left", padx=5) cmap_menu.bind("<<ComboboxSelected>>", self.update_cmap) # Color limit entry - print('Setup color limits') + self.message('Setup color limits') entry = tk.Entry(frame2, textvariable=self.color_limits) entry.pack(side="left", padx=5) entry.insert(0, "") # Default value @@ -252,14 +279,14 @@ class view(object): frame3.pack(fill="x", pady=pady_num) # res entry - print('Setup res dropdown') + self.message('Setup res dropdown') res_menu = ttk.Combobox(frame3, textvariable=self.selected_res, values=self.res_all, state="readonly") res_menu.pack(side="left", padx=5) res_menu.bind("<<ComboboxSelected>>", self.make_new_axis) # lon_lat_reg entry - print('Setup lon_reg') + self.message('Setup lon_reg') self.entry_lon_lat_reg = tk.Entry(frame3, textvariable=self.lon_lat_reg_tk) self.entry_lon_lat_reg.pack(side="left", padx=5) self.entry_lon_lat_reg.insert(0, "") # Default value @@ -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 @@ -294,21 +321,24 @@ class view(object): self.canvas.draw() # Start Tkinter loop - print('Go into mainloop') + self.message('Go into mainloop') root.mainloop() return + def message(self, message): + print(f'pyic_view: {message}') + # 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() # for saving the figure def save_figure(self, *args): fpath = self.save_fig_tk.get() - print(f'Saving figure {fpath}') + self.message(f'Saving figure {fpath}') #plt.savefig(fpath, dpi=300) plt.savefig(fpath, dpi=300, bbox_inches='tight') @@ -362,7 +392,7 @@ class view(object): self.canvas.mpl_disconnect(self.cid_motion) def load_data(self): - print('opening dataset') + self.message('opening dataset') mfdset_kwargs = dict( combine='nested', concat_dim='time', data_vars='minimal', coords='minimal', @@ -371,6 +401,8 @@ class view(object): ) + self.message('Data from these files is considered:') + self.message(self.flist) self.ds = xr.open_mfdataset( self.flist, **mfdset_kwargs, chunks=dict(time=1, depth=1, depth_2=1) @@ -389,7 +421,7 @@ class view(object): except: pass self.var_names = list(self.ds) - print(f"variables in data set: {self.var_names}") + self.message(f"variables in data set: {self.var_names}") self.var_name = self.var_names[0] # copy uuidOfHGrid to each variable as attribute @@ -413,7 +445,7 @@ class view(object): self.nzs[var] = 0 except: delvars.append(var) - print(f'Excluding the following variables from data set due to issues with dimension specification: {delvars}') + self.message(f'Excluding the following variables from data set due to issues with dimension specification: {delvars}') self.ds = self.ds.drop_vars(delvars) def plot_data(self): @@ -440,15 +472,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 +501,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 +509,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 +521,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): @@ -544,7 +592,7 @@ class view(object): # update depth slider self.slider_d.config(to=self.nzs[self.var_name]-1) - print(f'Updating data: {self.var_name}: it = {self.it}; iz = {self.iz}') + self.message(f'Updating data: {self.var_name}: it = {self.it}; iz = {self.iz}') # Get data and plot self.update_lon_lat_reg() @@ -578,7 +626,7 @@ class view(object): try: clim = self.get_clim(clim_str, self.dai) self.hm[0].set_clim(clim[0], clim[1]) - print(f'Updating clim to {tuple(map(float, clim))}') + self.message(f'Updating clim to {tuple(map(float, clim))}') self.canvas.draw() except ValueError: print(f'Invalid value for clim: {clim_str}') @@ -590,7 +638,7 @@ class view(object): def update_cmap(self, *args): # update cmap cmap = self.selected_cmap.get() - print(f"Updating cmap to {cmap}") + self.message(f"Updating cmap to {cmap}") if cmap.startswith('cmo'): cmap = cmap.split('.')[-1] cmap = getattr(cmocean.cm, cmap) @@ -604,7 +652,7 @@ class view(object): lon_lat_reg = str_to_array(lon_lat_reg_str) self.lon_reg = [lon_lat_reg[0], lon_lat_reg[1]] self.lat_reg = [lon_lat_reg[2], lon_lat_reg[3]] - print(f'Updating lon_reg to {tuple(map(float, self.lon_reg))} and lat_reg to {tuple(map(float, self.lat_reg))}') + self.message(f'Updating lon_reg to {tuple(map(float, self.lon_reg))} and lat_reg to {tuple(map(float, self.lat_reg))}') def get_clim(self, clim, data): # --- clim @@ -624,7 +672,7 @@ class view(object): def toggle_grid(self): if self.do_grid.get(): - print('Adding grid lines') + self.message('Adding grid lines') color='k' linewidth=0.5 if self.proj=="+proj=stere +lat_0=90 +lon_0=0": @@ -671,7 +719,7 @@ class view(object): self.hgs.append(hg) #print(f'nn = {nn}, {self.hgs}') else: - print('Removing grid lines') + self.message('Removing grid lines') for hg in self.hgs: hg.remove() self.canvas.draw() @@ -725,11 +773,7 @@ def main(): flist = iopts.fpath_data flist.sort() - print('Data from these files is considered:') - print(flist) - # Initial plot - print('Initialize plot') path_grid = 'auto' View = view(flist, path_grid=path_grid,