Skip to content

Commit 88fcdc8

Browse files
MesharouMesharou
committed
update for gigatl1
1 parent a2815c5 commit 88fcdc8

File tree

2 files changed

+242
-2
lines changed

2 files changed

+242
-2
lines changed

‎examples/datarmor/detection_gigatl1.py‎

Lines changed: 234 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,234 @@
1+
from datetime import datetime, timedelta
2+
from py_eddy_tracker import start_logger
3+
from py_eddy_tracker.dataset.grid import UnRegularGridDataset
4+
from numpy import zeros, arange, float
5+
from netCDF4 import Dataset
6+
7+
# for plotting
8+
from matplotlib import pyplot as plt
9+
10+
import sys,os,shutil
11+
12+
13+
class RomsDataset(UnRegularGridDataset):
14+
15+
__slots__ = list()
16+
17+
def init_speed_coef(self, uname, vname):
18+
# xi_u and eta_v must be specified because this dimension are not use in lon/lat
19+
print(self.indexs)
20+
u = self.grid(uname, indexs=self.indexs)
21+
v = self.grid(vname, indexs=self.indexs)
22+
print('u.shape',u.shape)
23+
24+
u = self.rho_2d(u.T).T
25+
v = self.rho_2d(v)
26+
self._speed_norm = (v ** 2 + u ** 2) ** .5
27+
28+
29+
@staticmethod
30+
def rho_2d(x):
31+
"""Transformation to have u or v on same grid than h
32+
"""
33+
M, Lp = x.shape
34+
new_x = zeros((M + 1, Lp))
35+
new_x[1:-1] = .5 * (x[:-1] + x[1:])
36+
new_x[0] = new_x[1]
37+
new_x[-1] = new_x[-2]
38+
return new_x
39+
40+
@staticmethod
41+
def psi2rho(var_psi):
42+
"""Transformation to have vrt on same grid than h
43+
"""
44+
M, L = var_psi.shape
45+
Mp=M+1; Lp=L+1
46+
Mm=M-1; Lm=L-1
47+
var_rho = zeros((Mp, Lp))
48+
var_rho[1:M,1:L]=0.25*(var_psi[0:Mm,0:Lm]+var_psi[0:Mm,1:L]+var_psi[1:M,0:Lm]+var_psi[1:M,1:L])
49+
var_rho[0,:]=var_rho[1,:]
50+
var_rho[Mp-1,:]=var_rho[M-1,:]
51+
var_rho[:,0]=var_rho[:,1]
52+
var_rho[:,Lp-1]=var_rho[:,L-1]
53+
return var_rho
54+
55+
##########################
56+
57+
varname = sys.argv[1]
58+
59+
print("varname id %s", varname)
60+
61+
#varname = 'zeta'
62+
63+
simul = 'gigatl6'
64+
65+
##########################
66+
67+
if __name__ == '__main__':
68+
start_logger().setLevel('DEBUG')
69+
70+
71+
72+
####
73+
# create case-specific folder
74+
folder = varname + '/' + sys.argv[2] + '/'
75+
76+
try:
77+
os.mkdir(folder)
78+
except OSError:
79+
print ("Directory %s already exists" % folder)
80+
81+
82+
# copy script in folder
83+
shutil.copy('detection_gigatl6.py',folder)
84+
85+
86+
87+
####
88+
# define simulation
89+
90+
if 'gigatl1' in simul:
91+
folder_simulation = '~/megatl/GIGATL1/GIGATL1_1h/SURF/final/'
92+
elif 'gigatl6' in simul:
93+
folder_simulation = './HIS/'
94+
95+
# Pick a depth/isopycnal
96+
97+
if varname == 'zeta':
98+
s_rho = -1
99+
depth = 0
100+
elif varname == 'ow':
101+
#isopycnal
102+
filename = './iso/gigatl6_1h_isopycnal_section.01440.nc'
103+
nc = Dataset(filename,'r')
104+
isopycnals = nc.variables['isopycnal'][:]
105+
nc.close()
106+
107+
s_rho = 1 #
108+
depth = isopycnals[s_rho]
109+
110+
111+
112+
# Time loop
113+
#for time in range(34800, 34801):
114+
for time in range(1440, 6000):
115+
116+
117+
# hourly data
118+
#dtfile = 1
119+
# 12-hourly
120+
dtfile = 12
121+
122+
tfile = 5*24//dtfile
123+
124+
###########
125+
realyear_origin = datetime(2004,1,15)
126+
date = realyear_origin + timedelta(days=float(time)*dtfile/24.)
127+
128+
infiletime = time%tfile
129+
filetime = time - time%tfile
130+
131+
132+
# Using times:
133+
#filename = './GIGATL6/gigatl6_1h_horizontal_section.' + '{0:05}'.format(filetime) + '.nc'
134+
135+
if varname =='ow':
136+
#filename = './GIGATL6/gigatl6_1h_horizontal_section.' + '{0:05}'.format(filetime) + '.nc'
137+
filename = './iso/gigatl6_1h_isopycnal_section.' + '{0:05}'.format(filetime) + '.nc'
138+
139+
140+
elif varname == 'zeta':
141+
142+
143+
# or using dates
144+
date1 = realyear_origin + timedelta(days=float(filetime)*dtfile/24.)
145+
date2 = date1 + timedelta(days=float(4.5))
146+
147+
if 'gigatl1' in simul:
148+
filedate = '{0:04}'.format(date1.year)+'-'+\
149+
'{0:02}'.format(date1.month) + '-'+\
150+
'{0:02}'.format(date1.day)
151+
152+
filename = folder_simulation + 'gigatl1_surf.' + filedate + '.nc'
153+
gridname = '~/megatl/GIGATL1/gigatl1_grd.nc'
154+
155+
elif 'gigatl6' in simul:
156+
filedate = '{0:04}'.format(date1.year)+'-'+\
157+
'{0:02}'.format(date1.month) + '-'+\
158+
'{0:02}'.format(date1.day)+ '-'+\
159+
'{0:04}'.format(date2.year)+'-'+\
160+
'{0:02}'.format(date2.month) + '-' +\
161+
'{0:02}'.format(date2.day)
162+
163+
print(filedate)
164+
165+
filename = './HIS/GIGATL6_1h_inst_surf_' + filedate + '.nc'
166+
gridname = filename
167+
168+
# Identification
169+
if varname=='zeta' and 'gigatl6' in simul:
170+
lon_name, lat_name = 'nav_lon_rho', 'nav_lat_rho'
171+
elif varname=='ow':
172+
lon_name, lat_name = 'lon', 'lat'
173+
else:
174+
lon_name, lat_name = 'lon_rho', 'lat_rho'
175+
176+
# domain grid points: x1, x2, y1, y2
177+
x1, x2 = 400, 1200
178+
y1, y2 = 400, 1200
179+
180+
x1, x2 = 0, 1600
181+
y1, y2 = 0, 2000
182+
183+
h = RomsDataset(filename, lon_name, lat_name, gridname = gridname,
184+
indexs=dict(time=infiletime,time_counter=infiletime,
185+
eta_rho=slice(y1, y2),
186+
xi_rho=slice(x1, x2),
187+
eta_v=slice(y1, y2-1),
188+
xi_u=slice(x1, x2-1),
189+
y_rho=slice(y1, y2),
190+
x_rho=slice(x1, x2),
191+
y_v=slice(y1, y2-1),
192+
y_u=slice(y1, y2),
193+
x_u=slice(x1, x2-1),
194+
x_v=slice(x1, x2),
195+
s_rho=s_rho)
196+
)
197+
198+
199+
200+
# Identification
201+
if varname=='zeta':
202+
z_min = -2 ; z_max = 1.5; step = 0.02
203+
elif varname=='ow':
204+
# ow is multiplied by 1e10
205+
z_min = -1; z_max = -0.01; step = 0.01
206+
207+
208+
a, c = h.eddy_identification(varname, 'u', 'v', date, z_min = z_min, z_max = z_max, step = step, pixel_limit=(10, 2000), shape_error=40, force_height_unit='m',force_speed_unit='m/s',vorticity_name='vrt')
209+
210+
211+
####
212+
213+
filename = 'gigatl6_1h_' + varname + '_'+ '{0:04}'.format(depth) + '_'
214+
215+
with Dataset(date.strftime(folder + 'Anticyclonic_' + filename + '%Y%m%d%H.nc'), 'w') as h:
216+
a.to_netcdf(h)
217+
with Dataset(date.strftime(folder + 'Cyclonic_' + filename + '%Y%m%d%H.nc'), 'w') as h:
218+
c.to_netcdf(h)
219+
220+
# PLOT
221+
222+
223+
fig = plt.figure(figsize=(15,7))
224+
ax = fig.add_axes([.03,.03,.94,.94])
225+
ax.set_title('Eddies detected -- Cyclonic(red) and Anticyclonic(blue)')
226+
#ax.set_ylim(-75,75)
227+
ax.set_xlim(250,360)
228+
ax.set_aspect('equal')
229+
a.display(ax, color='b', linewidth=.5)
230+
c.display(ax, color='r', linewidth=.5)
231+
ax.grid()
232+
fig.savefig(folder + 'eddies_' + date.strftime( filename + '%Y%m%d%H') +'.png')
233+
234+

