Creating depth profiles: ======================== The utility 'vertical_profiler.py' is used to generate basic statistics from ROMS/HOPS datasets that can be used by the isosurfacer. The statistics (mean and standard deviation of a given NetCDF variable) are gathered on a per-depth basis. The input files must be in rectilinear form (i.e. not in sigma-coordinates), and the depth coordinates should be the same across all datasteps. To use the vertical profiler, you pass it a variable name to gather statistics on, and a list of files to scan. The variable can be either a 3D or 4D variable. In the 4D case, the first dimension must be 'time'. The datasets are scanned, and a table is printed to stdout. Example: (save a temperature profile into the file 'temp_profile.txt' home> python vertical_profiler.py temp Aug2006_6hourly_ROMSoutput.nc > temp_profile.txt home> cat temp_profile.txt # depth mean stddev # ----- ---- ------ 0 15.67 1.742 5 15.63 1.126 10 15.4 1.22 15 14.95 1.362 20 14.33 1.428 30 13.33 1.397 40 12.63 1.328 50 12.09 1.24 60 11.69 1.174 <..snip..> The units in the table are the same as the units of the variable (in this case, depth is in meters, and temperature is in degrees Celsius). Note: ----- This script is a bit of a hack - it assumes that the depth variable will be named either 'depth' or 'ocean_depth', that time is called 'time', and that if the selected variable has a missing value sentinel, that it is named 'missing_value'. Note: ----- Unlike the isosurface extractor, vertical_profiler.py uses variable names directly from the NetCDF file, rather than "canonical" names. Using the isosurfacer roms_extractor.py: ======================================== (note - to see an example of full-fledged usage, skip to the end. The following examples build up to that usage) IMPORTANT NOTE: roms_extractor uses Python's own interpreter to evaluate both the settings file, and the expressions for isosurfacing and texture coordinates. In general, it is not possible to prevent python from executing arbitrary code on the machine running the isosurfacer. Therefore untrusted users should not be allowed to edit the settings file or enter expressions. In particular, if this code is to be used as part of a web portal or similar, the settings files and input expressions should be either externally validated, or limited to a fixed set of possible inputs. This shouldn't be too frightening though; users without malicuous intent will not be entering anything even remotely likely to cause any damage. The .settings file: ------------------- This code was written to handle output from different simulations (MBARI '03 ROMS and HOPS, MBARI 06 ROMS), and so needs to be told just a few things about which format it is getting data from, such as what the variable names are, and where the area of interest is. Hence the settings file. Here is an example for the MBARI '06 ROMS data: # A map of variable names in the dataset to canonical names. # The canonical names are : time, latitude, longitude, depth, salinity, # temperature, u, v, and w # (where u,v, and w are ocean current velocity) variable_map = { "time" : "time", "lat" : "latitude", "lon" : "longitude", "depth" : "depth", "salt" : "salinity", "temp" : "temperature", "u" : "u", "v" : "v", "w" : "w" } # A list of zero or more clip planes, for culling out masked portions of the dataset # Each item is a pair consisting of a point on the plane (as depth,latitude,longitude) # followed by the plane normal clip_planes=[( (0, 36.463, 236.92), (0, 0.41011, 0.91204)), ( (0, 35.185, 237.58), (0, 0.94733, -0.32027)), ( (0, 37.503, 236.94), (0, -0.96683, 0.25542))] This file is actually a python script that is evaluated when the isosurfacer is started. By default, the isosurfacer will try to read a settings file named 'default.settings', but that can be overridden with the --settings option. Input data format: ------------------ The isosurfacer takes NetCDF files representing 3D datasets as input, so any 4D datasets must be split up before using the isosurface (in our (NCSA's) case, the data went through a temporal interpolation phase first, to provide a longer, smoother animation) The input is mapped to a rectilinear cartesian grid, with uniform spacing in the latitude and longitudinal directions. Depth is mapped to the x coordinate, latitude to the y, and longitude to the z coordinate. For computing vertical vorticity, approximations to the horizontal and vertical spacing are made by taking the central latitude and longitude of the dataset, delta(latitude) and delta(longitude), and finding the spacing in meters using an oblate spheroid model of the earth. Using the isosurfacer: ---------------------- The basic usage to get an isosurface of a variable is: (assuming default.settings exists): vtkpython roms_extractor.py where is a mathematical expression involving temperature, salinity, u,v, w, mask (where grid points with missing values are set to 0, elsewhere 1), and vorticity_z (computed at runtime from u and v) or any other variable in the file. The is a comma-separated list of isosurface values to generate. The output is sent to standard out In general, if either or is more complex than a single variable or value, it should be quoted. * Example of basic usage: Generate an isosurface of salinity from hour 10 of the simulation at isovalue==34 PSUs, and save it in file sal_34.vtk: vtkpython roms_extractor.py Aug2006_6hourly.0010.nc \ salinity 34 >sal_34.vtk * Example of basic usage with expressions and multiple isovalues: Generate isosurfaces based on temperatures in Fahrenheit (for presentation to U.S. non-scientists, perhaps) vtkpython roms_extractor.py Aug2006_6hourly.0010.nc \ 'temperature * 5.0/9.0 + 32.0' \ '48,50,52' >temp_isos.vtk * Example of usage with multi-timestep dataset: When using a 4D dataset, the timestep to used is specified with the --timestep option. It's argument is a 0-based index into the time dimension (not the actual time values). So to generate isosurfaces from the 100th timestep: vtkpython roms_extractor.py Aug2006_6hourly_ROMSoutput.nc \ --timestep 99 \ salinity 34 >sal_34.vtk * Example using a depth profile: Generate isosurfaces where the temperature is +/- 1 standard deviation from the depth adjusted norm. Assume we have already generated a temperature profile named 'temperature_profile.txt' with the 'vertical_profiler' tool described above. Using the --profile option introduces two new symbols that can be used in the expression: mean and stddev vtkpython roms_extractor.py Aug2006_6hourly_ROMSoutput.nc \ --timestep 10 \ --profile temperature_profile.txt \ -- \ '(temperature-mean)/stddev' \ '-1,1' >temp_adjusted_isos.vtk In this case, we place '--' after all the optional command line arguments because otherwise the -1 value in the list of isosurface values would be considered an command-line switch. * Example using --sentinel: One problem we encountered was that VTK does not seem to have a good way of dealing with datasets with missing values. Unfortunately, this usually causes spurious false isosurfaces at the boundary between valid and missing values. We can use clip planes to eradicate these false surfaces in the open-water boundaries of the dataset, but this does not work for the ocean floor, especially near the shore. This can be seen if the result of the previous example is viewed in paraview. For isosurfaces which are pretty much closed, such as temperature or salinity perturbation (as opposed to layered, such as raw temperature or salinity), we can alleviate the problem somewhat by replacing missing values with a well chosen sentinel value. In general, for low-valued isosurfaces, the replacement value should be large (greater than the max value), and for high-valued isosurfaces the replacement value should be small. For this, we use the --sentinel option. Improving upon the previous example, we then have: vtkpython roms_extractor.py Aug2006_6hourly_ROMSoutput.nc \ --timestep 10 \ --sentinel 9999 \ --profile temperature_profile.txt \ -- \ '(temperature-mean)/stddev' \ -1 >temp_perturb-1.vtk vtkpython roms_extractor.py Aug2006_6hourly_ROMSoutput.nc \ --timestep 10 \ --sentinel -9999 \ --profile temperature_profile.txt \ -- \ '(temperature-mean)/stddev' \ 1 >temp_perturb+1.vtk *Example of using texture coordinates: You can also add texture coordinates to the isosurface, where the value of each u or v texture coordinate is an arbitrary expression (just like the isosurface expression). These are intended to be used as indices into a color lookup table, so that the surface can be colored by a different quantity than that which was used to produce the surface. The u and v expressions are specified by --u-texture and --v-texture respectively. Texture coordinates are linearly scaled to the range [0,1]. By default, roms_extractor maps min(u) and min(v) to 0, and max(u) and max(v) to 1. You can override this using the options --u-min, --u-max, --v-min, and --v-max. Values outside the given range will be clamped to [0,1]. We didn't actually use these texture coordinates to color the final version of our visualizations, as the scene was already fairly complex. However, here is an example where we take the temperature perturbation at 2 standard deviations above normal. Then we map actual temperature into the u texture coordinate, and salinity into the v texture coordinate. roms_extractor.py \ --sentinel -10 \ --profile temperature_profile.txt \ --u-tex temperature --u-min 6 --u-max 23 \ --v-tex salinity --v-min 32.111 --v-max 34.356 \ -- Aug2006_6hourly.0010.nc \ '(temperature-mean)/stddev' 2 \ > temp_perturb_with_tex" *Example: Real usage: This is the script we use to produce temperature perturbation isosurfaces. We generated 2 sets of isosurfaces, at +2 and -2 standard deviations, using the command lines: > gen_isos.sh +2 -9999 and > gen_isos.sh -2 9999 Note that this script uses the tools 'iota' (which takes a printf formatted string and an integral range. For each integer in the range, it outputs a line consisting of the input string formatted with the integer), and 'pardo' (given a list of commands with one command per line, performs the commands in parallel). These are both handy utilities, availible at http://dart.ncsa.uiuc.edu/tools. #! /bin/bash if [ $# -ne 2 ]; then echo "Usage: HOW.sh isoval sentinelval" exit 1 fi export CMD="\ /rd/mahall/ocean/scripts/roms_extractor.py \ --sentinel $2 \ --profile ../temperature_profile.txt \ --u-tex temperature --u-min -14 --u-max 20 \ --v-tex salinity --v-min 0.0 --v-max 56 \ -- /qd/mbari06/roms/interp8/Aug2006_6hourly.%04d.nc \ '(temperature-mean)/stddev' $1\ > /rd/mahall/ocean/mbari06/roms/temp_perturb/temp_$1.%04d.vtk" iota -n -f "$CMD" 0-975 |pardo -j2