| 1 | ### =========================================================================== |
|---|
| 2 | ### |
|---|
| 3 | ### Fill a coupler restart, ocean side |
|---|
| 4 | ### |
|---|
| 5 | ### =========================================================================== |
|---|
| 6 | ## |
|---|
| 7 | ## Warning, to install, configure, run, use any of Olivier Marti's |
|---|
| 8 | ## software or to read the associated documentation you'll need at least |
|---|
| 9 | ## one (1) brain in a reasonably working order. Lack of this implement |
|---|
| 10 | ## will void any warranties (either express or implied). |
|---|
| 11 | ## O. Marti assumes no responsability for errors, omissions, |
|---|
| 12 | ## data loss, or any other consequences caused directly or indirectly by |
|---|
| 13 | ## the usage of his software by incorrectly or partially configured |
|---|
| 14 | ## personal. |
|---|
| 15 | ## |
|---|
| 16 | ## SVN information |
|---|
| 17 | __Author__ = "$Author$" |
|---|
| 18 | __Date__ = "$Date$" |
|---|
| 19 | __Revision__ = "$Revision$" |
|---|
| 20 | __Id__ = "$Id$" |
|---|
| 21 | __HeadURL = "$HeadURL$" |
|---|
| 22 | ## |
|---|
| 23 | import netCDF4 |
|---|
| 24 | import numpy.ma as np |
|---|
| 25 | from scipy import ndimage |
|---|
| 26 | import nemo |
|---|
| 27 | import shutil, getopt, sys |
|---|
| 28 | |
|---|
| 29 | def MyFill (InputData=None) : |
|---|
| 30 | """ |
|---|
| 31 | Replace the value of masked 'InputData' cells by the value of the nearest valid data cell |
|---|
| 32 | From https://stackoverflow.com/questions/5551286/filling-gaps-in-a-numpy-array |
|---|
| 33 | |
|---|
| 34 | Input: InputData : numpy masked array of any dimension. |
|---|
| 35 | |
|---|
| 36 | Output: Return a filled array. |
|---|
| 37 | """ |
|---|
| 38 | |
|---|
| 39 | Invalid = np.where ( InputData[:,:].mask, True, False) |
|---|
| 40 | Indices = ndimage.distance_transform_bf(Invalid, metric='euclidean', sampling=None, return_distances=False, |
|---|
| 41 | return_indices=True, distances=None ) |
|---|
| 42 | FilledData = InputData[tuple(Indices)] |
|---|
| 43 | return FilledData |
|---|
| 44 | |
|---|
| 45 | def usage () : |
|---|
| 46 | texte = """%(prog)s usage : |
|---|
| 47 | python %(prog)s [-d] [-i <orca grid file>] [-n <perio>] |
|---|
| 48 | -d | --debug : debug |
|---|
| 49 | -i <file> | --input=<file> : input file (default: none) |
|---|
| 50 | -v <varlist> | --variable=<variables> : list of variable to fill (defautl: none) |
|---|
| 51 | -n <perio> | --perio=<perio> : periodicity type (default: try to guess) |
|---|
| 52 | # 1, 4, 6 : Cyclic on i dimension (generaly longitudes) |
|---|
| 53 | # 2 : Obsolete (was symmetric condition at southern boundary ?) |
|---|
| 54 | # 3, 4 : North fold T-point pivot (legacy ORCA2) |
|---|
| 55 | # 5, 6 : North fold F-point pivot (ORCA1, ORCA025, ORCA2 with new grid for paleo) |
|---|
| 56 | If <perio> is not specified, %(prog)s will try to guess it from the grid dimensions |
|---|
| 57 | example : |
|---|
| 58 | python %(prog)s -n 4 -i sstoce_ORCA2.3.nc |
|---|
| 59 | python %(prog)s -n 4 -i sstoce_ORCA2.3.nc -v O_SSTSST,OIceFrc |
|---|
| 60 | """ |
|---|
| 61 | print ( texte % { 'prog':sys.argv[0] } ) |
|---|
| 62 | |
|---|
| 63 | ## Default input parameters |
|---|
| 64 | InFile = None |
|---|
| 65 | OuFile = None |
|---|
| 66 | replace = False |
|---|
| 67 | ListVarName = None |
|---|
| 68 | nperio = None |
|---|
| 69 | Debug = False |
|---|
| 70 | ListExclude = None ; Exclude=False |
|---|
| 71 | |
|---|
| 72 | ## Command line options |
|---|
| 73 | try: |
|---|
| 74 | myopts, myargs = getopt.getopt ( sys.argv[1:], 'i:o:rv:xp:hd', [ 'input=', 'output=', 'replace', 'exclude', 'variable=', 'perio=', 'debug', 'help' ] ) |
|---|
| 75 | except getopt.GetoptError as cmdle : |
|---|
| 76 | print ( "Command line error : "+str(cmdle)+"\n" ) |
|---|
| 77 | usage () |
|---|
| 78 | sys.exit(1) |
|---|
| 79 | |
|---|
| 80 | for myopt, myval in myopts : |
|---|
| 81 | if myopt in [ '-h', '--help' ] : usage () ; sys.exit (0) ; |
|---|
| 82 | elif myopt in [ '-d', '--debug' ] : Debug = True ; |
|---|
| 83 | elif myopt in [ '-i', '--input' ] : InFile = myval ; |
|---|
| 84 | elif myopt in [ '-o', '--output' ] : |
|---|
| 85 | if replace : |
|---|
| 86 | print ("Error : you can not specify both -r|--replace and -o|--output" ) |
|---|
| 87 | usage () |
|---|
| 88 | sys.exit(1) |
|---|
| 89 | else : |
|---|
| 90 | if Debug : print ("Out file set to "+OuFile) |
|---|
| 91 | OuFile = myval ; |
|---|
| 92 | elif myopt in [ '-r', '--replace' ] : |
|---|
| 93 | if OuFile != None : |
|---|
| 94 | print ("Error : you can not specify both -r|--replace and -o|--output" ) |
|---|
| 95 | usage () |
|---|
| 96 | sys.exit(1) |
|---|
| 97 | else : |
|---|
| 98 | if Debug : print ("Out file set to input file") |
|---|
| 99 | OuFile = InFile ; s |
|---|
| 100 | elif myopt in [ '-p', '--perio' ] : nperio = int(myval) ; |
|---|
| 101 | elif myopt in [ '-v', '--variable' ] : |
|---|
| 102 | if Exclude : |
|---|
| 103 | ListExclude = myval.split(',') |
|---|
| 104 | else : |
|---|
| 105 | ListVarName = myval.split(',') |
|---|
| 106 | elif myopt in [ '-x', '--exclude' ] : |
|---|
| 107 | Exclude = True |
|---|
| 108 | if ListVarName != None : |
|---|
| 109 | ListExclude = ListVarName |
|---|
| 110 | ListVarName = None |
|---|
| 111 | |
|---|
| 112 | if OuFile == None : |
|---|
| 113 | OuFile = InFile.replace ( ".nc", "_filled.nc" ) |
|---|
| 114 | print ("Creates output file name: "+OuFile) |
|---|
| 115 | |
|---|
| 116 | # Copy the input file if needed |
|---|
| 117 | if OuFile != InFile : shutil.copyfile ( InFile, OuFile ) |
|---|
| 118 | |
|---|
| 119 | print ("Output file: " + OuFile ) |
|---|
| 120 | |
|---|
| 121 | # Open file |
|---|
| 122 | OuFile = netCDF4.Dataset ( OuFile , "r+" ) |
|---|
| 123 | |
|---|
| 124 | # Try to guess periodicity type |
|---|
| 125 | if nperio == None : |
|---|
| 126 | print ("Trying to guess nperio parameter") |
|---|
| 127 | jpoi = OuFile.dimensions["x"].size |
|---|
| 128 | jpoj = OuFile.dimensions["y"].size |
|---|
| 129 | print ("Grid dimensions: ("+str(jpoj)+", "+str(jpoi)+")") |
|---|
| 130 | if (jpoj, jpoi) == (149, 182) : |
|---|
| 131 | print ("ORCA 2 grid found: nperio may vary for this configuration") |
|---|
| 132 | print ("Choosen nperio=4") |
|---|
| 133 | nperio = 4 |
|---|
| 134 | if (jpoj, jpoi) == (332, 292) : |
|---|
| 135 | nperio = 6 |
|---|
| 136 | print ("ORCA1 grid found, nperio=6" ) |
|---|
| 137 | if (jpoj, jpoi) == (332, 362) : |
|---|
| 138 | nperio = 6 |
|---|
| 139 | print ("eORCA1 grid found, nperio=6" ) |
|---|
| 140 | |
|---|
| 141 | if nperio == None : |
|---|
| 142 | print ("%(prog)s could not guess the periodicity type of your file") |
|---|
| 143 | print ("Please specify -n|--nperio") |
|---|
| 144 | usage () |
|---|
| 145 | sys.exit(1) |
|---|
| 146 | |
|---|
| 147 | # Get variables from file is need |
|---|
| 148 | if ListVarName == None : ListVarName = list ( OuFile.variables.keys() ) |
|---|
| 149 | |
|---|
| 150 | # Exclude some var if needed |
|---|
| 151 | if Exclude : |
|---|
| 152 | for Var in ListExclude : |
|---|
| 153 | if Var in ListVarName : ListVarName.remove(Var) |
|---|
| 154 | |
|---|
| 155 | # Loop on variables |
|---|
| 156 | for VarName in ListVarName : |
|---|
| 157 | |
|---|
| 158 | Var = OuFile.variables[VarName] |
|---|
| 159 | if 'mask' in dir(Var[:,...]) : |
|---|
| 160 | print ( "Working on " + VarName ) |
|---|
| 161 | NewVar = MyFill ( InputData = Var[:,:] ) |
|---|
| 162 | NewVar = nemo.lbc (NewVar, nperio=nperio, cd_type='T', psgn=1.0) |
|---|
| 163 | OuFile.variables[VarName][:,:] = NewVar[:,:] |
|---|
| 164 | else : |
|---|
| 165 | print ( VarName + " is not masked" ) |
|---|
| 166 | |
|---|
| 167 | # Close file : writes update variables. |
|---|
| 168 | #OuFile.close() |
|---|
| 169 | |
|---|
| 170 | ## =========================================================================== |
|---|
| 171 | ## |
|---|
| 172 | ## That's all folk's !!! |
|---|
| 173 | ## |
|---|
| 174 | ## =========================================================================== |
|---|
