• Discoverability Visible
  • Join Policy Restricted
  • Created 01 Jul 2022

Regridding with CDO

Version 4
by (unknown)
Version 5
by (unknown)

Deletions or items before changed

Additions or items after changed

1 ==Introduction to Regridding==
2
3 In order to '''intercompare different models''', they must be interpolated from the native grid on which they were calculated to a common output grid. This procedure is known as '''regridding or remapping'''.
4
5 Fundamentally, '''regridding consists of two steps''':
6
7 # Creation of an interpolation matrix (weights file), which gives the amount that each point on the source grid contributes to each point on the target grid. These weights will depend on the type of interpolation desired. ISMIP6 has standardized on First-order Conservative Remapping (Jones, P.W. 1999, Monthly Weather Review, 127, 2204-2210), which does a better job preserving fluxes (or other integrals) than, eg., bilinear interpolation.
8 # For each point on the target grid, sum the data values on the source grid multiplied by the corresponding weights. The same weights file may be used for more than one input variable, provided all are defined on the same source grid.
9 Different software may do one or the other or both of these steps.
10
11 To create the weights file, you must have ''grid description files'' for your source and target grids. These files give the latitudes and longitudes of every grid point, along with the corners of the boundary of its surrounding cell. There are several ways to write these files, depending on the software you're using.
12
13 We have tested a number of software packages, and so far, the easiest way we have found to do regridding for '''structured''' (Heiko Goelzer) and '''unstructured triangular''' (NASA Goddard) grids is using the [https://code.mpimet.mpg.de/projects/cdo Climate Data Operators (CDO)]
14
15 [[Image(ice_sheet_image_thin.png)]]
16
17
18 ==Installing CDO==
19 '''CDO''' (Climate Data Operator) claims to run on any POSIX-compatible Unix-type OS, including Mac OSX, most Linux, IBM AIX, HP-UX, Solaris, BSD variants, and even Cygwin. For those systems for which they're provided, by far the easiest way to install CDO is to use a package manager, since these will also automatically install any other libraries CDO needs (its dependencies).
20 Available packages include [https://www.macports.org/ MacPorts] for OSX, .deb for Debian and Ubuntu, and .rpm for !RedHat and Fedora. We have successfully installed CDO on OSX with !MacPorts, but have not tested others (reports welcome!).
21 Links to these packages are on the [https://code.mpimet.mpg.de/projects/cdo/wiki#Installation-and-Supported-Platforms CDO Wiki] here. To install via a package manager, you will probably need administrator privileges (this is certainly the case with !MacPorts).
22 For other systems, CDO must be built from source code. Instructions for doing so are also at the CDO web site [https://code.mpimet.mpg.de/projects/cdo/wiki#Installation-and-Supported-Platforms here]. We haven't tested this.
23
24 At the very least, you will also need the [https://www.unidata.ucar.edu/software/netcdf/ NetCDF library] from Unidata so that CDO can read and write (classic) !NetCDF files. If your model files are in !NetCDF4 or HDF5 you'll also need the HDF5 library. The package managers should install these automatically if you don't already have them.
25
26 ----
27 ==CDO Grid Description Files==
28 The critical part of remapping with CDO is to create horizontal-grid description files describing both your input and output grids. These files provide the latitude and longitude of every grid point, and of the corners of its surrounding cell boundary.
29
30 A '''CDO grid file''' for the ISMIP6 standard 5 km polar stereographic grid for the Greenland ice sheet and the 8km polar stereographic grid for the Antarctic ice sheet can be obtained by emailing ismip6@gmail.com.
31
32 Details of the file format are in section 1.3 (especially 1.3.2.4) and '''Appendix D''' of the [https://code.mpimet.mpg.de/projects/cdo/wiki#Documentation CDO User's Guide].
33
34 Here is a simple '''example''' of a grid that has '''6 cells, 3 across and 2 up''':
35 {{{
36 # Simple 3x2 CDO grid file
37
38 gridtype = curvilinear
39 gridsize = 6
40 xsize = 3
41 ysize = 2
42
43 # Longitudes
44 xvals = 301.0 303.0 305.0
45 301.0 303.0 305.0
46
47 # Longitudes of cell corners
48 xbounds =
49 302.0 302.0 300.0 300.0
50 304.0 304.0 302.0 302.0
51 306.0 306.0 304.0 304.0
52 302.0 302.0 300.0 300.0
53 304.0 304.0 302.0 302.0
54 306.0 306.0 304.0 304.0
55
56 # Latitudes
57 yvals = 61.0 61.0 61.0
58 64.0 64.0 64.0
59
60 # Latitudes of cell corners
61 ybounds =
62 60.0 63.0 63.0 60.0
63 60.0 63.0 63.0 60.0
64 60.0 63.0 63.0 60.0
65 63.0 65.0 65.0 63.0
66 63.0 65.0 65.0 63.0
67 63.0 65.0 65.0 63.0
68 }}}
69
70 '''where:'''
71
72 `gridtype = curvilinear`
73
74 Several grid types are available. '''Curvilinear''' grids are 2-dimensional, but are not necessarily arranged along latitude and longitude lines. Because the standard polar-stereographic grid is rectangular in p-s x-y space but not lon-lat space, it is curvilinear. `gridtype = unstructured` is a simple list of grid points, with cells which may or may not be the same size, may be in any order, and may not necessarily be quadrilateral, and need not be adjoining. A Voroni triagulation is represented as an unstructured grid, as is an adaptive mesh. Notice that a curvilinear grid could as easily be described as unstructured, but the output variables would be 1D instead of 2D. For other options see the CDO User's Guide.
75
76 `gridsize = 6`
77 Total number of points in the grid. For input files, this MUST match the number of points in the variables. For curvilinear grids, gridsize = xsize * ysize. For unstructured grids, this is just the number of elements in your mesh.
78
79 `xsize = 3`
80 For curvilinear grids, the number of grid points in the longitude direction. Not used for unstructured grids.
81
82 `ysize = 2`
83 For curvilinear grids, the number of grid points in the latitude direction. Not used for unstructured grids.
84
85 `nvertex = 4`
86 Number of corners of the boundary of the cell surrounding each grid point. Optional for curvilinear grids (assumes 4), required for unstructured grids. A triangulation grid would have nvertex = 3.
87
88 `xvals = ...`
89 The longitude position of every grid cell.
90
91 `yvals = ...`
92 The latitude position of every grid cell.
93
94 `xbounds = ...`
95 The longitudes of the corners of the cell boundary surrounding each grid point. There must be (nvertex * gridsize) values. CDO requires that the corners be listed counterclockwise around the grid point.
96
97 `ybounds = ...`
98 The latitudes of the corners of the cell boundary surrounding each grid point. There must be (nvertex * gridsize) values.
99
100 Grid files may contain comment lines, which begin with "#", and blank lines.
101 In the example above, the grid point at (61.,302.) is the center of a cell with corners at (60.,302.), (63.,302.), (63.,300.), and (60.,300.).
102
103 -
----
+
104 ===Grid File Notes===
105 * Each discrete '''grid point''' (xval,yval) can be considered representative of some region immediately surrounding it. This region constitutes the cell. (For example, the value at the grid point could be the mean of the value across the entire cell.) In regular grids the grid point is often at the center of the '''cell''', although it doesn't need to be. In order to conserve flux when remapping from one grid to another, it is necessary to know the area of each cell, so the cells' boundaries must be described. This is done by specifying the position of each vertex along each cell boundary (xbounds,ybounds). CDO requires that cell-boundary vertices be given counterclockwise around the grid point.
106 * CDO identifies "x" with longitude and "y" with latitude (although there are options to rotate the pole). In curvilinear grids, cells are arranged with their longitudes (x) varying fastest, then latitudes (y). This matches the CF conventions. In '''unstructured grids''', cells may be in any order.
107 * '''Latitudes and longitudes''' may have any number of decimal places. They may even be integers (assuming that is appropriate for your grid).
108 * In the grid files' lists of latitudes and longitudes (xvals, yvals, xbounds, ybounds), values are separated by white space or newlines (not by commas). There can be any number of values on one line, however the total size of a line may not exceed 65,536 characters. ('''Error''' will be "readline Warning: end of line not found (maxlen = 65536)! Abort") This can easily happen if you create the file using a word processor in which newlines are paragraph separators.
109
110 In the grid files, the keywords must be at the beginnings of their own lines. Something like this is '''NOT''' permitted:
111 {{{
112 # WRONG!
113 xvals = 300 301 302 yvals = 60 61 62
114 }}}
115 * The grid files can be large! The standard Greenland 5km grid is 20 MB, and a Greenland 1 km grid took 490 MB.
116 * CDO is also supposed to accept SCRIP format grid files, but we haven't tested this.
117
118 [[Image(ice_sheet_image_thin_2.png)]]
119
120 ----
121 ==Running the Remapping (Conservative Method)==
122 ISMIP6 has standardized on First-order Conservative Remapping (Jones, P.W. 1999, Monthly Weather Review, 127, 2204-2210), which does a better job of preserving integrals of the data, like flux, between the source and target grids, although it may have a larger local interpolation error than other methods. The CDO operator [https://code.mpimet.mpg.de/projects/cdo/embedded/index.html#x1-5710002.12.6 remapcon] implements this method, and both generates the interpolation weights matrix and applies the weights:
123 {{{
124 cdo remapcon,outgrid.txt -setgrid,ingrid.txt infile.nc outfile.nc
125 }}}
126 '''where:'''
127
128 `infile.nc`
129 is the input netCDF file
130
131 `outfile.nc`
132 is the output netCDF file
133
134 `-setgrid,ingrid.txt`
135 (hyphen required) specifies the input grid description file. Note that this overrides any grid written in the input file. The total number of points `(gridsize=)` '''MUST''' match the number of points in the variables, but the structure (points in the x and y directions) can be different. (Formally, this runs the `cdo setgrid` operator to apply the grid description in ingrid.txt to infile.nc, before running the remapcon operator.)
136
137 `remapcon,outgrid.txt`
138 says to do the conservative first-order remapping onto the grid described by outgrid.txt.
139
140 CDO's syntax requires that there can be '''NO''' spaces following the commas after remapcon or setgrid.
141
142 `remapcon` is described in [https://code.mpimet.mpg.de/projects/cdo/embedded/index.html#x1-5370002.12 '''chapter 2.12''' of the CDO User's Guide]. There is further information about conservative remapping including an outline of the mathematical basis in the SCRIP User's Guide.
143
144 By default all variables in the input netCDF file will be remapped. To only do some of them, add the cdo selname operator:
145 {{{
146 cdo remapcon,outgrid.txt -selname,var1,var2 -setgrid,ingrid.txt infile.nc outfile.nc
147 }}}
148 which will select only the variables named var1 and var2 from the input file. You can specify as many as you like. See '''section 2.3.2''' of the CDO User's Guide for more options. As with `-setgrid`, there can be no spaces following the commas.
149
150
151 [[Image(ice_sheet_image_thin.png)]]
152
153 ----
154 ==Practical Examples of CDO and NCO Tools==
155 CDO is a suite of tools that enables the regridding of datasets as mentioned above. However, there is also a suite of tools called '''NCO''' (netCDF Operators). These tools can help prepare models for submissions by '''enabling the user to add/remove/change variables, attributes''' etc.
156
157 The easiest way to install the NCO tools is via '''homebrew, macports, or anaconda''''. If that is not possible, the tools can be built from source via git. For more detailed information, consult the user guide [https://nco.sourceforge.net/nco.html#Availability here].
158
159 Below are some practical examples of these tools that may be of use to the community:
160
161 ===Do a conservative regridding of a grid from one resolution to another===
162 Note that this is using the '''remapYcon''' for the regridding.
163 {{{
164 # Names of the grid description files. In this case the 10 & 5km versions.
165 ingrid="~/initMIP/grid_description_files/CDO/ISM_10km_CDO_file.txt"
166 outgrid="~initMIP/grid_description_files/CDO/ISM_05km_CDO_file.txt"
167
168 #names of input and output model datafiles
169 infile="dlithkdt_GIS_BGC_BISICLES3_init_10km.nc"
170 outfile="dlithkdt_GIS_BGC_BISICLES3_init_05km.nc"
171
172 #Execute the regridding to 5km resolution
173 cdo remapycon,$outgrid -setgrid,$ingrid $infile $outfile
174 }}}
175 ===Copy variables from one !NetCDF file to another===
176 This command copies the variables, "y, x, mapping" from: `mapping_B03_05.nc` to
177 `dlithkdt_GIS_BGC_BISICLES3_init_05km.nc`
178
179 {{{
180 xysrc="/mapping_B03_05.nc"
181 outfile="dlithkdt_GIS_BGC_BISICLES3_init_05km.nc"
182
183 ncks -A -v y,x,mapping $xysrc $outfile
184 }}}
185 ===Add an attribute to a variable in the !NetCDF===
186 This command operates on the file, file.nc. To this file, it will add an attribute called, `grid_mapping` to the variable, `topg`. The value of the attribute, ''grid_mapping'' will be `mapping`.
187 `-h` means to not add this operation to the history of all operations performed on this file.
188 `-a` specifies that this is an attribute
189 `o` specifies that it is to be overwritten, so if there is already an attribute named, `grid_mapping`, it will be overwritten with the value, `mapping`.
190
191 `c` Specifies to create the attribute, `grid_mapping` if it does not exist already.
192 {{{
193 ncatted -h -a grid_mapping,topg,o,c,"mapping" file.nc
194 }}}
195
196
197 ===Remove unwanted Variables===
198 Removes the variables, "lat_bnds,lon_bnds,lat,lon" and outputs a new file called, outfile.nc.
199 `-C` if output file, outfile.nc, doesn't exist then it will be created
200 `-O` output to a new file.
201 `-x` extract all variables from infile.nc and output to outfile.nc, '''EXCEPT''' for the variables specified
202 `-v` controls which variables are being removed
203 {{{
204 ncks -C -O -x -v lat_bnds,lon_bnds,lat,lon infile.nc outfile.nc
205 }}}
206
207
208 ===Remove unwanted Attributes===
209 This removes the attribute, `history` from the variable, `global`.
210 This operation is performed on the file, file.nc
211 {{{
212 ncatted -h -a history,global,d,, file.nc
213 }}}
214
215
216 ===Rename a variable===
217 This command will rename the variable, `lat` to `latitude`. It is operating on the file, file.nc
218 {{{
219 ncrename -v lat,latitude file.nc
220 }}}
221
222
223 ===Convert time from years to solar seconds===
224 Assuming the file, infile.nc, has a time variable, and it is an array of years (e.g. 0,1,2,3,4,...). This will multiply the seconds per year by the year number.
225 {{{
226 #1 year interval in seconds according to Heiko's grids - close to "solar" time
227 timeincr="31556926.0"
228
229 timefactor='time=time*'$timeincr';'
230 ncap2 -O -s $timefactor infile.nc outfile.nc
231 }}}
232
233
234 ===Convert Longitudes from E/W convention to 0-360 convention===
235 This reads in the file, infile.nc, and if it has a variable named, `lon` it converts the negative longitude to a 0-360 notation.
236 {{{
237 ncap2 -O -s 'where(lon<0) lon=360+lon' infile.nc outfile.nc
238 }}}
239
240
241 ===Crop data to a certain geographic region===
242 This command will crop data to just an area around Greenland. ncea still works, but it is being replaced by the command, `nces`. Note that the file must have variables, `latitude` and `longitude`.
243 {{{
244 ncea -d latitude,59.0,84.0 -d longitude,-95,-10 infile.nc outfile.nc
245 }}}
246
247
248 ===Regridding all variables in a model submission via a bash loop===
249 This is an example of how to loop through all the variables of a group's model and do the regridding.
250 {{{
251 #!/bin/bash
252
253 #Names of the grid description files. In this case the 10 & 5km versions
254 ingrid="~/initMIP/grid_description_files/CDO/ISM_10km_CDO_file.txt"
255 outgrid="~initMIP/grid_description_files/CDO/ISM_05km_CDO_file.txt"
256
257 #grid for getting x,y, mapping variables for the "Bamber" 5km Greenland grid.
258 xysrc="/~/initMIP/QGIS/mapping_B03_05.nc"
259
260 #temporary files that get removed
261 tmp1=$inpath"/tmp1.nc"
262
263 #in this example, looking at the group, MIROC, and their ICIES00 experiment.
264 ais=GIS
265 agrp=MIROC
266 amod=ICIES00
267 #1 year interval in seconds - close to "solar" time
268 timeincr="31556926.0"
269 #----------------------------------------------------------------------
270
271 #get a list of all the variables. Assumes variables are the same for
272 #asmb, init, and ctrl directories. you don't want the "scalar"
273 #variable.
274 vars=`ls $inpath/asmb|cut -d "_" -f 1`
275 delete=(scalar)
276 vars=${vars[@]/$delete}
277
278
279 for amod in $amod;do
280 for exp in asmb init ctrl; do
281
282 #loop through each of the variables..
283 for avar in $vars; do
284 infile=${inpath}/${exp}/${avar}_${ais}_${agrp}_${amod}_${exp}.nc
285 outfile=${outpath}/${exp}/${avar}_${ais}_${agrp}_${amod}_${exp}.nc
286
287 echo $infile
288 echo $outfile
289
290 #First, the basic remapping operation. This is the command
291 #that does the actual regridding. Note that it is writing it
292 #to the temp file for now.
293 cdo remapycon,$outgrid -setgrid,$ingrid $infile $tmp1
294
295 #Now add in X,Y, and mapping variables from a file
296 #that has a known good x,y grid, and the global Grid
297 #attribute. This operation copies the variables, "y, x, and
298 #mapping" from Heiko's file to the temp file. This helps the
299 #file be read in programs like QGIS, Panoply, etc.
300 ncks -A -v y,x,mapping $xysrc $tmp1
301
302 #for QGIS to recognise the mapping variable we just added,
303 #you need to add and attribute called, "grid_mapping" to
304 #each data variable($avar)
305 ncatted -h -a grid_mapping,$avar,o,c,"mapping" $tmp1
306
307 #this converts all the variables to floats. This helps save
308 #A LOT of space, so is worth doing if double precision is not
309 #really needed
310 var2float=$avar'=float('$avar')'
311
312 #-O does an overwrite
313 ncap2 -O -s $var2float $tmp1 $tmp1
314
315 #this removes unwanted variables. Note writing out to
316 #outfile name now
317 ncks -C -O -x -v lat_bnds,lon_bnds,lat,lon $tmp1 $outfile
318
319 #make adjustments to the mask variable's attributes....
320 #space before/after brackets is essential!
321
322 if [ ${avar} == "sfrgrf" ] || [ ${avar} == "sftflf" ] || [ ${avar} == "sftgif" ]; then
323 ncatted -h -a units,$avar,o,c,"unitless" $outfile
324 fi
325
326 #add global attribute telling about interpolation
327 ncatted -h -a Remapping,global,o,c,$avar" was conservatively regridded to a 5km grid using the 'cdo remapycon' utility. More details on the remapycon utility can be found here: https://code.zmaw.de/projects/cdo or by typing 'cdo -h remapycon' at the command prompt" $outfile
328
329 #remove the history and other global attributes you may not want:
330 ncatted -h -a history,global,d,, $outfile
331 ncatted -h -a CDI,global,d,, $outfile
332 ncatted -h -a NCO,global,d,, $outfile
333 ncatted -h -a CDO,global,d,, $outfile
334 ncatted -h -a nco_openmp_thread_number,global,d,, $outfile
335 ncatted -h -a history_of_appended_files,global,d,, $outfile
336
337 rm -f $tmp1
338
339 done
340 done
341 done
342 }}}
343
344 ----
345 ==Notes==
346 * It sounds from the CDO documentation as if it ought to be possible to read the input grid information directly from the input file, but we haven't gotten this to work yet. The '''error''' is "cdo remapcon (Abort): Unsupported grid type: generic". CDO may need particular !NetCDF attributes to be in the file for this to work.
347 * CDO `remapcon` seems to handle multiple time steps just fine. However, the time dimension needs to be named `time`.
348 * The same horizontal grid must apply to all timesteps. There appears to be no way to have different horizontal grids at different timesteps. If your native grid changes during a simulation, you will therefore need to create a new grid description file for each time step where the mesh has changed.
349 * If the input file already has variables called `"lon"` and `"lat"`, these will be remapped to the new grid as with all the other variables and renamed to `"lon_2"` and `"lat_2"`. The output file will have lon and lat variables as defined by the output grid-description file.