tnc2cdo.py - pism - [fork] customized build of PISM, the parallel ice sheet model (tillflux branch) (HTM) git clone git://src.adamsgaard.dk/pism (DIR) Log (DIR) Files (DIR) Refs (DIR) LICENSE --- tnc2cdo.py (10240B) --- 1 #!/usr/bin/env python3 2 3 # @package nc2cdo 4 # \author Andy Aschwanden, University of Alaska Fairbanks, USA 5 # \brief Script makes netCDF file ready for Climate Data Operators (CDO). 6 # \details Script adds attributes and variables to a netCDF which are required for 7 # post-processing with <a href="http://www.mad.zmaw.de/Pingo/post/post.cdo.home.html">CDO</a>. 8 # That is, the attribute 'coordinates = "lon lat"' is added to any gridded variable 9 # found in the file, except for lat and lon. For remapping methods other than bilinear, 10 # CDO additionally needs the corners of each grid cell (the lat/lon variables define the 11 # cell centers). We compute the corners of each grid cell on the x-y grid because this is 12 # relatively straightforward (compared to the lat-lon world), and then project the corners 13 # onto the lat-lon grid. The script attemps to retrieve the required mapping information from, 14 # if present, a global attribute called 'projection' which (hopefully) contains a valid 15 # Proj projection string. If this string is not available (this is true for standard PISM 16 # output), the script searches variables for the 'grid_mapping' attribute and translates information 17 # from the corresponding mapping variable into a Proj projection string. If neither is found, 18 # the script exists with an error message. 19 # 20 # Usage: 21 # 22 # \verbatim $ nc2cdo.py foo.nc \endverbatim 23 24 import sys 25 import numpy as np 26 from argparse import ArgumentParser 27 28 from pyproj import Proj 29 30 # try different netCDF modules 31 try: 32 from netCDF4 import Dataset as CDF 33 except: 34 print("netCDF4 is not installed!") 35 sys.exit(1) 36 37 # Set up the option parser 38 parser = ArgumentParser() 39 parser.description = '''Script makes netCDF file ready for Climate Data Operators (CDO). Either a global attribute "projection", a mapping variable, or a command-line proj string or a EPSG code must be given.''' 40 parser.add_argument("FILE", nargs=1) 41 parser.add_argument("--no_bounds", dest="bounds", action="store_false", 42 help="do not add lat/lon bounds.", default=True) 43 parser.add_argument("--srs", dest="srs", 44 help=''' 45 a valid proj string describing describing the projection 46 ''', default=None) 47 options = parser.parse_args() 48 args = options.FILE 49 srs = options.srs 50 bounds = options.bounds 51 52 if len(args) == 1: 53 nc_outfile = args[0] 54 else: 55 print('wrong number arguments, 1 expected') 56 parser.print_help() 57 exit(0) 58 59 # Get projection information from netCDF file 60 61 62 def get_projection_from_file(nc): 63 64 # First, check if we have a global attribute 'proj' 65 # which contains a Proj string: 66 try: 67 p = Proj(str(nc.proj)) 68 print( 69 'Found projection information in global attribute proj, using it') 70 except: 71 try: 72 p = Proj(str(nc.projection)) 73 print( 74 'Found projection information in global attribute projection, using it') 75 except: 76 try: 77 # go through variables and look for 'grid_mapping' attribute 78 for var in list(nc.variables.keys()): 79 if hasattr(nc.variables[var], 'grid_mapping'): 80 mappingvarname = nc.variables[var].grid_mapping 81 print( 82 'Found projection information in variable "%s", using it' % mappingvarname) 83 break 84 var_mapping = nc.variables[mappingvarname] 85 p = Proj(proj="stere", 86 ellps=var_mapping.ellipsoid, 87 datum=var_mapping.ellipsoid, 88 units="m", 89 lat_ts=var_mapping.standard_parallel, 90 lat_0=var_mapping.latitude_of_projection_origin, 91 lon_0=var_mapping.straight_vertical_longitude_from_pole, 92 x_0=var_mapping.false_easting, 93 y_0=var_mapping.false_northing) 94 except: 95 print('No mapping information found, exiting.') 96 sys.exit(1) 97 98 return p 99 100 101 if __name__ == "__main__": 102 103 # open netCDF file in 'append' mode 104 nc = CDF(nc_outfile, 'a') 105 106 # a list of possible x-dimensions names 107 xdims = ['x', 'x1'] 108 # a list of possible y-dimensions names 109 ydims = ['y', 'y1'] 110 111 # assign x dimension 112 for dim in xdims: 113 if dim in list(nc.dimensions.keys()): 114 xdim = dim 115 # assign y dimension 116 for dim in ydims: 117 if dim in list(nc.dimensions.keys()): 118 ydim = dim 119 120 # coordinate variable in x-direction 121 x_var = nc.variables[xdim] 122 # coordinate variable in y-direction 123 y_var = nc.variables[ydim] 124 125 # values of the coordinate variable in x-direction 126 easting = x_var[:] 127 # values of the coordinate variable in y-direction 128 northing = y_var[:] 129 130 # grid spacing in x-direction 131 de = np.diff(easting)[0] 132 # grid spacing in y-direction 133 dn = np.diff(northing)[0] 134 135 # number of grid points in x-direction 136 M = easting.shape[0] 137 # number of grid points in y-direction 138 N = northing.shape[0] 139 140 # number of grid corners 141 grid_corners = 4 142 # grid corner dimension name 143 grid_corner_dim_name = "nv4" 144 145 # array holding x-component of grid corners 146 gc_easting = np.zeros((M, grid_corners)) 147 # array holding y-component of grid corners 148 gc_northing = np.zeros((N, grid_corners)) 149 # array holding the offsets from the cell centers 150 # in x-direction (counter-clockwise) 151 de_vec = np.array([-de / 2, de / 2, de / 2, -de / 2]) 152 # array holding the offsets from the cell centers 153 # in y-direction (counter-clockwise) 154 dn_vec = np.array([-dn / 2, -dn / 2, dn / 2, dn / 2]) 155 # array holding lat-component of grid corners 156 gc_lat = np.zeros((N, M, grid_corners)) 157 # array holding lon-component of grid corners 158 gc_lon = np.zeros((N, M, grid_corners)) 159 160 if srs: 161 # use projection from command line 162 try: 163 proj = Proj(init=srs) 164 except: 165 proj = Proj(srs) 166 else: 167 # Get projection from file 168 proj = get_projection_from_file(nc) 169 170 if bounds: 171 for corner in range(0, grid_corners): 172 ## grid_corners in x-direction 173 gc_easting[:, corner] = easting + de_vec[corner] 174 # grid corners in y-direction 175 gc_northing[:, corner] = northing + dn_vec[corner] 176 # meshgrid of grid corners in x-y space 177 gc_ee, gc_nn = np.meshgrid( 178 gc_easting[:, corner], gc_northing[:, corner]) 179 # project grid corners from x-y to lat-lon space 180 gc_lon[:, :, corner], gc_lat[:, :, corner] = proj( 181 gc_ee, gc_nn, inverse=True) 182 183 # If it does not yet exist, create dimension 'grid_corner_dim_name' 184 if bounds and grid_corner_dim_name not in list(nc.dimensions.keys()): 185 nc.createDimension(grid_corner_dim_name, size=grid_corners) 186 187 var = 'lon_bnds' 188 # Create variable 'lon_bnds' 189 if not var in list(nc.variables.keys()): 190 var_out = nc.createVariable( 191 var, 'f', dimensions=(ydim, xdim, grid_corner_dim_name)) 192 else: 193 var_out = nc.variables[var] 194 # Assign units to variable 'lon_bnds' 195 var_out.units = "degreesE" 196 # Assign values to variable 'lon_nds' 197 var_out[:] = gc_lon 198 199 var = 'lat_bnds' 200 # Create variable 'lat_bnds' 201 if not var in list(nc.variables.keys()): 202 var_out = nc.createVariable( 203 var, 'f', dimensions=(ydim, xdim, grid_corner_dim_name)) 204 else: 205 var_out = nc.variables[var] 206 # Assign units to variable 'lat_bnds' 207 var_out.units = "degreesN" 208 # Assign values to variable 'lat_bnds' 209 var_out[:] = gc_lat 210 211 ee, nn = np.meshgrid(easting, northing) 212 lon, lat = proj(ee, nn, inverse=True) 213 214 var = 'lon' 215 # If it does not yet exist, create variable 'lon' 216 if not var in list(nc.variables.keys()): 217 var_out = nc.createVariable(var, 'f', dimensions=(ydim, xdim)) 218 else: 219 var_out = nc.variables[var] 220 # Assign values to variable 'lon' 221 var_out[:] = lon 222 # Assign units to variable 'lon' 223 var_out.units = "degreesE" 224 # Assign long name to variable 'lon' 225 var_out.long_name = "Longitude" 226 # Assign standard name to variable 'lon' 227 var_out.standard_name = "longitude" 228 if bounds: 229 # Assign bounds to variable 'lon' 230 var_out.bounds = "lon_bnds" 231 232 var = 'lat' 233 # If it does not yet exist, create variable 'lat' 234 if not var in list(nc.variables.keys()): 235 var_out = nc.createVariable(var, 'f', dimensions=(ydim, xdim)) 236 else: 237 var_out = nc.variables[var] 238 # Assign values to variable 'lat' 239 var_out[:] = lat 240 # Assign units to variable 'lat' 241 var_out.units = "degreesN" 242 # Assign long name to variable 'lat' 243 var_out.long_name = "Latitude" 244 # Assign standard name to variable 'lat' 245 var_out.standard_name = "latitude" 246 if bounds: 247 # Assign bounds to variable 'lat' 248 var_out.bounds = "lat_bnds" 249 250 # Make sure variables have 'coordinates' attribute 251 for var in list(nc.variables.keys()): 252 if (nc.variables[var].ndim >= 2): 253 nc.variables[var].coordinates = "lon lat" 254 255 # lat/lon coordinates must not have mapping and coordinate attributes 256 # if they exist, delete them 257 for var in ['lat', 'lon', 'lat_bnds', 'lon_bnds']: 258 if hasattr(nc.variables[var], 'grid_mapping'): 259 delattr(nc.variables[var], 'grid_mapping') 260 if hasattr(nc.variables[var], 'coordinates'): 261 delattr(nc.variables[var], 'coordinates') 262 263 # If present prepend history history attribute, otherwise create it 264 from time import asctime 265 histstr = asctime() + \ 266 ' : grid info for CDO added by nc2cdo.py, a PISM utility\n' 267 if 'History' in nc.ncattrs(): 268 nc.History = histstr + nc.History 269 elif 'history' in nc.ncattrs(): 270 nc.history = histstr + nc.history 271 else: 272 nc.history = histstr 273 274 for attr in ("projection", "proj"): 275 if hasattr(nc, attr): 276 delattr(nc, attr) 277 # Write projection attribute 278 nc.proj = proj.srs 279 # Close file 280 nc.close()