Merge branch 'master' into `non_lon_lat`

FindNetCDF_Fortran.cmake

@@ -5,22 +5,25 @@ if(TARGET NetCDF_Fortran::netcdff)
 else()
   # Find NetCDF dependency:
 
-  option(use_find_netcdf_module "Use the find module for NetCDF")
-
-  unset(extraArgs)
+  option(FIND_PACKAGE_PREFER_MODULE_netCDF
+    "Use directly the find module for NetCDF")
 
   if(${CMAKE_FIND_PACKAGE_NAME}_FIND_QUIETLY)
-    list(APPEND extraArgs QUIET)
+    set(maybe_quiet QUIET)
   endif()
   
   if(${CMAKE_FIND_PACKAGE_NAME}_FIND_REQUIRED)
-    list(APPEND extraArgs REQUIRED)
+    set(maybe_required REQUIRED)
   endif()
   
-  if(use_find_netcdf_module)
-    find_package(netCDF ${extraArgs})
+  if(FIND_PACKAGE_PREFER_MODULE_netCDF)
+    find_package(netCDF ${maybe_quiet} ${maybe_required})
   else()
-    find_package(netCDF CONFIG ${extraArgs})
+    find_package(netCDF CONFIG ${maybe_quiet})
+
+    if(NOT netCDF_FOUND)
+      find_package(netCDF ${maybe_quiet} ${maybe_required})
+    endif()
   endif()
 
   #-

Inst_eddies/Analysis/eddy_dump.py

@@ -15,25 +15,31 @@ import util_eddies
 
 parser = argparse.ArgumentParser()
 parser.add_argument("directory", help = "containing the three shapefiles")
+parser.add_argument("-i", "--ishape", type = int)
 args = parser.parse_args()
 
-reply = input("days_1950, eddy_index = ? ").split(",")
-days_1950 = int(reply[0])
-eddy_index = int(reply[1])
-
 handler = util_eddies.open_shpc(args.directory)
 
-# Find ishape:
-if handler["ishape_last"] is None:
-    print("ishape_last not found. Iterating through shapes...")
-    
-    for ishape, rec in enumerate(handler["readers"]["extremum"].iterRecords()):
-        if rec["days_1950"] == days_1950 and rec["eddy_index"] == eddy_index:
-            break
-    else:
-        sys.exit("Eddy with this date and index not found")
+if args.ishape:
+    ishape = args.ishape
 else:
-    ishape = util_eddies.comp_ishape(handler, days_1950, eddy_index)
+    reply = input("days_1950, eddy_index = ? ").split(",")
+    days_1950 = int(reply[0])
+    eddy_index = int(reply[1])
+
+    # Find ishape:
+    if handler["ishape_last"] is None:
+        print("ishape_last not found. Iterating through shapes...")
+
+        for ishape, rec \
+            in enumerate(handler["readers"]["extremum"].iterRecords()):
+            if rec["days_1950"] == days_1950 \
+               and rec["eddy_index"] == eddy_index:
+                break
+        else:
+            sys.exit("Eddy with this date and index not found")
+    else:
+        ishape = util_eddies.comp_ishape(handler, days_1950, eddy_index)
         
 for layer in ["extremum", "outermost_contour", "max_speed_contour", "centroid"]:
     if layer in handler["readers"]:

Inst_eddies/Analysis/plot_eddy_contours.py

