source: TOOLS/MOSAIX/Build_coordinates_mask.py @ 6173

Last change on this file since 6173 was 6124, checked in by snguyen, 2 years ago

solved bug for some bathymetry files by removing try section and always do except section when renaming y and x at beginning.

  • Property svn:keywords set to Date Author Revision Id SVN_Date
File size: 18.8 KB
Line 
1### ===========================================================================
2###
3### Build coordnates mask
4###
5### ===========================================================================
6# creates file coordinates_mask.nc used by mosaix from NEMO netcdf files
7# coordinates.nc and bathymetry.nc (on the same grid as coordinates.nc)
8
9import numpy  as np
10import xarray as xr
11import nemo
12import datetime, os, platform, argparse
13
14## SVN information
15__Author__   = "$Author$"
16__Date__     = "$Date$"
17__Revision__ = "$Revision$"
18__Id__       = "$Id$"
19__HeadURL__  = "$HeadURL: $"
20###
21
22### ===== Handling command line parameters ==================================================
23# Creating a parser
24parser = argparse.ArgumentParser (
25    description = """Read coordinates and bathymetry to build masks, grid bounds and areas for MOSAIX""",
26    epilog='-------- End of the help message --------',
27    formatter_class=argparse.ArgumentDefaultsHelpFormatter)
28
29# Adding arguments
30parser.add_argument ('--model'   , help='oce model name', type=str, default='eORCA1.2', choices=('paleORCA', 'ORCA2.3', 'ORCA2.4', 'eORCA1.2', 'eORCA025', 'eORCA025.1') )
31parser.add_argument ('--bathy'   , help='bathymetry file name', type=str, default='bathy_meter.nc' )
32parser.add_argument ('--coord'   , help='coordinates file name', type=str, default='coordinates.nc' )
33parser.add_argument ('--fout'    , help='output file name (given as an input to MOSAIX)', type=str, default='coordinates_mask.nc' )
34parser.add_argument ('--ocePerio', help='periodicity of ocean grid, if not given is imposed according to --model', type=int, default=0, choices=range(1, 7) )
35parser.add_argument ('--nemo4Ubug',help='reproduce NEMO lbc_lnk bug for U grid and periodicity 4', type=bool, default=False, choices=(True, False) )
36parser.add_argument ('--straits' , help='modify grid metrics for selected strait points (depends on model name)', type=bool, default=False, choices=(True, False) )
37parser.add_argument ('--coordPerio', help='impose periodicity of coordinates and areas through lbc', type=bool, default=False, choices=(True, False) )
38parser.add_argument ('--maskbathy', help='use the same formula as in NEMO to compute mask_T from the bathymetry file. Does not always give the same result as NEMO ???', type=bool, default=False, choices=(True, False) )
39
40# Parse command line
41myargs = parser.parse_args()
42
43# ocean model name
44model = myargs.model
45
46# bathymetry input file
47n_bathy = myargs.bathy
48# coordinates input file
49n_coord = myargs.coord
50# coordinates and mask output file
51n_out = myargs.fout
52
53# reproduce periodicity lbclnk bug for U and nperio = 4
54nemoUbug = myargs.nemo4Ubug
55# change metrics (and bathymetry) for straits
56straits = myargs.straits
57# periodicity imposed on coordinates, metrics and areas
58coordperio = myargs.coordPerio
59# function used to compute mask_T from bathymetry
60maskbathynemo = myargs.maskbathy
61
62# type of grid periodicity
63nperio = myargs.ocePerio
64
65# check of periodicity with type of grid
66if model in ('eORCA1.2', 'paleORCA', 'ORCA025', 'eORCA025', 'eORCA025.1') :
67    if nperio != 0 :
68        if nperio != 6 :
69            print(f'Warning ! model = {model} and ocePerio = {nperio} instead of 6 !')
70    else :
71        nperio = 6
72
73if model in ('ORCA2.3', 'ORCA2.4') :
74    if nperio != 0 :
75        if nperio != 4 :
76            print(f'Attention ! model = {model} and ocePerio = {nperio} instead of 4 !')
77    else :
78        nperio = 4
79
80print(f' model = {model}\n ocePerio = {nperio}\n bathy = {n_bathy}\n coord = {n_coord}\n fout  = {n_out}')
81print(f' Switchs : nemo4Ubug = {nemoUbug}, straits = {straits}, coordPerio = {coordperio}, maskbathy = {maskbathynemo}')
82
83##
84#!!! bathymetry and coordinates files
85##
86
87# open input files while removing the time dimension
88f_coord = xr.open_dataset (n_coord, decode_times=False).squeeze()
89f_bathy = xr.open_dataset (n_bathy, decode_times=False).squeeze()
90
91# Suppress time if necessary
92try    :
93    del f_coord['time']
94    print ('time successfully removed')
95except :
96    pass
97    print ('failed to suppress time')
98
99# rename latitude, longitude and grid variables in bathymetry
100Bathymetry = f_bathy['Bathymetry'].copy()
101
102nav_lon_grid_T = f_bathy['nav_lon'].data
103nav_lat_grid_T = f_bathy['nav_lat'].data
104Bathymetry = xr.DataArray (Bathymetry, coords = { "nav_lat_grid_T": (["y", "x"], nav_lat_grid_T),
105                                                  "nav_lon_grid_T": (["y", "x"], nav_lon_grid_T) } )
106   
107Bathymetry = Bathymetry.rename ({'y':'y_grid_T', 'x':'x_grid_T'})
108
109# modify Bathymetry in straits for select grids (might change computed mask_T)
110if straits :
111    # Open straits for usual grids
112    if model == 'ORCA2.3' :
113        # orca_r2: Gibraltar strait open
114        Bathymetry[101,139] = 284.
115        # orca_r2: Bab el Mandeb strait open
116        Bathymetry[87,159] = 137.
117
118# impose periodicity of bathymetry
119Bathymetry = nemo.lbc (Bathymetry, nperio=nperio, cd_type='T')
120# suppress ocean points at southernmost position only if nperio in {3, 4, 5, 6}
121if nperio in (3, 4, 5, 6) :
122    Bathymetry[0,:] = 0.0
123
124##
125#!!! Create masks from bathymetry
126##
127
128# Creation of mask_T following choosed option maskbathy
129if maskbathynemo :
130    # Use same formula as domzgr.
131    mask_T = xr.where (Bathymetry - 1. + 0.1 >= 0.0, 1, 0).astype (dtype='f4')
132else :
133    mask_T = xr.where (Bathymetry > 0.0, 1, 0).astype (dtype='f4')
134
135# Creation of U, V, W, F masks from mask_T
136mask_U = mask_T * mask_T.shift (x_grid_T=-1)
137mask_V = mask_T * mask_T.shift (y_grid_T=-1)
138mask_F = mask_T * mask_T.shift (y_grid_T=-1) * mask_T.shift (x_grid_T=-1) * mask_T.shift (y_grid_T=-1, x_grid_T=-1)
139mask_W = mask_T
140
141# loop on TUVFW to modify coordinates and attributes of mask_[TUVFW] and maskutil_[TUVFW]
142for cd_type in ['T', 'U', 'V', 'F', 'W'] :
143    MaskName = 'mask_'     + cd_type
144    UtilName = 'maskutil_' + cd_type
145    # impose periodicity of chosen grid model on masks
146    locals()[MaskName] = nemo.lbc (locals()[MaskName], nperio=nperio, cd_type=cd_type, nemo_4U_bug=nemoUbug).astype (dtype='f4')
147    # rename masks coordinates
148    if cd_type != 'T' :
149        locals()[MaskName] = locals()[MaskName].rename \
150                    ( {'y_grid_T'      :       'y_grid_'+cd_type, 'x_grid_T'      :       'x_grid_'+cd_type,
151                       'nav_lat_grid_T': 'nav_lat_grid_'+cd_type, 'nav_lon_grid_T': 'nav_lon_grid_'+cd_type} )
152
153    # create masks without duplicate points : maskutil_[TUVWF]
154    locals()[UtilName] = nemo.lbc_mask (locals()[MaskName].copy(), nperio=nperio, cd_type=cd_type)
155 
156    #set name attribute of mask dataset
157    locals()[MaskName].name = MaskName
158    locals()[UtilName].name = UtilName
159
160    # remove _FillVallue from mskutil
161    locals()[MaskName].encoding['_FillValue'] = None
162    locals()[UtilName].encoding['_FillValue'] = None
163
164    # add masks attributes
165    locals()[MaskName].attrs['cell_measures'] = 'area: area_grid_'+cd_type
166    locals()[UtilName].attrs['cell_measures'] = 'area: area_grid_'+cd_type
167
168##
169#!!! create grid coordinates from NEMO coordinates.nc file
170##
171
172angle = { 'lon' : 'glam', 'lat' : 'gphi' }
173gridv = { 'T' : 't', 'U' : 'u', 'V' : 'v', 'F' : 'f', 'W' : 't' }
174for cd_type in ['T', 'U', 'V', 'F', 'W'] :
175    for dir_type in ['lon', 'lat'] :
176        coord_name = 'nav_' + dir_type + '_grid_' + cd_type
177        dir_name = angle[dir_type]+gridv[cd_type]
178# impose or not periodicity on coordinates read from file
179        if coordperio :
180            locals()[coord_name] = nemo.lbc (f_coord[dir_name].copy(), nperio=nperio, cd_type=cd_type, nemo_4U_bug=nemoUbug)
181        else :
182            locals()[coord_name] = f_coord[dir_name].copy()
183
184        locals()[coord_name] = locals()[coord_name].rename( {'y':'y_grid_'+cd_type, 'x':'x_grid_'+cd_type} )
185
186# remove _FillValue and missing_value
187        locals()[coord_name].encoding['_FillValue'] = None
188        locals()[coord_name].encoding['missing_value'] = None
189
190# define name attribute
191        locals()[coord_name].name = coord_name
192
193# add coordinates attributes
194        locals()[coord_name].attrs['bounds'] ='bounds_' + dir_type + '_grid_' + cd_type
195
196    locals()['nav_lon_grid_'+cd_type].attrs['standard_name'] = 'longitude'
197    locals()['nav_lon_grid_'+cd_type].attrs['long_name']     = 'Longitude'
198    locals()['nav_lon_grid_'+cd_type].attrs['units']         = 'degrees_east'
199    locals()['nav_lat_grid_'+cd_type].attrs['standard_name'] = 'latitude'
200    locals()['nav_lat_grid_'+cd_type].attrs['long_name']     = 'Latitude'
201    locals()['nav_lat_grid_'+cd_type].attrs['units']         = 'degrees_north'
202
203##
204#!!! compute areas of cells at coordinate points
205##
206
207# create areas variables for grid cells from NEMO coordinates.nc file
208# create new variables e1 e2 to keep f_coord the same
209for cd_type in ['t', 'u', 'v', 'f'] :
210    for axis in ['1', '2'] :
211        coordName = 'e' + axis + cd_type
212        locals()[coordName]=f_coord[coordName].copy()
213        # remove zero values from areas
214        # need to be define for the extended grid south of -80S
215        # some point are undefined but you need to have e1 and e2 .NE. 0
216        locals()[coordName]=xr.where(locals()[coordName] == 0.0, 1.0e2, locals()[coordName])
217
218# Correct areas for straits
219if straits :
220    # ORCA R2 configuration
221    if model == 'ORCA2.3' :
222        # Gibraltar     : e2u reduced to 20 km
223        e2u[101,138:140] = 20.e3
224        # Bab el Mandeb : e2u reduced to 30 km
225        #                e1v reduced to 18 km
226        e1v[87,159] = 18.e3
227        e2u[87,159] = 30.e3
228        # Danish Straits: e2u reduced to 10 km
229        e2u[115,144:146] = 10.e3
230    # ORCA R1 configuration
231    if model == 'eORCA1.2' :
232        # Gibraltar : e2u reduced to 20 km
233        e2u[240,281:283] = 20.e3
234        # Bhosporus : e2u reduced to 10 km
235        e2u[247,313:315] =  10.e3
236        # Lombok : e1v reduced to 13 km
237        e1v[163:165,43] =  13.e3
238        # Sumba : e1v reduced to 8 km
239        e1v[163:165,47] =  8.e3
240        # Ombai : e1v reduced to 13 km
241        e1v[163:165,52] = 13.e3
242        # Timor Passage : e1v reduced to 20 km
243        #e1v[163:165,55] = 20.e3
244        # W Halmahera : e1v reduced to 30 km
245        e1v[180:182,54] = 30.e3
246        # E Halmahera : e1v reduced to 50 km
247        e1v[180:182,57] = 50.e3
248    # ORCA R05 configuration
249    if model == 'ORCA.05' :
250        # Reduced e2u at the Gibraltar Strait
251        e2u[326,562:564] =  20.e3
252        # Reduced e2u at the Bosphore Strait
253        e2u[342,626:628] =  10.e3
254        # Reduced e2u at the Sumba Strait
255        e2u[231,92:94] =  40.e3
256        # Reduced e2u at the Ombai Strait
257        e2u[231,102] =  15.e3
258        # Reduced e2u at the Palk Strait
259        e2u[269,14] =  10.e3
260        # Reduced e1v at the Lombok Strait
261        e1v[231:233,86] =  10.e3
262        # Reduced e1v at the Bab el Mandeb
263        e1v[275,661] =  25.e3
264
265# compute cells areas
266
267for cd_type in ['T', 'U', 'V', 'F', 'W'] :
268    areaName = 'area_grid_' + cd_type
269    if coordperio :
270        locals()[areaName] = nemo.lbc (locals()['e1'+gridv[cd_type]]*locals()['e2'+gridv[cd_type]],
271                                                           nperio=nperio, cd_type=cd_type, nemo_4U_bug=nemoUbug)
272    else :
273        locals()[areaName] = locals()['e1'+gridv[cd_type]]*locals()['e2'+gridv[cd_type]]
274
275    # rename indices
276    locals()[areaName] = locals()[areaName].rename ({'y':'y_grid_'+cd_type, 'x':'x_grid_'+cd_type})
277
278    # add attributes
279    locals()[areaName].name = areaName
280    locals()[areaName].attrs['standard_name'] = 'cell_area'
281    locals()[areaName].attrs['units']         = 'm2'
282
283    # remove fill values
284    locals()[areaName].encoding['_FillValue'] = None
285
286##
287#!!! compute grid bounds !!!
288##
289
290#---------------------------------------------------------------
291# function to generate grid bounds from NEMO coordinates.nc file
292def set_bounds (cdgrd) :
293    '''
294    Constructs lon/lat bounds
295    Bounds are numerated counter clockwise, from bottom left
296    See NEMO file OPA_SRC/IOM/iom.F90, ROUTINE set_grid_bounds, for more details
297    '''
298    # Define offset of coordinate representing bottom-left corner
299    if cdgrd in ['T', 'W'] : 
300        icnr = -1 ; jcnr = -1
301        corner_lon = f_coord['glamf'].copy() ; corner_lat = f_coord['gphif'].copy()
302        center_lon = f_coord['glamt'].copy() ; center_lat = f_coord['gphit'].copy()
303    if cdgrd == 'U'        : 
304        icnr =  0 ; jcnr = -1
305        corner_lon = f_coord['glamv'].copy() ; corner_lat = f_coord['gphiv'].copy()
306        center_lon = f_coord['glamu'].copy() ; center_lat = f_coord['gphiu'].copy()
307    if cdgrd == 'V'        : 
308        icnr = -1 ; jcnr =  0
309        corner_lon = f_coord['glamu'].copy() ; corner_lat = f_coord['gphiu'].copy()
310        center_lon = f_coord['glamv'].copy() ; center_lat = f_coord['gphiv'].copy()
311    if cdgrd == 'F'        : 
312        icnr = -1 ; jcnr = -1
313        corner_lon = f_coord['glamt'].copy() ; corner_lat = f_coord['gphit'].copy()
314        center_lon = f_coord['glamf'].copy() ; center_lat = f_coord['gphif'].copy()
315     
316    jpj, jpi = corner_lon.shape ;
317    nvertex = 4
318    dims = ['y_grid_' + cdgrd, 'x_grid_' + cdgrd, 'nvertex_grid_' + cdgrd]
319   
320    bounds_lon = xr.DataArray (np.zeros ((jpj, jpi, nvertex)), dims=dims)
321    bounds_lat = xr.DataArray (np.zeros ((jpj, jpi, nvertex)), dims=dims)
322   
323    idx = [(jcnr,icnr), (jcnr,icnr+1), (jcnr+1,icnr+1), (jcnr+1,icnr)]
324   
325    # Compute cell vertices that can be defined,
326    # and complete with periodicity
327    for nn in range (nvertex) :
328        tmp = np.roll (corner_lon, shift=tuple(-1*np.array(idx[nn])), axis=(-2,-1))
329        bounds_lon[1:jpj,1:jpi,nn] = tmp[1:jpj,1:jpi]
330        tmp = np.roll (corner_lat, shift=tuple(-1*np.array(idx[nn])), axis=(-2,-1))
331        bounds_lat[1:jpj,1:jpi,nn] = tmp[1:jpj,1:jpi]
332        bounds_lon[:,:,nn] = nemo.lbc (bounds_lon[:,:,nn], nperio=nperio, cd_type=cdgrd, nemo_4U_bug=nemoUbug)
333        bounds_lat[:,:,nn] = nemo.lbc (bounds_lat[:,:,nn], nperio=nperio, cd_type=cdgrd, nemo_4U_bug=nemoUbug)
334   
335    # Zero-size cells at closed boundaries if cell points provided,
336    # otherwise they are closed cells with unrealistic bounds
337    if not (nperio == 1 or nperio == 4 or nperio == 6) :
338        for nn in range (nvertex) : 
339            bounds_lon[:,0,nn]   = center_lon[:,0]   # (West or jpni = 1), closed E-W
340            bounds_lat[:,0,nn]   = center_lat[:,0]
341            bounds_lon[:,jpi,nn] = center_lon[:,jpi] # (East or jpni = 1), closed E-W
342            bounds_lat[:,jpi,nn] = center_lat[:,jpi]
343    if nperio != 2 :
344        for nn in range (nvertex) :
345            bounds_lon[0,:,nn] = center_lon[0,:] # (South or jpnj = 1), not symmetric
346            bounds_lat[0,:,nn] = center_lat[0,:]
347    if nperio < 3 :
348        for nn in range (nvertex) :
349            bounds_lon[jpj,:,nn] = center_lon[jpj,:] # (North or jpnj = 1), no north fold
350            bounds_lat[jpj,:,nn] = center_lat[jpj,:]
351   
352    # Rotate cells at the north fold
353    if nperio >= 3 :
354        # Working array for location of northfold
355        z_fld = nemo.lbc (np.ones ((jpj, jpi)), nperio=nperio, cd_type=cdgrd, psgn=-1., nemo_4U_bug=nemoUbug)
356        z_fld = np.repeat((z_fld == -1.0)[...,np.newaxis],4,axis=2)
357        # circular shift of 2 indices in bounds third index
358        bounds_lon_tmp = np.roll (bounds_lon, shift=-2, axis=2)
359        bounds_lat_tmp = np.roll (bounds_lat, shift=-2, axis=2)
360        bounds_lon[:,:,:]  = np.where (z_fld, bounds_lon_tmp[:,:,:] , bounds_lon[:,:,:] )
361        bounds_lat[:,:,:]  = np.where (z_fld, bounds_lat_tmp[:,:,:] , bounds_lat[:,:,:] )
362   
363    # Invert cells at the symmetric equator
364    if nperio == 2 :
365        bounds_lon_tmp = np.roll (bounds_lon, shift=-2, axis=2)
366        bounds_lat_tmp = np.roll (bounds_lat, shift=-2, axis=2)
367        bounds_lon[0,:,:] = bounds_lon[0,:,:]
368        bounds_lat[0,:,:] = bounds_lat[0,:,:]
369   
370    #bounds_lon.attrs['coordinates'] = 'nav_lat_grid_' + cdgrd + ' nav_lon_grid_' + cdgrd
371    #bounds_lat.attrs['coordinates'] = 'nav_lat_grid_' + cdgrd + ' nav_lon_grid_' + cdgrd
372    bounds_lon.attrs['units']       = 'degrees_east'
373    bounds_lat.attrs['units']       = 'degrees_north'
374    bounds_lon.name = 'bounds_lon_grid_' + cdgrd
375    bounds_lat.name = 'bounds_lat_grid_' + cdgrd
376    # remove _FillValue
377    bounds_lon.encoding['_FillValue'] = None
378    bounds_lat.encoding['_FillValue'] = None
379
380    return bounds_lon, bounds_lat
381#------------------------------------------------------------
382
383bounds_lon_grid_T, bounds_lat_grid_T = set_bounds ('T')
384bounds_lon_grid_U, bounds_lat_grid_U = set_bounds ('U')
385bounds_lon_grid_V, bounds_lat_grid_V = set_bounds ('V')
386bounds_lon_grid_W, bounds_lat_grid_W = set_bounds ('W')
387bounds_lon_grid_F, bounds_lat_grid_F = set_bounds ('F')
388
389# build xarray dataset to be saved
390ds = xr.Dataset ({
391    'mask_T'      : mask_T     ,
392    'mask_U'      : mask_U     ,
393    'mask_V'      : mask_V     ,
394    'mask_W'      : mask_W     ,
395    'mask_F'      : mask_F     ,
396    'area_grid_T' : area_grid_T,
397    'area_grid_U' : area_grid_U,
398    'area_grid_V' : area_grid_V,
399    'area_grid_W' : area_grid_W,
400    'area_grid_F' : area_grid_F,
401    'maskutil_T'  : maskutil_T ,
402    'maskutil_U'  : maskutil_U ,
403    'maskutil_V'  : maskutil_V ,
404    'maskutil_W'  : maskutil_W ,
405    'maskutil_F'  : maskutil_F ,
406    'bounds_lon_grid_T': bounds_lon_grid_T,
407    'bounds_lat_grid_T': bounds_lat_grid_T,
408    'bounds_lon_grid_U': bounds_lon_grid_U,
409    'bounds_lat_grid_U': bounds_lat_grid_U,
410    'bounds_lon_grid_V': bounds_lon_grid_V,
411    'bounds_lat_grid_V': bounds_lat_grid_V,
412    'bounds_lon_grid_W': bounds_lon_grid_W,
413    'bounds_lat_grid_W': bounds_lat_grid_W,
414    'bounds_lon_grid_F': bounds_lon_grid_F,
415    'bounds_lat_grid_F': bounds_lat_grid_F,
416})
417
418#replace nav_lon nav_lat with variables obtained from NEMO coordinates.nc file
419#by construction nav_lon nav_lat come from the bathymetry
420for cd_type in ['T', 'U', 'V', 'F', 'W'] :
421    for dir_type in ['lon', 'lat'] :
422        coord_name = 'nav_' + dir_type + '_grid_' + cd_type
423        ds.coords[coord_name]=locals()[coord_name]
424
425ds.attrs['name']         = 'coordinates_mask'
426ds.attrs['description']  = 'coordinates and mask for MOSAIX'
427ds.attrs['title']        = 'coordinates_mask'
428ds.attrs['source']       = 'IPSL Earth system model'
429ds.attrs['group']        = 'ICMC IPSL Climate Modelling Center'
430ds.attrs['Institution']  = 'IPSL https.//www.ipsl.fr'
431ds.attrs['Model']        = model
432ds.attrs['timeStamp']    = '{:%Y-%b-%d %H:%M:%S}'.format (datetime.datetime.now ())
433ds.attrs['history']      = 'Build from ' + n_coord + ' and ' + n_bathy
434ds.attrs['directory']    = os.getcwd     ()
435try:
436    ds.attrs['user']         = os.getlogin   ()
437except:
438    ds.attrs['user']         = 'NoUser'
439ds.attrs['HOSTNAME']     = platform.node ()
440ds.attrs['Python']       = 'Python version: ' +  platform.python_version ()
441ds.attrs['xarray']       = 'xarray version: ' +  xr.__version__
442ds.attrs['OS']           = platform.system  ()
443ds.attrs['release']      = platform.release ()
444ds.attrs['hardware']     = platform.machine ()
445
446ds.attrs['SVN_Author']   = "$Author$"
447ds.attrs['SVN_Date']     = "$Date$"
448ds.attrs['SVN_Revision'] = "$Revision$"
449ds.attrs['SVN_Id']       = "$Id$"
450ds.attrs['SVN_HeadURL']  = "$HeadURL: http://forge.ipsl.jussieu.fr/igcmg/svn/TOOLS/MOSAIX/Build_coordinates_mask.py $"
451
452# save to output file
453ds.to_netcdf (n_out)
Note: See TracBrowser for help on using the repository browser.