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()