Lesson 12 Using NetCDF Data

The netCDF of Uinidata is a popular binary format used by the scientific community to archive dense grid data because of its simple interface and multi-platform compatibility. A large numbers of meteorological and oceanic data in netCDF format can be downloaded from many servers of research institutes.

The utility library of zeGraph includes an object for reading and creating netCDF data. The object supports all variable types in reading and writing and all attribute types in reading, but only character and double types of attribute in writing. To read a unknown netCDF file, the first step is to find out its contents, i.e., dimensions, variable types, and attributes using the info() function of zeNetCDF. For example, assuming you have the file etopo120.cdf (search google and download if necessary), these lines of code.

require("register")
nc = zeUtl.new("cdf")
nc:open("etopo120.cdf")
nc:info()

prints out these:

Dimensions:
    ETOPO120X = 180
    ETOPO120Y = 90

Variables:
    doule ETOPO120X(ETOPO120X)
        units: degrees_east
        modulo:  
        point_spacing: even

    doule ETOPO120Y(ETOPO120Y)
        units: degrees_north
        point_spacing: even

    float ROSE(ETOPO120Y, ETOPO120X)
        missing_value:  -1.000000e+034
        _FillValue:  -1.000000e+034
        long_name: RELIEF OF THE SURFACE OF THE EARTH
        history: From etopo120
        units: METERS

Global attributes:
    history: FERRET V4.45 (GUI) 22-May-97

If you are familiar with netCDF specifications, you immediately know what are in the file. It is common for a netCDF file to include the attribute of valid range for all variables. This file does not. First, let's see data of the ETOPO120Y variable.

require("register")
nc, data = zeUtl.new("cdf", "double")
nc:open("etopo120.cdf")
nc:getvar("ETOPO120Y", data)
data:print()

--[[Outputs:
-8.900000e+001
-8.700000e+001
...
8.700000e+001
8.900000e+001

Well, they are just latitude values. Now, let's have a look at the longitude variable.

require("register")
nc, data = zeUtl.new("cdf", "double")
nc:open("etopo120.cdf")
nc:getvar("ETOPO120X", data)
data:print()

--[[Outputs:
2. 100000e+001
2.300000e+001
2.500000e+001
...
3.750000e+002
3.770000e+002
3.790000e+002
]]

You may be surprised that the longitude ranges from 21 degree to 379 degree, because such a arrangement makes using the data a bit more difficult. What we intend to show you here is how to use the earth relief data to create a global map of filld continents. Because of the strange longitudinal range, rearranging the relief data is necessary. Fortunately, you can do this with just one function call of zeArray object. If you want a better quality map, get relief data with higher resolutions from NOAA's server.

require("register")

nc, rose, data, xyz, x, y, z, idx
    = zeUtl.new("cdf", "float", "double",
      "double", "double", "double", "double", "char")

nc:open("x:\\geo-data\\etopo120.cdf")
nc:getvar("ROSE", rose)

-- Convert data from float to double

ptr, type = rose:getptr()
nele, nvec = rose:size()
data:setptr(ptr, type, nele * nvec)

data:reshape(180, 90)
data:transpose()

-- This lines rearrange the data so that the they
-- correspond to longitude of 1 to 359.

data:shift(-10) 

-- Scale the data for coloring later

max = data:max()
data:div(max)

-- Convert grid data to x, y and z

zeMake.grid2points(xyz, 1, data)
xyz:getarr(0, x)
xyz:getarr(1, y)
xyz:getarr(2, z)

-- find continent mask
      
z:ge(0, idx)

n = idx:count()
xyz:resize(n, 3)
x:getidx(idx, data)
xyz:setarr(0, data)

y:getidx(idx, data)
xyz:setarr(1, data)

z:getidx(idx, data)
xyz:setarr(2, data)

render, scene, node, point, vert
    = zeGrf.new("render", "scene", "node", "point", "vertex")

render:add(scene)
scene:set{node = node}
node:add(point)
point:set{vertex = vert, color = {0, 1, 0, 1}}
vert:add(xyz)

render:tofile("earth_relief.png")

earth_relief.png

How about creating your own netCDF file. Very simple. We suggest that you use the dumping tool of netCDF to verify the contents of the following example.

require("register")

nc, dat = zeUtl.new("cdf", "double")
nc:open("my.cdf", "w")

nc:defdim("lat", 181)
nc:defdim("lon", 361)

nc:defvar("lat", "double", 0)
nc:defvar("lon", "double", 1)
nc:defvar("sst", "double", 0, 1)

dat:range(-90, 1, 181)
nc:putvar("lat", dat)

dat:range(0, 1, 361)
nc:putvar("lon", dat)

dat:range(0, 1, 361 * 181)
nc:putvar("sst", dat)

nc:putatt("", "Copoyright", "zeGraph 2004")
nc:putatt("sst", "range", 0, 361 * 181 - 1)