pcolormesh: use shading='nearest' so no interpolation

example/plot_Efield.py

@@ -19,7 +19,7 @@ def plotVmaxx1it(ax: mpl.axes.Axes, V: xarray.DataArray) -> None:
         ax.set_ylabel("Potential [V]")
         ax.set_xlabel("mag. latitude [deg.]")
     elif V.ndim == 2:
-        hi = ax.pcolormesh(dat["mlon"], dat["mlat"], V, cmap="bwr")
+        hi = ax.pcolormesh(dat["mlon"], dat["mlat"], V, cmap="bwr", shading="nearest")
         ax.set_xlabel("mag. longitude [deg.]")
         ax.set_ylabel("mag. latitude [deg.]")
         ax.figure.colorbar(hi, ax=ax).set_label("potential [V]")

example/plotdata.py

@@ -14,7 +14,7 @@ def plotplasma(dat: xarray.Dataset):
             continue
 
         ax = figure().gca()
-        ax.pcolormesh(v.squeeze())
+        ax.pcolormesh(v.squeeze(), shading="nearest")
         ax.set_title(f"{k}  {t}")
 
 

scripts/plot_neutral.py

@@ -22,6 +22,6 @@
         ax = fg.gca()
 
         ax.set_title(k)
-        ax.pcolormesh(f[k][:].T)
+        ax.pcolormesh(f[k][:].transpose(), shading="nearest")
 
 plt.show()

src/gemini3d/compare/plot.py

@@ -132,16 +132,16 @@ def diff2d(A, B, name: str, fg, axs):
     bmin = min(A.min(), B.min())
     bmax = max(A.max(), B.max())
 
-    hi = axs[0].pcolormesh(A, cmap=cmap, vmin=bmin, vmax=bmax)
+    hi = axs[0].pcolormesh(A, cmap=cmap, vmin=bmin, vmax=bmax, shading="nearest")
     fg.colorbar(hi, ax=axs[0])
 
-    hi = axs[1].pcolormesh(B, cmap=cmap, vmin=bmin, vmax=bmax)
+    hi = axs[1].pcolormesh(B, cmap=cmap, vmin=bmin, vmax=bmax, shading="nearest")
     fg.colorbar(hi, ax=axs[1])
 
     dAB = A - B
     b = max(abs(dAB.min()), abs(dAB.max()))
 
-    hi = axs[2].pcolormesh(A - B, cmap="bwr", vmin=-b, vmax=b)
+    hi = axs[2].pcolormesh(A - B, cmap="bwr", vmin=-b, vmax=b, shading="nearest")
     fg.colorbar(hi, ax=axs[2])
 
     return abs(dAB).max().data

src/gemini3d/efield/__init__.py

@@ -196,7 +196,7 @@ def Efield_BCs(cfg: dict[str, T.Any], xg: dict[str, T.Any]) -> xarray.Dataset:
     # LEAVE THE SPATIAL AND TEMPORAL INTERPOLATION TO THE
     # FORTRAN CODE IN CASE DIFFERENT GRIDS NEED TO BE TRIED.
     # THE EFIELD DATA DO NOT TYPICALLY NEED TO BE SMOOTHED.
-    print("Size used for Efield input:  ",llon,llat);
+    print("Size used for Efield input:  ", llon, llat)
     write.Efield(E, cfg["E0dir"])
 
     return E

src/gemini3d/grid/tilted_dipole.py

@@ -500,7 +500,7 @@ def tilted_dipole3d_NUx2(cfg: dict[str, T.Any]) -> dict[str, T.Any]:
 
     ###########################################################################
     print(" Uniform grid lat lims:  ", xg["glat"].min(), xg["glat"].max())
-    
+
     print(" Generating non-uniform grid...")
     # Determine a target differential spacing based on user extents and number of
     #   grid points.
@@ -540,8 +540,7 @@ def tilted_dipole3d_NUx2(cfg: dict[str, T.Any]) -> dict[str, T.Any]:
     """
     xg = generate_tilted_dipole3d(q, pnew, phi)
 
-
-    print(" Shifting end L-shell by:  ", pnew[-3]-p[-3], pnew[-3], p[-3])
+    print(" Shifting end L-shell by:  ", pnew[-3] - p[-3], pnew[-3], p[-3])
     print(" Non-uniform grid lat lims:  ", xg["glat"].min(), xg["glat"].max())
-    
+
     return xg

src/gemini3d/magtools.py

@@ -95,7 +95,7 @@ def makegrid_full(
     write_grid: bool = False,
 ) -> dict[str, T.Any]:
     """
