source: TOOLS/MOSAIX/Build_coordinates_mask.py @ 6041

Last change on this file since 6041 was 5951, checked in by snguyen, 3 years ago

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