‎src/py_eddy_tracker/dataset/grid.py‎

Lines changed: 8 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -255,6 +255,7 @@ class GridDataset(object):
255255
"y_dim",
256256
"coordinates",
257257
"filename",
258+
"gridname",
258259
"dimensions",
259260
"indexs",
260261
"variables_description",
@@ -270,12 +271,13 @@ class GridDataset(object):
270271
N = 1
271272

272273
def __init__(
273-
self, filename, x_name, y_name, centered=None, indexs=None, unset=False
274+
self, filename, x_name, y_name, gridname=None, centered=None, indexs=None, unset=False
274275
):
275276
"""
276277
:param str filename: Filename to load
277278
:param str x_name: Name of longitude coordinates
278279
:param str y_name: Name of latitude coordinates
280+
:param str gridname: Name of file containing lon,lat (filename is used if not defined)
279281
:param bool,None centered: Allow to know how coordinates could be used with pixel
280282
:param dict indexs: A dictionary that sets indexes to use for non-coordinate dimensions
281283
:param bool unset: Set to True to create an empty grid object without file
@@ -292,6 +294,10 @@ def __init__(
292294
self.centered = centered
293295
self.contours = None
294296
self.filename = filename
297+
if gridname is not None:
298+
self.gridname = gridname
299+
else:
300+
self.gridname = filename
295301
self.coordinates = x_name, y_name
296302
self.vars = dict()
297303
self.indexs = dict() if indexs is None else indexs
@@ -384,7 +390,7 @@ def load(self):
384390
Get coordinates and setup coordinates function
385391
"""
386392
x_name, y_name = self.coordinates
387-
with Dataset(self.filename) as h:
393+
with Dataset(self.gridname) as h:
388394

389395
self.x_dim = h.variables[x_name].dimensions
390396
self.y_dim = h.variables[y_name].dimensions

0 commit comments

Comments
 (0)