-     Make a field point to cover the entire mlat/mlon range for a given simulation grid
+    Make a field point to cover the entire mlat/mlon range for a given simulation grid
     """
 
     direc = Path(direc).expanduser()
@@ -119,10 +119,10 @@ def makegrid_full(
 
     # TABULATE THE SOURCE OR GRID CENTER LOCATION
     if "sourcemlon" not in cfg.keys():
-        thdist = xg["theta"].mean()
+        # thdist = xg["theta"].mean()
         phidist = xg["phi"].mean()
     else:
-        thdist = np.pi / 2 - np.radians(cfg["sourcemlat"])
+        # thdist = np.pi / 2 - np.radians(cfg["sourcemlat"])
         # zenith angle of source location
         phidist = np.radians(cfg["sourcemlon"])
 
@@ -134,14 +134,14 @@ def makegrid_full(
     # Computer a buffer region around the grid
     thmin = xg["theta"].min()
     thmax = xg["theta"].max()
-    dtheta=thmax-thmin
-    thmin=thmin-dtheta/5
-    thmax=thmax+dtheta/5
+    dtheta = thmax - thmin
+    thmin = thmin - dtheta / 5
+    thmax = thmax + dtheta / 5
     phimin = xg["phi"].min()
     phimax = xg["phi"].max()
-    dphi=phimax-phimin
-    phimin=phimin-dphi/5
-    phimax=phimax+dtheta/5
+    dphi = phimax - phimin
+    phimin = phimin - dphi / 5
+    phimax = phimax + dtheta / 5
 
     theta = np.linspace(thmin, thmax, ltheta)
     phi = phidist if flag2D else np.linspace(phimin, phimax, lphi)
@@ -169,7 +169,6 @@ def makegrid_full(
     return mag
 
 
-
 def magframe(
     filename: str | Path,
     *,
@@ -260,8 +259,8 @@ def magframe(
     with h5py.File(filename, "r") as f:
         for k in {"Br", "Btheta", "Bphi"}:
             if flatlist:
-                dat[k] = f[f"/magfields/{k}"][:]                
+                dat[k] = f[f"/magfields/{k}"][:]
             else:
                 dat[k] = f[f"/magfields/{k}"][:].reshape((lr, ltheta, lphi))
-                
+
     return dat

src/gemini3d/model.py

@@ -164,7 +164,7 @@ def equilibrium(cfg: dict[str, T.Any]):
 
     # %% Equilibrium input generation
     dat = equilibrium_state(cfg, xg)
-    
+
     write.state(cfg["indat_file"], dat)
 
 

src/gemini3d/particles/__init__.py

@@ -50,10 +50,10 @@ def particles_BCs(cfg: dict[str, T.Any], xg: dict[str, T.Any]):
     else:
         Qfunc = str2func("gemini3d.particles.gaussian2d")
 
-    Qtmp=Qfunc(pg, cfg["Qprecip"], cfg["Qprecip_background"])
+    Qtmp = Qfunc(pg, cfg["Qprecip"], cfg["Qprecip_background"])
 
     for i in range(i_on, i_off):
-        pg["Q"][i, :, :] = Qtmp[i,:,:]
+        pg["Q"][i, :, :] = Qtmp[i, :, :]
         pg["E0"][i, :, :] = cfg["E0precip"]
 
     assert np.isfinite(pg["Q"]).all(), "Q flux must be finite"

src/gemini3d/patch/plot.py

@@ -92,6 +92,7 @@ def patch(indir: Path, var: set[str]) -> None:
                     dat[k][ix1, :, :].transpose(),
                     vmin=clim[k][0],
                     vmax=clim[k][1],
+                    shading="nearest",
                 )
                 ax.text(
                     x2[x2.size // 2],

src/gemini3d/plot/cartesian.py

@@ -171,7 +171,7 @@ def plot_interp(
             # hack for pcolormesh to put labels in center of pixel
             wl = kwargs["wavelength"] + [""]
             hi = ax.pcolormesh(
-                np.arange(len(wl)), yp / 1e3, f(yp)[:, i].T, shading="nearest"
+                np.arange(len(wl)), yp / 1e3, f(yp)[:, i].transpose(), shading="nearest"
             )
             ax.set_xticks(np.arange(len(wl)) + 0.5)
             ax.set_xticklabels(wl)