source: Swollen/swwreader/swwreader.cpp @ 103

Last change on this file since 103 was 103, checked in by darran, 20 years ago
  • xllcorner,yllcorner offset working
File size: 14.5 KB
Line 
1/*
2    SWWReader
3
4    Reader of SWW files from within OpenSceneGraph.
5
6    copyright (C) 2004-2005 Geoscience Australia
7*/
8
9#include <SWWReader.h>
10#include <string>
11#include <fstream>
12#include <netcdf.h>
13#include <osg/Notify>
14
15#include <osgDB/Registry>
16#include <osgDB/ReadFile>
17
18#include <stdlib.h>
19#include <stdio.h>
20#include <gdal_priv.h>
21
22
23// only constructor, requires netcdf file
24SWWReader::SWWReader(const std::string& filename)
25{
26    // assume failure until end of constructor
27    _valid = false;
28
29    // netcdf filename
30    _filename = new std::string(filename);
31
32    // netcdf open
33    _status.push_back( nc_open(_filename->c_str(), NC_NOWRITE, &_ncid) );
34    if (this->_statusHasError()) return;
35
36    // dimension ids
37    _status.push_back( nc_inq_dimid(_ncid, "number_of_volumes", &_nvolumesid) );
38    _status.push_back( nc_inq_dimid(_ncid, "number_of_vertices", &_nverticesid) );
39    _status.push_back( nc_inq_dimid(_ncid, "number_of_points", &_npointsid) );
40    _status.push_back( nc_inq_dimid(_ncid, "number_of_timesteps", &_ntimestepsid) );
41    if (this->_statusHasError()) return;
42
43    // dimension values
44    _status.push_back( nc_inq_dimlen(_ncid, _nvolumesid, &_nvolumes) );
45    _status.push_back( nc_inq_dimlen(_ncid, _nverticesid, &_nvertices) );
46    _status.push_back( nc_inq_dimlen(_ncid, _npointsid, &_npoints) );
47    _status.push_back( nc_inq_dimlen(_ncid, _ntimestepsid, &_ntimesteps) );
48    if (this->_statusHasError()) return;
49
50    // variable ids
51    _status.push_back( nc_inq_varid (_ncid, "x", &_xid) );
52    _status.push_back( nc_inq_varid (_ncid, "y", &_yid) );
53    _status.push_back( nc_inq_varid (_ncid, "z", &_zid) );
54    _status.push_back( nc_inq_varid (_ncid, "volumes", &_volumesid) );
55    _status.push_back( nc_inq_varid (_ncid, "time", &_timeid) );
56    _status.push_back( nc_inq_varid (_ncid, "stage", &_stageid) );
57    if (this->_statusHasError()) return;
58
59    // allocation of variable arrays, destructor responsible for cleanup
60    _px = new float[_npoints];
61    _py = new float[_npoints];
62    _pz = new float[_npoints];
63    _ptime = new float[_ntimesteps];
64    _pvolumes = new unsigned int[_nvertices * _nvolumes];
65    _pstage = new float[_npoints];
66
67    // loading variables from netcdf file
68    _status.push_back( nc_get_var_float (_ncid, _xid, _px) );  // x vertices
69    _status.push_back( nc_get_var_float (_ncid, _yid, _py) );  // y vertices
70    _status.push_back( nc_get_var_float (_ncid, _zid, _pz) );  // bedslope heights
71    _status.push_back( nc_get_var_float (_ncid, _timeid, _ptime) );  // time array
72    _status.push_back( nc_get_var_int (_ncid, _volumesid, (int *) _pvolumes) );  // triangle indices
73    if (this->_statusHasError()) return;
74
75
76    // sww file can optionally contain bedslope texture image filename
77    size_t attlen; // length of text attribute (if it exists)
78    if( nc_inq_attlen(_ncid, NC_GLOBAL, "texture", &attlen) != NC_ENOTATT )
79    {
80       char* filename;
81       if( (filename = (char*) malloc(attlen+1)) != NULL){
82          int status;
83          status = nc_get_att_text(_ncid, NC_GLOBAL, "texture", filename);
84          if( status == NC_NOERR )
85          {
86             filename[attlen] = '\0';  // ensure string is terminated, not a requirement for netcdf attributes
87             osg::notify(osg::INFO) << "[SWWReader] embedded image filename: " << filename <<  std::endl;
88             setBedslopeTexture( std::string(filename) );
89          }
90          free( filename );
91       }
92    }
93
94    // sww file can optionally contain georeference offset
95    int status1 = nc_get_att_float(_ncid, NC_GLOBAL, "xllcorner", &_xllcorner);
96    int status2 = nc_get_att_float(_ncid, NC_GLOBAL, "yllcorner", &_yllcorner);
97    if( status1 == NC_NOERR && status2 == NC_NOERR)
98    {
99       osg::notify(osg::INFO) << "[SWWReader] xllcorner: " << _xllcorner <<  std::endl;
100       osg::notify(osg::INFO) << "[SWWReader] yllcorner: " << _yllcorner <<  std::endl;
101    }
102    else
103    {
104       _xllcorner = 0.0;  // default value
105       _yllcorner = 0.0;
106    }
107
108
109    // loop index
110    size_t iv;
111
112    // vertex indices
113    unsigned int v1index, v2index, v3index;
114
115    // bedslope range and resultant scale factors (note, these don't take into
116    // account any xllcorner, yllcorner offsets - used during texture coord assignment)
117    float xmin, xmax, xrange;
118    float ymin, ymax, yrange;
119    float zmin, zmax, zrange;
120    float aspect_ratio;
121    xmin = _px[0]; 
122    xmax = _px[0];
123    ymin = _py[0]; 
124    ymax = _py[0];
125    zmin = _pz[0]; 
126    zmax = _pz[0];
127    for( iv=1; iv < _npoints; iv++ )
128    {
129       if( _px[iv] < xmin ) xmin = _px[iv];
130       if( _px[iv] > xmax ) xmax = _px[iv];
131       if( _py[iv] < ymin ) ymin = _py[iv];
132       if( _py[iv] > ymax ) ymax = _py[iv];
133       if( _pz[iv] < zmin ) zmin = _pz[iv];
134       if( _pz[iv] > zmax ) zmax = _pz[iv];
135    }
136
137    xrange = xmax - xmin;
138    yrange = ymax - ymin;
139    zrange = zmax - zmin;
140    aspect_ratio = xrange/yrange;
141    _xscale = 1.0 / xrange;
142    _yscale = 1.0 / yrange;
143   
144    // handle case of a flat bed that doesn't necessarily pass through z=0
145    _zscale = (zrange == 0.0) ? 1.0 : 1.0/zrange;
146
147    _xoffset = xmin;
148    _yoffset = ymin;
149    _zoffset = zmin;
150    _xcenter = 0.5;
151    _ycenter = 0.5;
152    _zcenter = 0.0;
153
154    if( aspect_ratio > 1.0 )
155    {
156        _scale = 1.0/xrange;
157        _ycenter /= aspect_ratio;
158    }
159    else
160    {
161        _scale = 1.0/yrange;
162        _xcenter /= aspect_ratio;
163    }
164
165
166    osg::notify(osg::INFO) << "[SWWReader] xscale: " << _xscale <<  std::endl;
167    osg::notify(osg::INFO) << "[SWWReader] yscale: " << _yscale <<  std::endl;
168    osg::notify(osg::INFO) << "[SWWReader] zmin: " << zmin <<  std::endl;
169    osg::notify(osg::INFO) << "[SWWReader] zmax: " << zmax <<  std::endl;
170    osg::notify(osg::INFO) << "[SWWReader] zscale: " << _zscale <<  std::endl;
171
172
173
174    // alpha-scaling defaults
175    _alphamin = 0.8;       
176    _alphamax = 1.0;
177    _heightmin = 0.0;
178    _heightmax = 1.0;   
179
180
181    // compute triangle connectivity, a list (indexed by vertex number)
182    // of lists (indices of triangles sharing this vertex)
183    _connectivity = std::vector<triangle_list>(_npoints);
184    for (iv=0; iv < _nvolumes; iv++)
185    {
186        v1index = _pvolumes[3*iv+0];
187        v2index = _pvolumes[3*iv+1];
188        v3index = _pvolumes[3*iv+2];
189        _connectivity.at(v1index).push_back(iv);
190        _connectivity.at(v2index).push_back(iv);
191        _connectivity.at(v3index).push_back(iv);
192    }
193
194    // bedslope vertex array, shifting and scaling vertices to unit cube
195    // centred about the origin
196    _bedslopevertices = new osg::Vec3Array;
197    for (iv=0; iv < _npoints; iv++)
198        _bedslopevertices->push_back( osg::Vec3( (_px[iv]-_xoffset)*_scale - _xcenter, 
199                                                 (_py[iv]-_yoffset)*_scale - _ycenter, 
200                                                 (_pz[iv]-_zoffset)*_scale - _zcenter) );
201
202    // bedslope index array, pvolumes array indexes into x, y and z
203    _bedslopeindices = new osg::DrawElementsUShort(osg::PrimitiveSet::TRIANGLES, _nvolumes*_nvertices);
204    for (iv=0; iv < _nvolumes*_nvertices; iv++)
205        (*_bedslopeindices)[iv] = _pvolumes[iv];
206
207
208    // calculate bedslope primitive normal and centroid arrays
209    _bedslopenormals = new osg::Vec3Array;
210    _bedslopecentroids = new osg::Vec3Array;
211    osg::Vec3 v1, v2, v3, side1, side2, nrm;
212    for (iv=0; iv < _nvolumes; iv++)
213    {
214        v1index = _pvolumes[3*iv+0];
215        v2index = _pvolumes[3*iv+1];
216        v3index = _pvolumes[3*iv+2];
217
218        v1 = _bedslopevertices->at(v1index);
219        v2 = _bedslopevertices->at(v2index);
220        v3 = _bedslopevertices->at(v3index);
221
222        side1 = v2 - v1;
223        side2 = v3 - v2;
224        nrm = side1^side2;
225        nrm.normalize();
226
227        _bedslopenormals->push_back( nrm );
228        _bedslopecentroids->push_back( (v1+v2+v3)/3.0 );
229    }
230
231    // success!
232    _valid = true;
233
234}
235
236
237
238SWWReader::~SWWReader()
239{
240    _status.push_back( nc_close(_ncid) );
241    delete[]( _px, _py, _pz, _ptime, _pvolumes, _pstage );
242}
243
244
245
246void SWWReader::setBedslopeTexture( std::string filename )
247{
248    osg::notify(osg::INFO) << "[SWWReader::setBedslopetexture] filename: " << filename <<  std::endl;
249   _bedslopetexture = new std::string(filename);
250   _bedslopegeodata.hasData = false;
251
252   // GDAL Geospatial Image Library
253   GDALDataset *pGDAL;
254   GDALAllRegister();
255   pGDAL = (GDALDataset *) GDALOpen( filename.c_str(), GA_ReadOnly );
256   if( pGDAL == NULL )
257      osg::notify(osg::INFO) << "[SWWReader::setBedslopetexture] GDAL unable to open file: " << filename <<  std::endl;
258   else
259   {
260      osg::notify(osg::INFO) << "[SWWReader::setBedslopetexture] GDAL Driver: "
261                             << pGDAL->GetDriver()->GetDescription() << "/" 
262                             << pGDAL->GetDriver()->GetMetadataItem( GDAL_DMD_LONGNAME ) 
263                             << std::endl;
264
265      _bedslopegeodata.xresolution = pGDAL->GetRasterXSize();
266      _bedslopegeodata.yresolution = pGDAL->GetRasterYSize();
267          osg::notify(osg::INFO) << "[SWWReader::setBedslopetexture] xresolution: " << 
268          _bedslopegeodata.xresolution <<  std::endl;
269          osg::notify(osg::INFO) << "[SWWReader::setBedslopetexture] yresolution: " << 
270          _bedslopegeodata.yresolution <<  std::endl;
271
272      // check if image contains geodata
273      double geodata[6];
274      if( pGDAL->GetGeoTransform( geodata ) == CE_None )
275      {
276          _bedslopegeodata.xorigin = geodata[0];
277          _bedslopegeodata.yorigin = geodata[3];
278          _bedslopegeodata.rotation = geodata[2];
279          _bedslopegeodata.xpixel = geodata[1];
280          _bedslopegeodata.ypixel = geodata[5];
281          _bedslopegeodata.hasData = true;
282
283          osg::notify(osg::INFO) << "[SWWReader::setBedslopetexture] xorigin: " << _bedslopegeodata.xorigin <<  std::endl;
284          osg::notify(osg::INFO) << "[SWWReader::setBedslopetexture] yorigin: " << _bedslopegeodata.yorigin <<  std::endl;
285          osg::notify(osg::INFO) << "[SWWReader::setBedslopetexture] xpixel: " << _bedslopegeodata.xpixel <<  std::endl;
286          osg::notify(osg::INFO) << "[SWWReader::setBedslopetexture] ypixel: " << _bedslopegeodata.ypixel <<  std::endl;
287          osg::notify(osg::INFO) << "[SWWReader::setBedslopetexture] rotation: " << _bedslopegeodata.rotation <<  std::endl;
288      }
289   }
290
291}
292
293
294
295osg::Image* SWWReader::getBedslopeTexture()
296{
297   return osgDB::readImageFile( _bedslopetexture->c_str() );
298}
299
300
301
302osg::ref_ptr<osg::Vec2Array> SWWReader::getBedslopeTextureCoords()
303{
304   _bedslopetexcoords = new osg::Vec2Array;
305    if( _bedslopegeodata.hasData )  // georeferenced bedslope texture
306    {
307        float xorigin = _bedslopegeodata.xorigin;
308        float yorigin = _bedslopegeodata.yorigin;
309        float xrange = _bedslopegeodata.xpixel * _bedslopegeodata.xresolution;
310        float yrange = _bedslopegeodata.ypixel * _bedslopegeodata.yresolution;
311        for( unsigned int iv=0; iv < _npoints; iv++ )
312            _bedslopetexcoords->push_back( osg::Vec2( (_px[iv] + _xllcorner - xorigin)/xrange, 
313                                                      (_py[iv] + _yllcorner - yorigin)/yrange ) );
314    }
315    else
316    {
317        // decal'd using (x,y) locations scaled by extents into range [0,1]
318        for( unsigned int iv=0; iv < _npoints; iv++ )
319            _bedslopetexcoords->push_back( osg::Vec2( (_px[iv]-_xoffset)*_xscale, (_py[iv]-_yoffset)*_yscale ) );
320    }
321
322   return _bedslopetexcoords;
323}
324
325
326
327osg::ref_ptr<osg::Vec3Array> SWWReader::getStageVertexArray(unsigned int index)
328{
329    size_t start[2], count[2], iv;
330    const ptrdiff_t stride[2] = {1,1};
331    start[0] = index;
332    start[1] = 0;
333    count[0] = 1;
334    count[1] = _npoints;
335
336    // stage heights from netcdf file (x and y are same as bedslope)
337    int status = nc_get_vars_float (_ncid, _stageid, start, count, stride, _pstage);
338
339    // load stage vertex array, scaling and shifting vertices to lie in the unit cube
340    _stagevertices = new osg::Vec3Array;
341    for (iv=0; iv < _npoints; iv++)
342        _stagevertices->push_back( osg::Vec3( (_px[iv]-_xoffset)*_scale - _xcenter, 
343                                              (_py[iv]-_yoffset)*_scale - _ycenter, 
344                                              (_pstage[iv]-_zoffset)*_scale - _zcenter) );
345
346    // stage index, per primitive normal and centroid arrays
347    _stageprimitivenormals = new osg::Vec3Array;
348    _stagecentroids = new osg::Vec3Array;
349    osg::Vec3 v1b, v2b, v3b;
350    osg::Vec3 v1s, v2s, v3s;
351    osg::Vec3 side1, side2, nrm;
352    unsigned int v1index, v2index, v3index;
353
354    // over all stage triangles
355    for (iv=0; iv < _nvolumes; iv++)
356    {
357        v1index = _pvolumes[3*iv+0];
358        v2index = _pvolumes[3*iv+1];
359        v3index = _pvolumes[3*iv+2];
360
361        v1s = _stagevertices->at(v1index);
362        v2s = _stagevertices->at(v2index);
363        v3s = _stagevertices->at(v3index);
364        v1b = _bedslopevertices->at(v1index);
365        v2b = _bedslopevertices->at(v2index);
366        v3b = _bedslopevertices->at(v3index);
367
368        // current triangle primitive normal
369        side1 = v2s - v1s;
370        side2 = v3s - v2s;
371        nrm = side1^side2;
372        nrm.normalize();
373
374        // store centroid and primitive normal
375        _stagecentroids->push_back( (v1s+v2s+v3s)/3.0 );
376        _stageprimitivenormals->push_back( nrm );
377    }
378
379
380    // stage height above bedslope mapped as alpha value
381    //      alpha = min( a(h-hmin) + alphamin, alphamax),  h >= hmin
382    //      alpha = 0,                                     h < hmin
383    // where a = (alphamax-alphamin)/(hmax-hmin)
384    float alpha, height, alphascale;
385    alphascale = (_alphamax - _alphamin) / (_heightmax - _heightmin);
386    _stagecolors = new osg::Vec4Array;
387    for (iv=0; iv < _npoints; iv++)
388    {
389        height = _stagevertices->at(iv).z() - _bedslopevertices->at(iv).z();
390       
391        if (height < _heightmin)
392        {
393            alpha = 0.0;
394        } 
395        else 
396        {
397            alpha = alphascale * (height - _heightmin) + _alphamin;
398            if( alpha > _alphamax ) 
399                alpha = _alphamax;       
400        }
401        _stagecolors->push_back( osg::Vec4( 1.0, 1.0, 1.0, alpha ) );
402    }
403
404
405    // per-vertex normals calculated as average of primitive normals
406    // from contributing triangles
407    _stagevertexnormals = new osg::Vec3Array;
408    int num_shared_triangles, triangle_index;
409    for (iv=0; iv < _npoints; iv++)
410    {
411        nrm.set(0,0,0);
412        num_shared_triangles = _connectivity.at(iv).size();
413        for (int i=0; i<num_shared_triangles; i++ )
414        {
415            triangle_index = _connectivity.at(iv).at(i);
416            nrm += _stageprimitivenormals->at(triangle_index);
417        }
418
419        nrm = nrm / num_shared_triangles;  // average
420        nrm.normalize();
421        _stagevertexnormals->push_back(nrm);
422    }
423
424    return _stagevertices;
425}
426
427
428
429
430bool SWWReader::_statusHasError()
431{
432    bool haserror = false;  // assume success, trap failure
433
434    std::vector<int>::iterator iter;
435    for (iter=_status.begin(); iter != _status.end(); iter++)
436    {
437        if (*iter != NC_NOERR)
438        {
439            osg::notify(osg::WARN) << nc_strerror(*iter) <<  std::endl;
440            haserror = true;
441            nc_close(_ncid);
442            break;
443        }
444    }
445
446    // on return start gathering result values afresh
447    _status.clear();
448
449    return haserror;
450}
Note: See TracBrowser for help on using the repository browser.