Longitude and latitude axes instead of indices in examine_ssh_values.py.

Tests/examine_ssh_values.py

@@ -3,26 +3,35 @@
 import netCDF4
 import matplotlib.pyplot as plt
 import numpy as np
+import math
 
 with netCDF4.Dataset("h.nc") as f:
     ssh = f.variables["adt"][:].squeeze()
+    longitude = f.variables["lon"][:]
+    latitude = f.variables["lat"][:]
+
+lon_center = float(input("longitude of center = ? "))
+lat_center = float(input("latitude of center = ? "))
+
+# Base 0 indices, assuming regular grid::
+i_center = round(float((lon_center - longitude[0])
+                       / (longitude[1] - longitude[0])))
+j_center = round(float((lat_center - latitude[0]) / (latitude[1] - latitude[0])))
+# (float to convert from numpy.float64, so that round produces an int)
 
-print("Give base 0 indices below:")
-i_min = int(input("imin = ? "))
-i_max = int(input("imax = ? "))
-j_min = int(input("jmin = ? "))
-j_max = int(input("jmax = ? "))
 level = float(input("level = ? "))
 
 plt.figure()
-plt.contour(ssh, (level,))
-i_array = np.arange(i_min,i_max + 1)
-j_array = np.arange(j_min,j_max + 1)
-i_array_2d, j_array_2d = np.meshgrid(i_array, j_array)
-plt.plot(i_array_2d.reshape(-1), j_array_2d.reshape(-1), "+")
-
-for i in i_array:
-    for j in j_array:
-        plt.annotate(ssh[j,i], (i, j))
+plt.contour(longitude, latitude, ssh, (level,))
+delta = 4
+lon_2d, lat_2d = np.meshgrid(longitude[i_center - delta:i_center + delta + 1],
+                             latitude[j_center - delta:j_center + delta + 1])
+plt.plot(lon_2d.reshape(-1), lat_2d.reshape(-1), "+")
+plt.plot(longitude[i_center], latitude[j_center], "o", markersize = 10,
+         fillstyle = "none")
+
+for i in np.arange(i_center - delta,i_center + delta + 1):
+    for j in np.arange(j_center - delta,j_center + delta + 1):
+        plt.annotate(ssh[j,i], (longitude[i], latitude[j]))
 
 plt.show()

Tests/plot_snapshot.py

@@ -29,11 +29,11 @@ lon_2d, lat_2d = np.meshgrid(longitude, latitude)
 plt.figure()
 plt.plot(lon_2d.reshape(-1), lat_2d.reshape(-1), "+", color="gray")
 
-reader = shapefile.Reader("extremum_1")
+reader_extr = shapefile.Reader("extremum_1")
 reader_outer = shapefile.Reader("outermost_contour_1")
 reader_m_s = shapefile.Reader("max_speed_contour_1")
 
-for shape_rec_extr, shape_outer, shape_m_s in zip(reader.iterShapeRecords(),
+for shape_rec_extr, shape_outer, shape_m_s in zip(reader_extr.iterShapeRecords(),
                                                   reader_outer.iterShapes(),
                                                   reader_m_s.iterShapes()):
     points = shape_rec_extr.shape.points[0]