Skip to content
Snippets Groups Projects
Commit 0ce6b1b4 authored by Lionel GUEZ's avatar Lionel GUEZ
Browse files

If an outermost contour has exactly the minimum amplitude, it is

acceptable. So define flat_extr by a strict inequality in procedure
get_snapshot.

Procedure set_outermost_contour no longer has a dummy argument
noise_around. We no longer call outermost_possible_level so we do not
have an artificial discontinuity of results with and without maximum
amplitude. We abandon the idea of reducing the amplitude of initial
level_good compared to innermost_level: problem of consistence with
get_snapshot, and it does not seem worth the trouble, just abandon
those problematic extrema. Allow for null outermost_contour instead of
aborting.

In procedure get_snapshot, s%list_vis(i)%suff_amp is first defined
only if flat_extr(i). Also we are no longer sure that
set_outermost_contour finds an outermost contour so we have to test
this to define s%list_vis(i)%suff_amp. noise_around is now defined
only for flat_extr(i) and s%list_vis(i)%suff_amp. Also, since we may
not find an outermost contour even if not flat_extr(i), we update
s%extr_map after the second call to set_outermost_contour.

In procedure local_extrema, we use the mask diff_center instead of
mask_center. So we no longer need procedure construct_mask_center. In
procedure local_extrema, we ignore one in two adjacent degenerate
extrema instead of aborting.

In plot_snapshot.py, color extrema and allow for null outermost
contour.

