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

In test_mean_speed.py and in program test_set_max_speed, it is hard to

remember the order of positional arguments so we rather use fixed
input file names and no arguments.
parent e2c549ec
No related branches found
No related tags found
No related merge requests found
&MAIN_NML IND_TARG_EXTR=19,11 /
#!/usr/bin/env python3
import netCDF4
from mpl_toolkits import basemap
from scipy import interpolate
import shapefile
import matplotlib.pyplot as plt
import numpy as np
with netCDF4.Dataset("uv_2006_01_01.nc") as f:
u = f.variables["u"][0, 214:231, 20:49]
v = f.variables["v"][0, 214:231, 20:49]
lon = 5.125 + np.arange(29) * 0.25
lat = -36.375 + np.arange(17) * 0.25
lon2d, lat2d = np.meshgrid(lon, lat)
lon_rad = lon * np.pi / 180
lat_rad = lat * np.pi / 180
lon_0 = 9.625
lat_0 = - 33.875
lon_0_rad = lon_0 * np.pi/180
lat_0_rad = lat_0 * np.pi/180
m = basemap.Basemap(projection="stere", lon_0 = lon_0, lat_0 = lat_0,
llcrnrlon = 8, llcrnrlat = - 35, urcrnrlon = 11,
urcrnrlat=-32.5)
u_rot, v_rot = m.rotate_vector(u, v, lon, lat)
values = np.stack((u, v), axis = -1)
sf = shapefile.Reader("test_get_eddy")
plt.figure()
# Wind on the contour only
ax_wind_on_contours = plt.gca()
plt.figure()
# The whole wind field
ax_whole_wind = plt.gca()
for ax in [ax_wind_on_contours, ax_whole_wind]:
m.drawparallels(range(-35,-32), labels = [1, 0, 0, 0], ax = ax)
m.drawmeridians(range(8,12), labels = [0, 0, 0, 1], ax = ax)
m.scatter(lon_0, lat_0, latlon = True, ax = ax)
for p in sf.iterShapes():
p_array = np.array(p.points)
lon_interp_rad = p_array[:, 0]
lat_interp_rad = p_array[:, 1]
lon_interp = lon_interp_rad * 180 / np.pi
lat_interp = lat_interp_rad * 180 / np.pi
xi = np.column_stack((lat_interp_rad, lon_interp_rad))
values_x = interpolate.interpn((lat_rad, lon_rad), values, xi)
ui = values_x[:, 0]
vi = values_x[:, 1]
ui_rot, vi_rot = m.rotate_vector(ui, vi, lon_interp, lat_interp)
x = np.cos(lat_0_rad) * (lon_interp_rad - lon_0_rad)
y = lat_interp_rad - lat_0_rad
v_azim = (x * vi - y * ui) / np.sqrt(x**2 + y**2)
print("mean azimuthal speed = ", np.mean(v_azim), "m s-1")
m.plot(lon_interp, lat_interp, latlon=True, ax = ax_wind_on_contours)
m.quiver(lon_interp, lat_interp, ui_rot, vi_rot, latlon = True,
ax = ax_wind_on_contours)
m.plot(lon_interp, lat_interp, latlon=True, ax = ax_whole_wind)
m.quiver(lon2d, lat2d, u_rot, v_rot, latlon = True, ax = ax_whole_wind)
plt.show()
......@@ -9,14 +9,12 @@ import numpy as np
import sys
import math
if len(sys.argv) != 3: sys.exit("Required arguments: velocity-file shapefile")
with netCDF4.Dataset(sys.argv[1]) as f:
with netCDF4.Dataset("uv.nc") as f:
u = f.variables["u"][:]
v = f.variables["v"][:]
lon = f.variables["lon"][:]
lat = f.variables["lat"][:]
lon2d, lat2d = np.meshgrid(lon, lat)
lon_rad = lon * np.pi / 180
......@@ -24,7 +22,7 @@ lat_rad = lat * np.pi / 180
values = np.stack((u, v), axis = -1)
sf = shapefile.Reader(sys.argv[2])
sf = shapefile.Reader("test_set_max_speed")
p = sf.shape(0)
p = np.array(p.points)
......@@ -40,8 +38,8 @@ values_x = interpolate.interpn((lat_rad, lon_rad), values, xi)
ui = values_x[:, 0]
vi = values_x[:, 1]
lon_0 = 9.625
lat_0 = - 33.875
lon_0 = float(input("longitude extremum = ? "))
lat_0 = float(input("latitude extremum = ? "))
lon_0_rad = lon_0 * np.pi/180
lat_0_rad = lat_0 * np.pi/180
......
......@@ -4,11 +4,10 @@ program test_set_max_speed
use, intrinsic:: ISO_c_binding
use derived_types, only: eddy
use jumble, only: get_command_arg_dyn
use netcdf, only: nf90_nowrite
use netcdf95, only: nf95_open, nf95_close, nf95_inq_varid, nf95_get_var, &
nf95_gw_var
use nr_util, only: assert, pi
use nr_util, only: pi
use set_max_speed_m, only: set_max_speed
use shapelib, only: shpt_polygon, shpfileobject, ftdouble, shpclose, &
shpobject, dbfreadattribute
......@@ -17,7 +16,6 @@ program test_set_max_speed
implicit none
character(len = :), allocatable:: adt_file, velocity_file, extr_map_file
integer ncid, varid
real, allocatable:: ssh(:, :) ! sea-surface height, in m
real, allocatable:: u(:, :), v(:, :) ! wind, in m s-1
......@@ -35,26 +33,20 @@ program test_set_max_speed
integer ifield, ishape
logical:: cyclone = .true.
real:: ssh_extremum = 0.2759
real:: coord_extr(2) = [9.625, - 33.875] * deg_over_rad
real:: coord_extr(2) = [9.625, - 33.875]
real(kind = c_double) outermost_contour_ssh
namelist /main_nml/ cyclone, ssh_extremum, coord_extr, ind_targ_extr
!----------------------------------------------------------------
call assert(COMMAND_ARGUMENT_COUNT() == 3, &
"Required arguments: ADT-file velocity-file extr-map-file")
call get_command_arg_dyn(1, adt_file)
call get_command_arg_dyn(2, velocity_file)
call get_command_arg_dyn(3, extr_map_file)
write(unit = error_unit, nml = main_nml)
write(unit = error_unit, fmt = *) "Enter namelist main_nml."
read(unit = *, nml = main_nml)
write(unit = *, nml = main_nml)
print *, "Reading from ", adt_file, "..."
call nf95_open(adt_file, nf90_nowrite, ncid)
print *, "Reading from h.nc..."
call nf95_open("h.nc", nf90_nowrite, ncid)
call nf95_inq_varid(ncid, "adt", varid)
call nf95_gw_var(ncid, varid, ssh)
......@@ -67,8 +59,8 @@ program test_set_max_speed
call nf95_close(ncid)
print *, "Reading from ", velocity_file, "..."
call nf95_open(velocity_file, nf90_nowrite, ncid)
print *, "Reading from uv.nc..."
call nf95_open("uv.nc", nf90_nowrite, ncid)
call nf95_inq_varid(ncid, "u", varid)
call nf95_gw_var(ncid, varid, u)
......@@ -78,8 +70,8 @@ program test_set_max_speed
call nf95_close(ncid)
print *, "Reading from ", extr_map_file, "..."
call nf95_open(extr_map_file, nf90_nowrite, ncid)
print *, "Reading from extr_map.nc..."
call nf95_open("extr_map.nc", nf90_nowrite, ncid)
call nf95_inq_varid(ncid, "extr_map", varid)
call nf95_gw_var(ncid, varid, extr_map)
......@@ -95,7 +87,7 @@ program test_set_max_speed
e%outermost_contour%ssh = outermost_contour_ssh
e%cyclone = cyclone
e%ssh_extremum = ssh_extremum
e%coord_extr = coord_extr
e%coord_extr = coord_extr * deg_over_rad
e%outermost_contour%n_points = psobject%nvertices
e%outermost_contour%closed = .true.
allocate(e%outermost_contour%points(2, e%outermost_contour%n_points))
......
......@@ -80,25 +80,26 @@
"input" : "&main_nml /\n"
},
{
"args": ["$compil_prod_dir/test_set_max_speed",
"$input_dir/h_outermost.nc", "$input_dir/uv_outermost.nc",
"$input_dir/extr_map_outermost.nc"],
"args": "$compil_prod_dir/test_set_max_speed",
"title": "Set_max_speed",
"required": ["$tests_old_dir/Set_outermost_contour/test_set_outermost_contour.shp",
"required": [["$input_dir/h_outermost.nc", "h.nc"],
["$input_dir/uv_outermost.nc", "uv.nc"],
["$input_dir/extr_map_outermost.nc", "extr_map.nc"],
"$tests_old_dir/Set_outermost_contour/test_set_outermost_contour.shp",
"$tests_old_dir/Set_outermost_contour/test_set_outermost_contour.shx",
"$tests_old_dir/Set_outermost_contour/test_set_outermost_contour.dbf"],
"input" : "&main_nml /\n"
},
{
"args": ["$compil_prod_dir/test_set_max_speed",
"$input_dir/h_test_region.nc",
"$input_dir/uv_test_region.nc",
"$input_dir/extr_map_negative_2_8.nc"],
"args": "$compil_prod_dir/test_set_max_speed",
"title": "Set_max_speed_noise",
"required": ["$tests_old_dir/Set_outermost_contour_noise_2_8/test_set_outermost_contour.shp",
"required": [["$input_dir/h_test_region.nc", "h.nc"],
["$input_dir/uv_test_region.nc", "uv.nc"],
["$input_dir/extr_map_negative_2_8.nc", "extr_map.nc"],
"$tests_old_dir/Set_outermost_contour_noise_2_8/test_set_outermost_contour.shp",
"$tests_old_dir/Set_outermost_contour_noise_2_8/test_set_outermost_contour.shx",
"$tests_old_dir/Set_outermost_contour_noise_2_8/test_set_outermost_contour.dbf"],
"stdin": "$stdin_dir/set_max_speed_noise_nml.txt"
"input": "&MAIN_NML IND_TARG_EXTR=19,11 /\n"
},
{
"args": "$compil_prod_dir/test_get_snapshot",
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment