1 //*******************************************************************
2 //
3 // License: See top level LICENSE.txt file.
4 //
5 // Author: Ken Melero (kmelero@imagelinks.com)
6 // Orginally written by Jamie Moyers (jmoyers@geeks.com)
7 // Adapted from the package KDEM.
8 // Description: This class defines a DEM grid.
9 //
10 //********************************************************************
11 // $Id$
12
13 #include <ossim/support_data/ossimDemGrid.h>
14 #include <ossim/base/ossimNotify.h>
15 #include <ossim/base/ossimIoStream.h>
16 #include <ossim/support_data/ossimDemPoint.h>
17 #include <ossim/support_data/ossimDemUtil.h>
18
19
ossimDemGrid(ossim_float32 missingDataValue)20 ossimDemGrid::ossimDemGrid(ossim_float32 missingDataValue)
21 : _missDataVal(missingDataValue),
22 _width(0),
23 _height(0),
24 _grid(0),
25 _firstTime(true),
26 _curProfile(0)
27 {
28 }
29
~ossimDemGrid()30 ossimDemGrid::~ossimDemGrid()
31 {
32 if(_grid)
33 {
34 delete [] _grid;
35 _grid =0;
36 }
37 }
38
39 long
read(ossim::istream & dem,bool incrementalRead)40 ossimDemGrid::read(ossim::istream& dem, bool incrementalRead)
41 {
42 if (_firstTime)
43 {
44 dem >> _header;
45 }
46
47 long retval;
48 if (_header.getGroundRefSysCode() == 0) // Geographic
49 retval = fillGeographic(dem,incrementalRead);
50 else
51 retval = fillUTM(dem,incrementalRead); // This may not work if it's
52 // really in State Plane.
53
54 if (_firstTime)
55 _firstTime = false;
56
57 return retval;
58 }
59
60 ossimDemHeader const&
getHeader() const61 ossimDemGrid::getHeader() const
62 {
63 return _header;
64 }
65
66 long
getWidth() const67 ossimDemGrid::getWidth() const
68 {
69 return _width;
70 }
71
72 long
getHeight() const73 ossimDemGrid::getHeight() const
74 {
75 return _height;
76 }
77
getElevation(long x,long y) const78 ossim_float32 ossimDemGrid::getElevation(long x, long y) const
79 {
80 if (_grid == 0)
81 return _missDataVal;
82
83 if ((x < 0) || (y < 0) || (x >= _width) || (y >= _height))
84 return _missDataVal;
85
86 return _grid[(y * _width) + x];
87 }
88
getMissingDataValue() const89 ossim_float32 ossimDemGrid::getMissingDataValue() const
90 {
91 return _missDataVal;
92 }
93
fillGeographic(ossim::istream & dem,bool incrementalRead)94 long ossimDemGrid::fillGeographic(ossim::istream& dem,bool incrementalRead)
95 {
96 if (_firstTime) {
97 _curProfile = 0;
98 _width = _header.getProfileColumns();
99 }
100
101 while (_curProfile < _width)
102 {
103 _profiles.push_back(ossimDemProfile());
104 dem >> _profiles.back();
105 _curProfile++;
106 if (incrementalRead)
107 return _width - _curProfile + 1;
108 }
109
110 // Assume all profiles have as many elevations as the first.
111 _height = _profiles[0].getNumberOfElevations();
112 if(_grid)
113 {
114 delete [] _grid;
115 _grid =0;
116 }
117 _grid = new ossim_float32[_width * _height];
118
119 ossimDemPoint sw_corner = _profiles[0].getProfileLocation();
120 _northwest_x = sw_corner.getX();
121 _northwest_y = sw_corner.getY()
122 + ((_profiles[0].getNumberOfElevations() - 1) * _header.getSpatialResY());
123
124
125 unsigned int i,j;
126 for (i = 0; (int)i < _width; i++)
127 {
128 ossimDemElevationVector const& elev = _profiles[i].getElevations();
129 for (j = 0; j < elev.size(); j++)
130 {
131 setElevation(i, _height - j - 1, elev[j]);
132 }
133 }
134
135 _profiles.erase(_profiles.begin(), _profiles.end());
136
137 return 0;
138 }
139
140 long
fillUTM(ossim::istream & dem,bool incrementalRead)141 ossimDemGrid::fillUTM(ossim::istream& dem, bool incrementalRead)
142 {
143 // 7.5 UTM DEMs are small enough we can get away with doing this stupid...
144
145 unsigned int i;
146 unsigned int x,y;
147
148 if (_firstTime)
149 {
150 _curProfile = 0;
151 _width = _header.getProfileColumns();
152 }
153
154
155 while (_curProfile < _width)
156 {
157 _profiles.push_back(ossimDemProfile());
158 dem >> _profiles.back();
159 _curProfile++;
160 if (incrementalRead)
161 return _width - _curProfile + 1;
162 }
163
164 double dy = _header.getSpatialResY();
165
166 // Determine min and max Y values.
167 // Some DEMs can have profiles which do not have any
168 // elevations, and erroneous (x,y) values. I suspect
169 // these are probably illegal DEMs, but we'll try to
170 // do the right thing anyway.
171 ossimDemPoint curpoint;
172 i = 0;
173 while ((_profiles[i].getNumberOfElevations() == 0) &&
174 (i < _profiles.size()))
175 i++;
176 if (i < _profiles.size())
177 curpoint = _profiles[i].getProfileLocation();
178 else
179 {
180 // XXX This isn't the best way to handle this...
181 ossimNotify(ossimNotifyLevel_FATAL) << "FATAL ossimDemGrid::fillUTM: Yikes! All profiles have zero elevations!\n";
182 return(1);
183 }
184
185 double miny, maxy;
186 miny = curpoint.getY();
187 maxy = miny;
188 double profymin, profymax; // Min and max y values for current profile.
189 for (i = 0; i < _profiles.size(); i++)
190 {
191 if (_profiles[i].getNumberOfElevations() > 0)
192 {
193 curpoint = _profiles[i].getProfileLocation();
194 profymin = curpoint.getY();
195 profymax = profymin + ((_profiles[i].getNumberOfElevations() - 1) * dy);
196
197 if (profymin < miny)
198 miny = profymin;
199 if (profymax > maxy)
200 maxy = profymax;
201 }
202 }
203
204 // We now have minimum and maximum y values over all profiles in the DEM.
205 // Allocate a rectangular array large enough to hold them.
206
207 _height = static_cast<long>(((maxy - miny) / dy) + 1);
208 if(_grid)
209 {
210 delete [] _grid;
211 _grid =0;
212 }
213 _grid = new ossim_float32[_width * _height];
214
215 // Fill grid with the "missing data" value.
216 for (i = 0; (int)i < _width * _height; i++)
217 _grid[i] = _missDataVal;
218
219 ossimDemPoint sw_corner = _profiles[0].getProfileLocation();
220 _northwest_x = sw_corner.getX();
221 _northwest_y = maxy;
222
223
224 // Now, insert the elevations in the profiles in the appropriate place in the grid.
225
226 long startpos;
227 for (x = 0; (int)x < _width; x++)
228 {
229 ossimDemElevationVector const& elev = _profiles[x].getElevations();
230 curpoint = _profiles[x].getProfileLocation();
231 startpos = static_cast<long>((curpoint.getY() - miny) / dy);
232 for (y = 0; y < elev.size(); y++)
233 {
234 setElevation(x, _height - startpos - y - 1, elev[y]);
235 }
236 }
237
238 _profiles.erase(_profiles.begin(), _profiles.end());
239
240 return 0;
241 }
242
243 void
getGroundCoords(long x,long y,double & ground_x,double & ground_y)244 ossimDemGrid::getGroundCoords(long x, long y, double& ground_x,
245 double& ground_y)
246 {
247 ground_x = _northwest_x + (x * _header.getSpatialResX());
248 ground_y = _northwest_y - (y * _header.getSpatialResY());
249 }
250