@@ -131,15 +131,10 @@ def snapshot(ax, ishape_list, handler, *, dashed = False, light = False,
                 if not light: lines[0].set_marker("o")
                 if dashed: lines[0].set_linestyle("dashed")
 
-def plot_grid_bb(shpc_dir, ax):
+def plot_grid_bb(grid_nml, ax):
     """Grid bounding box."""
-    
-    file = path.join(shpc_dir, "grid_nml.txt")
-    try:
-        grid_nml = f90nml.read(file)["grid_nml"]
-    except FileNotFoundError:
-        print("grid_nml.txt not found. Will not plot bounding box.")
-    else:
+
+    if grid_nml is not None:
         step = grid_nml["STEP_DEG"]
         rect = patches.Rectangle(grid_nml["corner_deg"],
                                  (grid_nml["nlon"] - 1) * step[0],
@@ -153,10 +148,6 @@ if __name__ == "__main__":
     import netCDF4
 
     parser = argparse.ArgumentParser()
-    parser.add_argument("-v", "--velocity", help = "plot velocity field",
-                        action = "store_true")
-    parser.add_argument("-s", "--scale", default = 20, type = float,
-                        help = "scale of arrows for the velocity field")
     parser.add_argument("-d", "--date", type = int,
                         help = "date, as days since 1950-1-1")
     parser.add_argument("-g", "--grid", help = "plot grid",
@@ -175,17 +166,24 @@ if __name__ == "__main__":
                         help = "Save file to specified format")
     args = parser.parse_args()
 
-    if args.grid or args.velocity:
-        with netCDF4.Dataset("h.nc") as f:
-            if "lon" in f.variables:
-                lon = "lon"
-                lat = "lat"
-            else:
-                lon = "longitude"
-                lat = "latitude"
-                
-            longitude = f[lon][:]
-            latitude = f[lat][:]
+    if args.grid or args.window is None:
+        file = path.join(args.shpc_dir[0], "grid_nml.txt")
+    
+        try:
+            grid_nml = f90nml.read(file)["grid_nml"]
+        except FileNotFoundError:
+            print("grid_nml.txt not found. Will not plot bounding box.")
+            grid_nml = None
+
+    if args.grid:
+        width = (grid_nml["NLON"] - 1) * grid_nml["STEP_DEG"][0]
+        height = (grid_nml["NLat"] - 1) * grid_nml["STEP_DEG"][1]
+        longitude = np.linspace(grid_nml["CORNER_DEG"][0],
+                                grid_nml["CORNER_DEG"][0] + width,
+                                grid_nml["NLON"])
+        latitude = np.linspace(grid_nml["CORNER_DEG"][1],
+                               grid_nml["CORNER_DEG"][1] + height,
+                               grid_nml["NLat"])
     
     if args.window is not None:
         llcrnrlon, llcrnrlat, urcrnrlon, urcrnrlat = args.window
@@ -193,14 +191,14 @@ if __name__ == "__main__":
         if urcrnrlon - llcrnrlon > 360:
             sys.exit("bad values of urcrnrlon and llcrnrlon")
 
-        if args.grid or args.velocity:
+        if args.grid:
             longitude += np.ceil((llcrnrlon - longitude) / 360) * 360
             # (in [llcrnrlon, llcrnrlon + 2 pi[)
 
             lon_mask = longitude <= urcrnrlon
             lat_mask = np.logical_and(latitude >= llcrnrlat,
                                       latitude <= urcrnrlat)
-    elif args.grid or args.velocity:
+    elif args.grid:
         lon_mask = np.ones(len(longitude), dtype = bool)
         lat_mask = np.ones(len(latitude), dtype = bool)
 
@@ -216,19 +214,7 @@ if __name__ == "__main__":
         ax.plot(lon_2d.reshape(-1), lat_2d.reshape(-1), transform = src_crs,
                 marker = "+", color = "gray", linestyle = "None")
 
-    if args.window is None: plot_grid_bb(args.shpc_dir[0], ax)
-
-    if args.velocity:
-        with netCDF4.Dataset("uv.nc") as f:
-            quiver_return = ax.quiver(longitude[lon_mask],
-                                      latitude[lat_mask],
-                                      f["ugos"][0, lat_mask][:, lon_mask],
-                                      f["vgos"][0, lat_mask][:, lon_mask],
-                                      scale = args.scale,
-                                      scale_units = "width",
-                                      transform = src_crs)
-        plt.quiverkey(quiver_return, 0.9, 0.9, 1, r"1 m s$^{-1}$",
-                      coordinates = "figure")
+    if args.window is None: plot_grid_bb(grid_nml, ax)
 
     for shpc_dir in args.shpc_dir:
         handler = util_eddies.open_shpc(shpc_dir)

Inst_eddies/Analysis/plot_velocity.py

+#!/usr/bin/env python3
+
+import numpy as np
+import cartopy.crs as ccrs
+import sys
+import matplotlib.pyplot as plt
+import argparse
+import netCDF4
+
+parser = argparse.ArgumentParser()
+parser.add_argument("-s", "--scale", default = 20, type = float,
+                    help = "scale of arrows for the velocity field")
+parser.add_argument("-w", "--window", help = "choose a limited plot window",
+                    type = float, nargs = 4,
+                    metavar = ("llcrnrlon", "llcrnrlat", "urcrnrlon",
+                               "urcrnrlat"))
+parser.add_argument("--save", metavar = "format",
+                    help = "Save file to specified format")
+parser.add_argument("-u", "--undefined", action = "store_true",
+                    help = "plot points where velocity is not defined")
+parser.add_argument("input_file", help = "NetCDF file containing velocity")
+args = parser.parse_args()
+
+with netCDF4.Dataset(args.input_file) as f:
+    if "lon" in f.variables:
+        lon = "lon"
+        lat = "lat"
+    else:
+        lon = "longitude"
+        lat = "latitude"
+
+    longitude = f[lon][:]
+    latitude = f[lat][:]
+
+if args.window is not None:
+    llcrnrlon, llcrnrlat, urcrnrlon, urcrnrlat = args.window
+
+    if urcrnrlon - llcrnrlon > 360:
+        sys.exit("bad values of urcrnrlon and llcrnrlon")
+
+    longitude += np.ceil((llcrnrlon - longitude) / 360) * 360
+    # (in [llcrnrlon, llcrnrlon + 2 pi[)
+
+    lon_mask = longitude <= urcrnrlon
+    lat_mask = np.logical_and(latitude >= llcrnrlat,
+                              latitude <= urcrnrlat)
+else:
+    lon_mask = np.ones(len(longitude), dtype = bool)
+    lat_mask = np.ones(len(latitude), dtype = bool)
+
+fig = plt.figure()
+src_crs = ccrs.PlateCarree()
+projection = ccrs.PlateCarree()
+##projection = ccrs.NorthPolarStereo()
+ax = plt.axes(projection = projection)
+
+with netCDF4.Dataset(args.input_file) as f:
+    if args.undefined:
+        undef_velocity = np.logical_or(f["ugos"][0].mask, f["vgos"][0].mask)
+        lon_2d, lat_2d = np.meshgrid(longitude[lon_mask],
+                                     latitude[lat_mask])
+        ax.plot(lon_2d[undef_velocity].reshape(-1),
+                lat_2d[undef_velocity].reshape(-1), transform = src_crs,
+                marker = "*", color = "violet", linestyle = "None")
+    else:
+        quiver_return = ax.quiver(longitude[lon_mask],
+                                  latitude[lat_mask],
+                                  f["ugos"][0, lat_mask][:, lon_mask],
+                                  f["vgos"][0, lat_mask][:, lon_mask],
+                                  scale = args.scale,
+                                  scale_units = "width",
+                                  transform = src_crs)
+        plt.quiverkey(quiver_return, 0.9, 0.9, 1, r"1 m s$^{-1}$",
+                      coordinates = "figure")
+
+ax.gridlines(draw_labels = True)
+ax.coastlines()
+
+if args.save:
+    plt.savefig(f"plot_velocity.{args.save}")
+    print(f'Created "plot_velocity.{args.save}".')
+else:
+    plt.show()

Inst_eddies/Analysis/tests.json

@@ -59,18 +59,12 @@
 	"description": "With window option."
     },
     {
-	"title": "Plot_eddy_contours_vel",
+	"title": "Plot_velocity",
 	"command":
 	[
-	    "$src_dir/Inst_eddies/Analysis/plot_eddy_contours.py",
-	    "$tests_old_dir/Extraction_eddies_region_4/SHPC_anti",
-	    "--velocity", "--save=png"
-	],
-	"description": "With velocity option.",
-	"required":
-	[
-	    ["$src_dir/Inst_eddies/Tests/Input/h_region_4.nc", "h.nc"],
-	    ["$src_dir/Inst_eddies/Tests/Input/uv_region_4.nc", "uv.nc"]
+	    "$src_dir/Inst_eddies/Analysis/plot_velocity.py",
+	    "$src_dir/Inst_eddies/Tests/Input/uv_region_4.nc",
+	    "--save=png"
 	]
     },
     {

Inst_eddies/Tests/CMakeLists.txt

@@ -100,9 +100,22 @@ add_executable(test_set_max_speed
   ${PROJECT_SOURCE_DIR}/Common/read_field_indices.F90
   ${PROJECT_SOURCE_DIR}/Common/read_eddy.f90
   ${PROJECT_SOURCE_DIR}/Common/shpc_create.f90
-  ${PROJECT_SOURCE_DIR}/Common/write_eddy.f90)
+  ${PROJECT_SOURCE_DIR}/Common/write_eddy.f90 get_var.f90)
 
 target_link_libraries(test_set_max_speed Geometry::geometry
   Numer_Rec_95::numer_rec_95 NetCDF95::netcdf95 Shapelib_03::shapelib_03
   Contour_531::contour_531 Jumble::jumble NR_util::nr_util
   NetCDF_Fortran::netcdff gpc_f)
+
+foreach(my_target IN ITEMS test_get_1_outerm test_set_all_outerm
+    test_good_contour test_inside_4 test_mean_speed
+    test_max_speed_contour_ssh test_nearby_extr test_local_extrema
+    test_set_max_speed)
+
+  set_target_properties(${my_target} PROPERTIES Fortran_MODULE_DIRECTORY
+    ${CMAKE_CURRENT_BINARY_DIR}/${my_target}_modules)
+
+  target_include_directories(${my_target} PRIVATE
+    ${PROJECT_BINARY_DIR}/${my_target}_modules>)
+
+endforeach()

Inst_eddies/Tests/Input/SHPC_degenerate/orientation.txt

+anticyclones

Inst_eddies/Tests/short_tests.json

@@ -270,6 +270,27 @@
 	    ]
 	]
     },