Synthesize output in test_local_extrema.py.
parent 2c5912b1
No related branches found
No related tags found
No related merge requests found
......@@ -17,8 +17,9 @@ module derived_types
logical cyclone
type(ssh_contour) outermost_contour, max_speed_contour
real max_speed ! average on max_speed_contour, in m s-1
logical suff_amp ! sufficient amplitude of ssh, between extremum
! and outermost contour
logical suff_amp ! outermost_contour found and sufficient
! amplitude of ssh between extremum and
! outermost contour
integer:: delta_in = huge(0)
! Maximum difference in time index where there is a direct
......@@ -36,9 +37,10 @@ module derived_types
integer number_vis_eddies ! number of visible eddies
integer, allocatable:: extr_map(:, :) ! (nlon, nlat) At a point
! of extremum SSH with sufficient amplitude: identification
! number or this extremum. At a point of extremum SSH with
! insufficient amplitude: the opposite of the identification
! of extremum SSH with outermost_contour found, of sufficient
! amplitude: identification number or this extremum. At a point
! of extremum SSH with no outermost_contour or outermost_contour
! of insufficient amplitude: the opposite of the identification
! number or this extremum. 0 at other points.
integer, allocatable:: ind_extr(:, :) ! (2, number_vis_eddies)
......
......@@ -93,7 +93,7 @@ contains
ind_targ_extr(2, s%number_vis_eddies), &
flat_extr(s%number_vis_eddies), noise_around(s%number_vis_eddies))
do i = 1, s%number_vis_eddies
forall (i = 1:s%number_vis_eddies)
s%list_vis(i)%coord_extr = corner + (s%ind_extr(:, i) - 1) * step
s%list_vis(i)%ssh_extremum = ssh(s%ind_extr(1, i), s%ind_extr(2, i))
s%list_vis(i)%cyclone = cyclone(i)
......@@ -103,12 +103,12 @@ contains
urc(:, i) = min(s%ind_extr(:, i) + max_radius, [nlon, nlat])
corner_window(:, i) = corner + (llc(:, i) - 1) * step
ind_targ_extr(:, i) = s%ind_extr(:, i) - llc(:, i) + 1
end do
end forall
if (min_amp == 0.) then
flat_extr = .false.
else
flat_extr = abs(innermost_level - s%list_vis%ssh_extremum) <= min_amp
flat_extr = abs(innermost_level - s%list_vis%ssh_extremum) < min_amp
end if
do i = 1, s%number_vis_eddies
......@@ -117,41 +117,58 @@ contains
innermost_level(i), &
s%extr_map(llc(1, i):urc(1, i), llc(2, i):urc(2, i)), &
ssh(llc(1, i):urc(1, i), llc(2, i):urc(2, i)), &
corner_window(:, i), step, noise_around = .false.)
s%list_vis(i)%suff_amp = abs(s%list_vis(i)%outermost_contour%ssh &
- s%list_vis(i)%ssh_extremum) >= min_amp
else
s%list_vis(i)%suff_amp = .true.
corner_window(:, i), step)
if (s%list_vis(i)%outermost_contour%n_points == 0) then
s%list_vis(i)%suff_amp = .false.
else
s%list_vis(i)%suff_amp &
= abs(s%list_vis(i)%outermost_contour%ssh &
- s%list_vis(i)%ssh_extremum) >= min_amp
end if
end if
end do
! We must modify s%extr_map in this loop, separate from the
! previous loop, so we do not influence the batch of
! set_outermost_contour.
do i = 1, s%number_vis_eddies
if (.not. s%list_vis(i)%suff_amp) s%extr_map(s%ind_extr(1, i), &
s%ind_extr(2, i)) = - s%extr_map(s%ind_extr(1, i), &
s%ind_extr(2, i))
if (flat_extr(i) .and. .not. s%list_vis(i)%suff_amp) &
s%extr_map(s%ind_extr(1, i), s%ind_extr(2, i)) &
= - s%extr_map(s%ind_extr(1, i), s%ind_extr(2, i))
end do
! (no forall because of gfortran internal compiler error)
! Define noise_around:
if (min_amp == 0.) then
noise_around = .false.
else
do i = 1, s%number_vis_eddies
if (s%list_vis(i)%suff_amp) noise_around(i) &
= any(s%extr_map(llc(1, i):urc(1, i), llc(2, i):urc(2, i)) &
< 0)
end do
end if
if (min_amp /= 0.) &
forall &
(i = 1:s%number_vis_eddies, &
flat_extr(i) .and. s%list_vis(i)%suff_amp) &
noise_around(i) &
= any(s%extr_map(llc(1, i):urc(1, i), llc(2, i):urc(2, i)) < 0)
do i = 1, s%number_vis_eddies
if (s%list_vis(i)%suff_amp .and. noise_around(i) &
.or. .not. flat_extr(i)) &
call set_outermost_contour(s%list_vis(i), ind_targ_extr(:, i), &
innermost_level(i), &
s%extr_map(llc(1, i):urc(1, i), llc(2, i):urc(2, i)), &
ssh(llc(1, i):urc(1, i), llc(2, i):urc(2, i)), &
corner_window(:, i), step, noise_around(i))
.or. .not. flat_extr(i)) then
call set_outermost_contour(s%list_vis(i), ind_targ_extr(:, i), &
innermost_level(i), &
s%extr_map(llc(1, i):urc(1, i), llc(2, i):urc(2, i)), &
ssh(llc(1, i):urc(1, i), llc(2, i):urc(2, i)), &
corner_window(:, i), step)
s%list_vis(i)%suff_amp &
= s%list_vis(i)%outermost_contour%n_points /= 0
end if
end do
! We must modify s%extr_map in this loop, separate from the
! previous loop, so we do not influence the batch of
! set_outermost_contour.
do i = 1, s%number_vis_eddies
if (.not. flat_extr(i) .and. .not. s%list_vis(i)%suff_amp) &
s%extr_map(s%ind_extr(1, i), s%ind_extr(2, i)) &
= - s%extr_map(s%ind_extr(1, i), s%ind_extr(2, i))
end do
! (no forall because of gfortran internal compiler error)
! Done with outermost contours, now let us take care of the
! max-speed contours.
......
......@@ -2,8 +2,6 @@ module local_extrema_m
implicit none
logical, save, private:: mask_center(3, 3)
contains
subroutine local_extrema(field, extr_map, ind_extr, innermost_level, &
......@@ -33,12 +31,13 @@ contains
logical, intent(out), allocatable:: local_min(:) ! (n_extr)
! Local:
integer n_extr, i, j, k, m, n
integer n_extr, i, j, k, m, n, count_diff
real innermost_map(size(field, 1), size(field, 2)) ! (m, n) map of
! levels for innermost contours, defined only for extremum points
real max_around, min_around
logical diff_center(3, 3) ! values different from the center value
!--------------------------------------------------------------------
......@@ -55,27 +54,51 @@ contains
do j = 2, n - 1
do i = 2, m - 1
max_around = maxval(field(i - 1:i + 1, j - 1:j + 1), mask_center)
diff_center = field(i - 1:i + 1, j - 1:j + 1) /= field(i, j)
count_diff = count(diff_center)
if (field(i, j) > max_around) then
n_extr = n_extr + 1
extr_map(i, j) = n_extr
innermost_map(i, j) = max_around
if (count_diff == 0) then
extr_map(i, j) = 0
else
min_around = minval(field(i - 1:i + 1, j - 1:j + 1), mask_center)
if (field(i, j) == max_around .or. field(i, j) == min_around) then
print *, "local_extrema: we require strict extrema"
print *, "field(i = ", i, ", j = ", j, ") = ", field(i, j)
print *, "max_around = ", max_around
print *, "min_around = ", min_around
stop 1
else if (field(i, j) < min_around) then
n_extr = n_extr + 1
extr_map(i, j) = - n_extr
innermost_map(i, j) = min_around
! {count_diff >= 1}
max_around = maxval(field(i - 1:i + 1, j - 1:j + 1), diff_center)
if (field(i, j) > max_around) then
! {local maximum, possibly degenerate}
if (count_diff <= 7 .and. any([extr_map(i - 1:i + 1, j - 1), &
extr_map(i - 1, j)] > 0)) then
! {degenerate maximum and one of previously examined
! values is a maximum, so field(i, j) equals one of
! those previously examined values}
extr_map(i, j) = 0
else
n_extr = n_extr + 1
extr_map(i, j) = n_extr
innermost_map(i, j) = max_around
end if
else
extr_map(i, j) = 0
! {field(i, j) < max_around}
min_around = minval(field(i - 1:i + 1, j - 1:j + 1), &
diff_center)
if (field(i, j) < min_around) then
! {local minimum, possibly degenerate}
if (count_diff <= 7 .and. &
any([extr_map(i - 1:i + 1, j - 1), &
extr_map(i - 1, j)] < 0)) then
! {degenerate minimum and one of previously examined
! values is a minimum, so field(i, j) equals one of
! those previously examined values}
extr_map(i, j) = 0
else
n_extr = n_extr + 1
extr_map(i, j) = - n_extr
innermost_map(i, j) = min_around
end if
else
! {min_around < field(i, j) < max_around}
extr_map(i, j) = 0
end if
end if
end if
end do
......@@ -89,12 +112,12 @@ contains
i = 1
j = 2
do k = 1, n_extr
do
! next (i, j)
i = i + 1
if (i == m) then
! Go to next column
i = 2
......@@ -122,13 +145,4 @@ contains
end subroutine local_extrema
!***********************************************************************
subroutine construct_mask_center
mask_center = .true.
mask_center(2, 2) = .false.
end subroutine construct_mask_center
end module local_extrema_m
......@@ -5,7 +5,7 @@ module set_outermost_contour_m
contains
subroutine set_outermost_contour(e, ind_targ_extr, innermost_level, &
extr_map, ssh, corner, step, noise_around)
extr_map, ssh, corner, step)
! Defines e%outermost_contour. Method: dichotomy on level of ssh.
......@@ -36,12 +36,10 @@ contains
real, intent(in):: step(:) ! (2)
! longitude and latitude steps, in rad
logical, intent(in):: noise_around
! Local:
integer i ! identifying number of the eddy
real, allocatable:: nearby_extr(:, :) ! (2, :) longitude and
! latitude, in rad, of all the significant extrema, except the
! target extremum
......@@ -53,56 +51,52 @@ contains
!-----------------------------------------------------------------
i = extr_map(ind_targ_extr(1), ind_targ_extr(2))
if (noise_around) then
nearby_extr = convert_to_reg_coord(argwhere(extr_map > 0 &
.and. extr_map /= i), corner, step)
level_try = merge(maxval(ssh), minval(ssh), e%cyclone)
nearby_extr = convert_to_reg_coord(argwhere(extr_map > 0 &
.and. extr_map /= i), corner, step)
e%outermost_contour%polyline = good_contour(corner, step, ssh, &
innermost_level, e%coord_extr, nearby_extr)
if (e%outermost_contour%n_points == 0) then
e%outermost_contour%ssh = e%ssh_extremum
e%outermost_contour%area = 0.
else
nearby_extr = convert_to_reg_coord(pack_indices(extr_map, &
excluded = [0, i]), corner, step)
level_try = outermost_possible_level(ssh, ind_targ_extr, e%cyclone)
end if
level_good = e%ssh_extremum &
+ (innermost_level - e%ssh_extremum) * (1. - 100. * epsilon(0.))
e%outermost_contour%polyline = good_contour(corner, step, ssh, level_good, &
e%coord_extr, nearby_extr)
call assert(e%outermost_contour%n_points /= 0, &
"set_outermost_contour innermost_level")
call assert((level_try - e%ssh_extremum) / (level_good - e%ssh_extremum) &
> 1., "set_outermost_contour level_try")
tentative_contour = good_contour(corner, step, ssh, level_try, &
e%coord_extr, nearby_extr)
if (tentative_contour%n_points /= 0) then
e%outermost_contour%polyline = tentative_contour
level_good = level_try
else
! {no good contour at level_try}
level_bad = level_try
do while (abs(level_bad - level_good) > precision)
level_try = (level_good + level_bad) / 2.
tentative_contour = good_contour(corner, step, ssh, level_try, &
e%coord_extr, nearby_extr)
if (tentative_contour%n_points /= 0) then
level_good = level_try
e%outermost_contour%polyline = tentative_contour
else
level_bad = level_try
end if
! {e%outermost_contour%polyline == good_contour(. . ., level_good) and
! e%outermost_contour%%n_points /= 0 and no good contour at level_bad}
end do
level_good = innermost_level
level_try = merge(maxval(ssh), minval(ssh), e%cyclone)
call assert((level_try - e%ssh_extremum) &
/ (level_good - e%ssh_extremum) > 1., &
"set_outermost_contour level_try")
tentative_contour = good_contour(corner, step, ssh, level_try, &
e%coord_extr, nearby_extr)
if (tentative_contour%n_points /= 0) then
e%outermost_contour%polyline = tentative_contour
level_good = level_try
else
! {no good contour at level_try}
level_bad = level_try
do while (abs(level_bad - level_good) > precision)
level_try = (level_good + level_bad) / 2.
tentative_contour = good_contour(corner, step, ssh, level_try, &
e%coord_extr, nearby_extr)
if (tentative_contour%n_points /= 0) then
level_good = level_try
e%outermost_contour%polyline = tentative_contour
else
level_bad = level_try
end if
! {e%outermost_contour%polyline == good_contour(. . .,
! level_good) and e%outermost_contour%%n_points /= 0 and
! no good contour at level_bad}
end do
end if
e%outermost_contour%ssh = level_good
e%outermost_contour%area &
= spherical_polygon_area(e%outermost_contour%polyline)
end if
e%outermost_contour%ssh = level_good
e%outermost_contour%area &
= spherical_polygon_area(e%outermost_contour%polyline)
end subroutine set_outermost_contour
end module set_outermost_contour_m
&main_nml NOISE_AROUND=t /
#!/usr/bin/env python3
import netCDF4
import matplotlib.pyplot as plt
import numpy as np
with netCDF4.Dataset("h.nc") as f:
ssh = f.variables["adt"][:].squeeze()
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.show()
......@@ -9,6 +9,7 @@ import shapefile
import numpy as np
import matplotlib.pyplot as plt
import netCDF4
import matplotlib.markers
parser = argparse.ArgumentParser()
parser.add_argument("-v", "--velocity", help = "plot velocity field",
......@@ -37,14 +38,14 @@ for shape_rec in reader.iterShapeRecords():
suff_amp_list.append(shape_rec.record[5])
points = np.array(points)
plt.plot(points[:, 0], points[:, 1], "o", markersize = 10, fillstyle = "none",
color = "black")
colors = np.choose(cyclone, ("red", "blue"))
plt.scatter(points[:, 0], points[:, 1], s = 100, edgecolors = colors,
facecolors = "none")
for i, e in enumerate(points, start = 1):
plt.annotate(str(i), e, xytext = (3, 3), textcoords = "offset points")
colors = np.choose(cyclone, ("red", "blue"))
# Outermost and max-speed contours:
reader_outer = shapefile.Reader("outermost_contour_1")
......@@ -53,20 +54,21 @@ reader_m_s = shapefile.Reader("max_speed_contour_1")
for shape_outer, shape_m_s, my_color, suff_amp in zip(reader_outer.iterShapes(),
reader_m_s.iterShapes(),
colors, suff_amp_list):
points = np.array(shape_outer.points)
lines = plt.plot(points[:, 0], points[:, 1], color = my_color)
if suff_amp == 0:
lines[0].set_marker("s")
lines[0].set_fillstyle("none")
elif shape_m_s.shapeType == shapefile.NULL:
lines[0].set_marker("x")
else:
lines[0].set_marker("o")
if shape_m_s.shapeType != shapefile.NULL:
points = np.array(shape_m_s.points)
plt.plot(points[:, 0], points[:, 1], color = my_color, marker = "o")
if shape_outer.shapeType != shapefile.NULL:
points = np.array(shape_outer.points)
lines = plt.plot(points[:, 0], points[:, 1], color = my_color)
if suff_amp == 0:
lines[0].set_marker("s")
lines[0].set_fillstyle("none")
elif shape_m_s.shapeType == shapefile.NULL:
lines[0].set_marker("x")
else:
lines[0].set_marker("o")
if shape_m_s.shapeType != shapefile.NULL:
points = np.array(shape_m_s.points)
plt.plot(points[:, 0], points[:, 1], color = my_color, marker = "o")
if args.velocity:
with netCDF4.Dataset("uv.nc") as f:
......
......@@ -3,7 +3,6 @@ program test_get_snapshot
use derived_types, only: snapshot
use dispatch_snapshot_m, only: dispatch_snapshot
use get_snapshot_m, only: get_snapshot
use local_extrema_m, only: construct_mask_center
use jumble, only: new_unit
use netcdf, only: nf90_nowrite
use netcdf95, only: nf95_open, find_coord, nf95_inquire_dimension, &
......@@ -57,7 +56,6 @@ program test_get_snapshot
call nf95_close(ncid)
call construct_mask_center
call get_snapshot(s, min_amp, m = 1, n_proc = 1, k_end = 1, max_delta = 4, &
nlon = nlon, nlat = nlat, k = 1, max_radius = [20, 20], &
corner = [longitude(1), latitude(1)] * deg_over_rad, &
......
program test_local_extrema
use jumble, only: new_unit
use local_extrema_m, only: local_extrema, construct_mask_center
use local_extrema_m, only: local_extrema
use netcdf, only: nf90_nowrite, nf90_clobber, NF90_FLOAT, NF90_INT
use netcdf95, only: nf95_open, nf95_close, nf95_inq_varid, nf95_get_var, &
nf95_put_var, nf95_create, nf95_def_dim, nf95_def_var, nf95_put_att, &
......@@ -54,7 +54,6 @@ program test_local_extrema
call nf95_close(ncid)
call construct_mask_center
call local_extrema(ssh, extr_map, ind_extr, innermost_level, cyclone)
call new_unit(unit)
......
......@@ -30,23 +30,22 @@ my_min[:,0] = False
my_min[:,-1] = False
ilat_max, ilon_max = np.where(my_max)
print("Maxima (longitude index, latitude index) (1-based):")
print(np.column_stack((ilon_max + 1, ilat_max + 1)))
print("Maxima (longitude, latitude):")
print(np.column_stack((longitude[ilon_max], latitude[ilat_max])))
ilat_min, ilon_min = np.where(my_min)
print("Minima (longitude index, latitude index) (1-based):")
print(np.column_stack((ilon_min + 1, ilat_min + 1)))
print("Minima:")
print(np.column_stack((longitude[ilon_min], latitude[ilat_min])))
print("Values around maxima:")
print("Maxima:")
print("(longitude index, latitude index) (1-based) -- (longitude, latitude):")
print("Values around maxima\n")
for j,i in np.argwhere(my_max):
print("({}, {}) -- ({}, {})".format(i + 1, j + 1, longitude[i],
latitude[j]))
print(ssh[j-1:j+2, i-1:i+2], "\n")
print("Values around minima:")
print("Minima:")
print("(longitude index, latitude index) (1-based) -- (longitude, latitude):")
print("Values around minima\n")
for j,i in np.argwhere(my_min):
print("({}, {}) -- ({}, {})".format(i + 1, j + 1, longitude[i],
latitude[j]))
print(ssh[j-1:j+2, i-1:i+2], "\n")
f.close()
......
......@@ -34,10 +34,9 @@ program test_set_outermost_contour
real, parameter:: step = 0.25 / 180. * pi ! in rad
TYPE(shpfileobject) shphandle
integer field_number, shape_number
logical:: noise_around = .false.
namelist /main_nml/ ilon_llc, ilat_llc, ilon_urc, ilat_urc, ind_targ_extr, &
innermost_level, cyclone, noise_around
innermost_level, cyclone
!----------------------------------------------------------------
......@@ -82,8 +81,7 @@ program test_set_outermost_contour
e%coord_extr = [longitude(ind_targ_extr(1)), latitude(ind_targ_extr(2))]
e%cyclone = cyclone
call set_outermost_contour(e, ind_targ_extr, innermost_level, extr_map, &
ssh, corner = [longitude(1), latitude(1)], step = [step, step], &
noise_around = noise_around)
ssh, corner = [longitude(1), latitude(1)], step = [step, step])
print *, "e%outermost_contour%closed = ", e%outermost_contour%closed
print *, "Radius of disk of equal area: ", &
......
......@@ -46,7 +46,7 @@
"$input_dir/extr_map_negative_2_8.nc"],
"stdout" : "test_set_outermost_contour_stdout.txt",
"directory" : "Tests_new/Set_outermost_contour_noise_2_8",
"stdin" : "$stdin_dir/set_outermost_contour_nml.txt"
"stdin" : "$stdin_dir/default_main_nml.txt"
},
{
"args" : ["$compil_prod_dir/test_set_outermost_contour",
......@@ -54,7 +54,7 @@
"$input_dir/extr_map_negative_2.nc"],
"stdout" : "test_set_outermost_contour_stdout.txt",
"directory" : "Tests_new/Set_outermost_contour_noise_2",
"stdin" : "$stdin_dir/set_outermost_contour_nml.txt"
"stdin" : "$stdin_dir/default_main_nml.txt"
},
{
"args" : ["$compil_prod_dir/test_max_speed_contour_ssh",
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment