Changeset 4151
- Timestamp:
- 11/22/18 11:45:17 (6 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
TOOLS/MOSAIX/cotes_etal.py
r4148 r4151 1 1 # -*- Mode: python -*- 2 import netCDF4 3 import numpy as np 4 import nemo 5 from scipy import ndimage 6 import sys, os, platform, argparse, textwrap, time 7 8 zero = np.dtype('float64').type(0.0) 9 zone = np.dtype('float64').type(1.0) 10 epsfrac = np.dtype('float64').type(1.0E-10) 11 pi = np.pi 12 rad = pi/np.dtype('float64').type(180.0) 13 ra = np.dtype('float64').type(6400000.0) # Earth radius 14 15 def geodist (plon1, plat1, plon2, plat2) : 16 """Distance between two points (on sphere)""" 17 zs = np.sin (rad*plat1) * np.sin (rad*plat2) + np.cos (rad*plat1) * np.cos (rad*plat2) * np.cos(rad*(plon2-plon1)) 18 zs = np.maximum (-zone, np.minimum (zone, zs)) 19 geodist = np.arccos (zs) 20 return geodist 21 2 ### =========================================================================== 22 3 ### 23 ### A partir d'un fichier de poids, complete avec les rivieres et 24 ### les poids pour le run off 25 ### On met a zero les poids sur les points non cotiers, et on renormalise 4 ### Compute runoff weights. 5 ### For LMDZ only. Not suitable for DYNAMICO 26 6 ### 27 28 29 ### Lecture de la ligne de comande 7 ### =========================================================================== 8 ## 9 ## Warning, to install, configure, run, use any of Olivier Marti's 10 ## software or to read the associated documentation you'll need at least 11 ## one (1) brain in a reasonably working order. Lack of this implement 12 ## will void any warranties (either express or implied). 13 ## O. Marti assumes no responsability for errors, omissions, 14 ## data loss, or any other consequences caused directly or indirectly by 15 ## the usage of his software by incorrectly or partially configured 16 ## personal. 17 ## 30 18 # SVN information 31 19 __Author__ = "$Author$" … … 35 23 __HeadURL__ = "$HeadURL$" 36 24 __SVN_Date__ = "$SVN_Date: $" 37 38 ### 39 40 ### ===== Handling command line parameters ================================================== 25 ## 26 27 import netCDF4 28 import numpy as np 29 import nemo 30 from scipy import ndimage 31 import sys, os, platform, argparse, textwrap, time 32 33 ## Userful constants 34 zero = np.dtype('float64').type(0.0) 35 zone = np.dtype('float64').type(1.0) 36 epsfrac = np.dtype('float64').type(1.0E-10) 37 pi = np.pi 38 rad = pi/np.dtype('float64').type(180.0) 39 ra = np.dtype('float64').type(6400000.0) # Earth radius 40 41 ## Functions 42 def geodist (plon1, plat1, plon2, plat2) : 43 """Distance between two points (on sphere)""" 44 zs = np.sin (rad*plat1) * np.sin (rad*plat2) + np.cos (rad*plat1) * np.cos (rad*plat2) * np.cos(rad*(plon2-plon1)) 45 zs = np.maximum (-zone, np.minimum (zone, zs)) 46 geodist = np.arccos (zs) 47 return geodist 48 49 50 51 52 ### ===== Reading command line parameters ================================================== 41 53 # Creating a parser 42 54 parser = argparse.ArgumentParser ( 43 55 description = """Compute calving weights""", 44 epilog='-------- This is the end of the help message --------')56 epilog='-------- End of the help message --------') 45 57 46 58 # Adding arguments … … 48 60 parser.add_argument ('--atm' , help='atm model name (LMD*)', type=str, default='LMD9695' ) 49 61 parser.add_argument ('--atmCoastWidth', help='width of the coastal band in atmosphere (in grid points)', default=1 ) 50 parser.add_argument ('--oceCoastWidth', help='width of the coastal band in ocean (in grid points)' , default=1)62 parser.add_argument ('--oceCoastWidth', help='width of the coastal band in ocean (in grid points)' , default=2 ) 51 63 parser.add_argument ('--grids', help='grids file name', default='grids.nc' ) 52 64 parser.add_argument ('--areas', help='masks file name', default='areas.nc' ) … … 55 67 parser.add_argument ('--output', help='output rmp file name', default='rmp_tlmd_to_torc_runoff.nc' ) 56 68 57 58 69 # Parse command line 59 70 myargs = parser.parse_args() … … 68 79 atm_Name = myargs.atm 69 80 oce_Name = myargs.oce 70 # Largeur de la bande cotiere cote atm81 # Width of the coastal band (land points) in the atmopshere 71 82 atmCoastWidth = myargs.atmCoastWidth 72 # Largeur de la bande cotiere cote oce83 # Width of the coastal band (ocean points) in the ocean 73 84 oceCoastWidth = myargs.oceCoastWidth 74 85 … … 94 105 atm_grid_dims = atm_grid_area.shape 95 106 (atm_nvertex, atm_jpj, atm_jpi) = atm_grid_corner_lat.shape 96 if atm_jpi == 1 : atm_perio = 0 # Icosahedre 97 else : atm_perio = 1 # 107 atm_perio = 0 98 108 atm_grid_pmask = nemo.lbc_mask (atm_grid_imask, 'T', atm_perio) 99 109 atm_address = np.reshape ( np.arange(atm_jpj*atm_jpi), (atm_jpj, atm_jpi) ) … … 115 125 oce_grid_size = oce_jpj*oce_jpi 116 126 117 print ('Calculs points terre/oce/cote sur grille atmosphere' ) 118 119 # Fill closed sea with image processing library 127 ## Fill closed sea with image processing library 120 128 oce_grid_imask = nemo.lbc_mask ( 1-ndimage.binary_fill_holes (1-nemo.lbc(oce_grid_imask, nperio=oce_perio, cd_type='T')), nperio=oce_perio, cd_type='T' ) 121 129 130 ## 122 131 print ("Determination d'une bande cotiere ocean") 123 132 124 oceLand = np.where (oce_grid_imask[:] < 0.5, True, False) 125 oceOcean = np.where (oce_grid_imask[:] > 0.5, True, False) 126 127 oceCoast = np.where (ndimage.uniform_filter(oceOcean.astype(float), size=1+2*oceCoastWidth)<1.0,True,False) & oceOcean 133 oceLand = np.where (oce_grid_pmask[:] < 0.5, True, False) 134 oceOcean = np.where (oce_grid_pmask[:] > 0.5, True, False) 135 136 NNocean = 1+2*oceCoastWidth 137 oceOceanFiltered = ndimage.uniform_filter(oceOcean.astype(float), size=NNocean) 138 oceCoast = np.where (oceOceanFiltered<(1.0-0.5/(NNocean**2)),True,False) & oceOcean 128 139 oceCoast = nemo.lbc_mask (oceCoast, oce_perio, 'T') 129 140 130 print ('N ombre de pointoceLand : ' + str(oceLand.sum()) )131 print ('N ombre de pointoceOcean : ' + str(oceOcean.sum()) )132 print ('N ombre de pointoceCoast : ' + str(oceCoast.sum()) )133 134 # Faire une liste de ces points141 print ('Number of points in oceLand : ' + str(oceLand.sum()) ) 142 print ('Number of points in oceOcean : ' + str(oceOcean.sum()) ) 143 print ('Number of points in oceCoast : ' + str(oceCoast.sum()) ) 144 145 # Arrays with coastal points only 135 146 oceCoast_grid_center_lon = oce_grid_center_lon[oceCoast] 136 147 oceCoast_grid_center_lat = oce_grid_center_lat[oceCoast] … … 142 153 143 154 print ("Determination d'une bande cotiere atmosphere " ) 144 atmLand = np.where (o2aFrac[:] < epsfrac, True, False) 145 atmFrac = np.where (o2aFrac[:] < epsfrac, True, False) & np.where (o2aFrac[:] < (zone-epsfrac), True, False) 146 atmOcean = np.where (o2aFrac[:] < (zone-epsfrac), True, False) 147 148 atmCoast = np.where (ndimage.uniform_filter(atmLand.astype(float), size=1+2*atmCoastWidth)<1.0,True,False) & atmLand 155 atmLand = np.where (o2aFrac[:] < epsfrac , True, False) 156 atmLandFrac = np.where (o2aFrac[:] < zone-epsfrac , True, False) 157 atmFrac = np.where (o2aFrac[:] > epsfrac , True, False) & np.where (o2aFrac[:] < (zone-epsfrac), True, False) 158 atmOcean = np.where (o2aFrac[:] < (zone-epsfrac), True, False) 159 atmOceanFrac = np.where (o2aFrac[:] > epsfrac , True, False) 160 161 NNatm = 1+2*atmCoastWidth 162 atmLandFiltered = ndimage.uniform_filter(atmLand.astype(float), size=NNatm) 163 atmCoast = np.where (atmLandFiltered<(1.0-0.5/(NNatm**2)),True,False) & atmLandFrac 149 164 atmCoast = nemo.lbc_mask (atmCoast, 1, 'T') 150 165 151 print ('N ombre de pointatmLand : ' + str(atmLand.sum()) )152 print ('N ombre de pointatmOcean : ' + str(atmOcean.sum()) )153 print ('N ombre de pointatmCoast : ' + str(atmCoast.sum()) )154 155 # Faire une liste de ces points166 print ('Number of points in atmLand : ' + str(atmLand.sum()) ) 167 print ('Number of points in atmOcean : ' + str(atmOcean.sum()) ) 168 print ('Number of points in atmCoast : ' + str(atmCoast.sum()) ) 169 170 # Arrays with coastal points only 156 171 atmCoast_grid_center_lon = atm_grid_center_lon[atmCoast] 157 172 atmCoast_grid_center_lat = atm_grid_center_lat[atmCoast] … … 165 180 oce_address = np.empty ( shape=(0), dtype=np.int32 ) 166 181 182 ## Loop on atmosphere coastal points 167 183 for ja in np.arange(len(atmCoast_grid_pmask)) : 168 184 z_dist = geodist ( atmCoast_grid_center_lon[ja], atmCoast_grid_center_lat[ja], oceCoast_grid_center_lon, oceCoast_grid_center_lat) … … 185 201 186 202 187 print (' fini')203 print ('End of loop') 188 204 189 205 num_links = remap_matrix.shape[0] 190 206 191 # Output file207 ### Output file 192 208 runoff = myargs.output 193 209 f_runoff = netCDF4.Dataset ( runoff, 'w', format='NETCDF3_64BIT' ) … … 203 219 f_runoff.associatedFiles = grids + " " + areas + " " + masks 204 220 f_runoff.directory = os.getcwd () 205 f_runoff.description = "Generated with ???"221 f_runoff.description = "Generated with cotes_etal.py" 206 222 f_runoff.title = runoff 207 223 f_runoff.Program = "Generated by " + sys.argv[0] + " with flags " + str(sys.argv[1:]) … … 242 258 243 259 v_remap_matrix[:] = remap_matrix 244 v_atm_address [:] = atm_address 245 v_oce_address [:] = oce_address 260 v_atm_address [:] = atm_address + 1 # OASIS uses Fortran style indexing 261 v_oce_address [:] = oce_address + 1 246 262 247 263 v_atm_grid_dims = f_runoff.createVariable ( 'src_grid_dims' , 'i4', ('src_grid_rank',) ) … … 256 272 v_atm_grid_area = f_runoff.createVariable ( 'src_grid_area' , 'f8', ('src_grid_size',) ) 257 273 v_atm_grid_area.long_name="Grid area" ; v_atm_grid_area.standard_name="cell_area" ; v_atm_grid_area.units="m2" 258 v_atm_grid_imask = f_runoff.createVariable ( 'src_grid_imask' , ' f8', ('src_grid_size',) )274 v_atm_grid_imask = f_runoff.createVariable ( 'src_grid_imask' , 'i4', ('src_grid_size',) ) 259 275 v_atm_grid_imask.long_name="Land-sea mask" ; v_atm_grid_imask.units="Land:1, Ocean:0" 276 v_atm_grid_pmask = f_runoff.createVariable ( 'src_grid_pmask' , 'i4', ('src_grid_size',) ) 277 v_atm_grid_pmask.long_name="Land-sea mask (periodicity removed)" ; v_atm_grid_pmask.units="Land:1, Ocean:0" 260 278 261 279 v_atm_grid_dims [:] = atm_grid_dims … … 266 284 v_atm_grid_area [:] = atm_grid_area[:].ravel() 267 285 v_atm_grid_imask [:] = atm_grid_imask[:].ravel() 286 v_atm_grid_pmask [:] = atm_grid_pmask[:].ravel() 268 287 269 288 # -- … … 274 293 v_oce_grid_center_lon.units='degrees_east' ; v_oce_grid_center_lon.long_name='Longitude' ; v_oce_grid_center_lon.long_name='longitude' ; v_oce_grid_center_lon.bounds="dst_grid_corner_lon" 275 294 v_oce_grid_center_lat.units='degrees_north' ; v_oce_grid_center_lat.long_name='Latitude' ; v_oce_grid_center_lat.long_name='latitude' ; v_oce_grid_center_lat.bounds="dst_grid_corner_lat" 276 v_oce_grid_corner_lon = f_runoff.createVariable ( ' oce_grid_corner_lon', 'f8', ('dst_grid_size', 'dst_grid_corners') )295 v_oce_grid_corner_lon = f_runoff.createVariable ( 'dst_grid_corner_lon', 'f8', ('dst_grid_size', 'dst_grid_corners') ) 277 296 v_oce_grid_corner_lat = f_runoff.createVariable ( 'dst_grid_corner_lat', 'f8', ('dst_grid_size', 'dst_grid_corners') ) 278 297 v_oce_grid_corner_lon.units="degrees_east" 279 298 v_oce_grid_corner_lat.units="degrees_north" 280 v_oce_grid_area = f_runoff.createVariable ( 'dst_grid_area' , 'f8', ('dst_grid_size',))299 v_oce_grid_area = f_runoff.createVariable ( 'dst_grid_area' , 'f8', ('dst_grid_size',) ) 281 300 v_oce_grid_area.long_name="Grid area" ; v_oce_grid_area.standard_name="cell_area" ; v_oce_grid_area.units="m2" 282 v_oce_grid_imask = f_runoff.createVariable ( 'dst_grid_imask' , ' f8', ('dst_grid_size',) )301 v_oce_grid_imask = f_runoff.createVariable ( 'dst_grid_imask' , 'i4', ('dst_grid_size',) ) 283 302 v_oce_grid_imask.long_name="Land-sea mask" ; v_oce_grid_imask.units="Land:1, Ocean:0" 303 v_oce_grid_pmask = f_runoff.createVariable ( 'dst_grid_pmask' , 'i4', ('dst_grid_size',) ) 304 v_oce_grid_pmask.long_name="Land-sea mask (periodicity removed)" ; v_oce_grid_pmask.units="Land:1, Ocean:0" 284 305 285 306 v_oce_grid_dims [:] = oce_grid_dims … … 290 311 v_oce_grid_area [:] = oce_grid_area[:].ravel() 291 312 v_oce_grid_imask [:] = oce_grid_imask[:].ravel() 292 293 # For diags 294 v_atm_lon_addressed = f_runoff.createVariable ( 'src_lon_addressed', 'f8', ('num_links',) ) 295 v_atm_lat_addressed = f_runoff.createVariable ( 'src_lat_addressed', 'f8', ('num_links',) ) 296 v_oce_lon_addressed = f_runoff.createVariable ( 'dst_lon_addressed', 'f8', ('num_links',) ) 297 v_oce_lat_addressed = f_runoff.createVariable ( 'dst_lat_addressed', 'f8', ('num_links',) ) 313 v_oce_grid_pmask [:] = oce_grid_pmask[:].ravel() 314 315 ## For diag, not used by OASIS 316 v_atm_lon_addressed = f_runoff.createVariable ( 'src_lon_addressed' , 'f8', ('num_links',) ) 317 v_atm_lat_addressed = f_runoff.createVariable ( 'src_lat_addressed' , 'f8', ('num_links',) ) 318 v_atm_area_addressed = f_runoff.createVariable ( 'src_area_addressed' , 'f8', ('num_links',) ) 319 v_atm_imask_addressed = f_runoff.createVariable ( 'src_imask_addressed', 'i4', ('num_links',) ) 320 v_atm_pmask_addressed = f_runoff.createVariable ( 'src_pmask_addressed', 'i4', ('num_links',) ) 321 322 v_oce_lon_addressed = f_runoff.createVariable ( 'dst_lon_addressed' , 'f8', ('num_links',) ) 323 v_oce_lat_addressed = f_runoff.createVariable ( 'dst_lat_addressed' , 'f8', ('num_links',) ) 324 v_oce_area_addressed = f_runoff.createVariable ( 'dst_area_addressed' , 'f8', ('num_links',) ) 325 v_oce_imask_addressed = f_runoff.createVariable ( 'dst_imask_addressed', 'i4', ('num_links',) ) 326 v_oce_pmask_addressed = f_runoff.createVariable ( 'dst_pmask_addressed', 'i4', ('num_links',) ) 298 327 299 328 v_atm_lon_addressed.long_name="Longitude" ; v_atm_lon_addressed.standard_name="longitude" ; v_atm_lon_addressed.units="degrees_east" 300 v_atm_lat_addressed.long_name="Latitude" ; v_atm_lat_addressed.standard_name="latitude" ; v_atm_lat_addressed.units="degrees_east" 301 v_atm_lon_addressed [:] = atm_grid_center_lon.ravel()[atm_address].ravel() 302 v_atm_lat_addressed [:] = atm_grid_center_lat.ravel()[atm_address].ravel() 329 v_atm_lat_addressed.long_name="Latitude" ; v_atm_lat_addressed.standard_name="latitude" ; v_atm_lat_addressed.units="degrees_north" 330 v_atm_lon_addressed [:] = atm_grid_center_lon.ravel()[atm_address].ravel() 331 v_atm_lat_addressed [:] = atm_grid_center_lat.ravel()[atm_address].ravel() 332 v_atm_area_addressed [:] = atm_grid_area.ravel()[atm_address].ravel() 333 v_atm_imask_addressed[:] = 1-atm_grid_imask.ravel()[atm_address].ravel() 334 v_atm_pmask_addressed[:] = 1-atm_grid_pmask.ravel()[atm_address].ravel() 303 335 304 336 v_oce_lon_addressed.long_name="Longitude" ; v_oce_lon_addressed.standard_name="longitude" ; v_oce_lon_addressed.units="degrees_east" 305 v_oce_lat_addressed.long_name="Latitude" ; v_oce_lat_addressed.standard_name="latitude" ; v_oce_lat_addressed.units="degrees_east" 306 v_oce_lon_addressed [:] = oce_grid_center_lon.ravel()[oce_address].ravel() 307 v_oce_lat_addressed [:] = oce_grid_center_lat.ravel()[oce_address].ravel() 337 v_oce_lat_addressed.long_name="Latitude" ; v_oce_lat_addressed.standard_name="latitude" ; v_oce_lat_addressed.units="degrees_north" 338 v_oce_lon_addressed [:] = oce_grid_center_lon.ravel()[oce_address].ravel() 339 v_oce_lat_addressed [:] = oce_grid_center_lat.ravel()[oce_address].ravel() 340 v_oce_area_addressed [:] = oce_grid_area.ravel()[oce_address].ravel() 341 v_oce_imask_addressed[:] = 1-oce_grid_imask.ravel()[oce_address].ravel() 342 v_oce_pmask_addressed[:] = 1-oce_grid_pmask.ravel()[oce_address].ravel() 343 344 v_atmLand = f_runoff.createVariable ( 'atmLand' , 'i4', ('src_grid_size',) ) 345 v_atmLandFiltered = f_runoff.createVariable ( 'atmLandFiltered', 'f4', ('src_grid_size',) ) 346 v_atmLandFrac = f_runoff.createVariable ( 'atmLandFrac' , 'i4', ('src_grid_size',) ) 347 v_atmFrac = f_runoff.createVariable ( 'atmFrac' , 'i4', ('src_grid_size',) ) 348 v_atmOcean = f_runoff.createVariable ( 'atmOcean' , 'i4', ('src_grid_size',) ) 349 v_atmOceanFrac = f_runoff.createVariable ( 'atmOceanFrac' , 'i4', ('src_grid_size',) ) 350 v_atmCoast = f_runoff.createVariable ( 'atmCoast' , 'i4', ('src_grid_size',) ) 351 352 v_atmLand[:] = atmLand.ravel() 353 v_atmLandFrac[:] = atmLandFrac.ravel() 354 v_atmLandFiltered[:] = atmLandFiltered.ravel() 355 v_atmFrac[:] = atmFrac.ravel() 356 v_atmOcean[:] = atmOcean.ravel() 357 v_atmOceanFrac[:] = atmOceanFrac.ravel() 358 v_atmCoast[:] = atmCoast.ravel() 359 360 v_oceLand = f_runoff.createVariable ( 'oceLand' , 'i4', ('dst_grid_size',) ) 361 v_oceOcean = f_runoff.createVariable ( 'oceOcean' , 'i4', ('dst_grid_size',) ) 362 v_oceOceanFiltered = f_runoff.createVariable ( 'oceOceanFiltered', 'f4', ('dst_grid_size',) ) 363 v_oceCoast = f_runoff.createVariable ( 'oceCoast' , 'i4', ('dst_grid_size',) ) 364 365 v_oceLand[:] = oceLand.ravel() 366 v_oceOcean[:] = oceOcean.ravel() 367 v_oceOceanFiltered[:] = oceOceanFiltered.ravel() 368 v_oceCoast[:] = oceCoast.ravel() 308 369 309 370 f_runoff.close ()
Note: See TracChangeset
for help on using the changeset viewer.