+    {
+        "required":
+	[
+            ["$src_dir/Inst_eddies/Tests/Input/degenerated_SSH.nc", "h.nc"],
+            ["$src_dir/Inst_eddies/Tests/Input/degenerated_SSH.nc", "uv.nc"],
+            [
+		"$src_dir/Inst_eddies/Tests/Input/empty_outside_points.csv",
+		"outside_points.csv"
+	    ]
+        ],
+        "title": "Set_max_speed_degenerate",
+	"description": "The SSH value of the extremum is exactly equal to the SSH value of the point next to the extremum, eastward. This points is selected as the point with maximum azimuthal speed. Then good_contour returns a null polyline.",
+        "commands":
+	[
+	    ["mkdir", "SHPC"],
+	    [
+		"$build_dir/Inst_eddies/test_set_max_speed",
+		"$src_dir/Inst_eddies/Tests/Input/SHPC_degenerate"
+	    ]
+	]
+    },
     {
         "create_file": ["main_nml.txt", "&main_nml min_amp = 0./\n"],
 	"input": "20454\n",

Inst_eddies/Tests/test_set_max_speed.f90

 program test_set_max_speed
 
   use, intrinsic:: ISO_FORTRAN_ENV
+  use, intrinsic:: ieee_arithmetic, only: ieee_value, IEEE_QUIET_NAN
+
 
   ! Libraries:
   use jumble, only: new_unit, count_lines, get_command_arg_dyn
@@ -10,6 +12,7 @@ program test_set_max_speed
   use nr_util, only: deg_to_rad
 
   use derived_types, only: eddy, shpc
+  use get_var_m, only: get_var
   use read_eddy_m, only: read_eddy
   use set_max_speed_m, only: set_max_speed
   use shpc_close_m, only: shpc_close
@@ -57,15 +60,16 @@ program test_set_max_speed
 
   print *, "Reading from uv.nc..."
   call nf95_open("uv.nc", nf90_nowrite, ncid)
-
   allocate(u(nlon, nlat), v(nlon, nlat))
-
-  call nf95_inq_varid(ncid, "ugos", varid)
-  call nf95_get_var(ncid, varid, u)
-
-  call nf95_inq_varid(ncid, "vgos", varid)
-  call nf95_get_var(ncid, varid, v)
-
+  call get_var(periodic = .false., max_rad_lon = 0, values = u, ncid = ncid, &
+       nlon = nlon, name = "ugos", &
+       new_fill_value = ieee_value(0., IEEE_QUIET_NAN))
+  call get_var(periodic = .false., max_rad_lon = 0, values = v, ncid = ncid, &
+       nlon = nlon, name = "vgos", &
+       new_fill_value = ieee_value(0., IEEE_QUIET_NAN))
+  ! (We will need quiet NaNs rather the original fill values for u and
+  ! v when we compute the max-speed contours and when we search the
+  ! ssh of max-speed contours.)
   call nf95_close(ncid)
 
   print *, "Reading from shapefiles..."

Inst_eddies/max_speed_contour_ssh.f90

@@ -33,7 +33,7 @@ contains
 
     !----------------------------------------------------------------------
 
-    ! Assuming that ascending indices in u and v mean correspond to
+    ! Assuming that ascending indices in u and v correspond to
     ! ascending longitude and latitude:
     v_azim(east, :) = v(ind_extr(1) + 1:ind_extr(1) + radius - 1, &
          ind_extr(2))
@@ -50,7 +50,12 @@ contains
        ! outermost contour as max-speed contour.
        max_speed_contour_ssh = missing_ssh
     else
-       l = maxloc(abs(sum(v_azim, dim = 1) / 4.), dim = 1)
+       if (radius >= 3) then
+          l = maxloc(abs(sum(v_azim, dim = 1) / 4.), dim = 1)
+       else
+          l = 1
+       end if
+       
        direction = maxloc(abs(v_azim(:, l)), dim = 1)
 
        select case (direction)