1 : /******************************************************************************
2 : * $Id: ntf_raster.cpp 10645 2007-01-18 02:22:39Z warmerdam $
3 : *
4 : * Project: NTF Translator
5 : * Purpose: Handle UK Ordnance Survey Raster DTM products. Includes some
6 : * raster related methods from NTFFileReader and the implementation
7 : * of OGRNTFRasterLayer.
8 : * Author: Frank Warmerdam, warmerdam@pobox.com
9 : *
10 : ******************************************************************************
11 : * Copyright (c) 1999, Frank Warmerdam
12 : *
13 : * Permission is hereby granted, free of charge, to any person obtaining a
14 : * copy of this software and associated documentation files (the "Software"),
15 : * to deal in the Software without restriction, including without limitation
16 : * the rights to use, copy, modify, merge, publish, distribute, sublicense,
17 : * and/or sell copies of the Software, and to permit persons to whom the
18 : * Software is furnished to do so, subject to the following conditions:
19 : *
20 : * The above copyright notice and this permission notice shall be included
21 : * in all copies or substantial portions of the Software.
22 : *
23 : * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
24 : * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
25 : * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
26 : * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
27 : * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
28 : * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
29 : * DEALINGS IN THE SOFTWARE.
30 : ****************************************************************************/
31 :
32 : #include "ntf.h"
33 :
34 : CPL_CVSID("$Id: ntf_raster.cpp 10645 2007-01-18 02:22:39Z warmerdam $");
35 :
36 : /************************************************************************/
37 : /* ==================================================================== */
38 : /* NTFFileReader Raster Methods */
39 : /* ==================================================================== */
40 : /************************************************************************/
41 :
42 : /************************************************************************/
43 : /* IsRasterProduct() */
44 : /************************************************************************/
45 :
46 31226 : int NTFFileReader::IsRasterProduct()
47 :
48 : {
49 : return GetProductId() == NPC_LANDRANGER_DTM
50 31226 : || GetProductId() == NPC_LANDFORM_PROFILE_DTM;
51 : }
52 :
53 : /************************************************************************/
54 : /* EstablishRasterAccess() */
55 : /************************************************************************/
56 :
57 0 : void NTFFileReader::EstablishRasterAccess()
58 :
59 : {
60 : /* -------------------------------------------------------------------- */
61 : /* Read the type 50 record. */
62 : /* -------------------------------------------------------------------- */
63 : NTFRecord *poRecord;
64 :
65 0 : while( (poRecord = ReadRecord()) != NULL
66 : && poRecord->GetType() != NRT_GRIDHREC
67 : && poRecord->GetType() != NRT_VTR )
68 : {
69 0 : delete poRecord;
70 : }
71 :
72 0 : if( poRecord->GetType() != NRT_GRIDHREC )
73 : {
74 0 : delete poRecord;
75 : CPLError( CE_Failure, CPLE_AppDefined,
76 : "Unable to find GRIDHREC (type 50) record in what appears\n"
77 0 : "to be an NTF Raster DTM product." );
78 0 : return;
79 : }
80 :
81 : /* -------------------------------------------------------------------- */
82 : /* Parse if LANDRANGER_DTM */
83 : /* -------------------------------------------------------------------- */
84 0 : if( GetProductId() == NPC_LANDRANGER_DTM )
85 : {
86 0 : nRasterXSize = atoi(poRecord->GetField(13,16));
87 0 : nRasterYSize = atoi(poRecord->GetField(17,20));
88 :
89 : // NOTE: unusual use of GeoTransform - the pixel origin is the
90 : // bottom left corner!
91 0 : adfGeoTransform[0] = atoi(poRecord->GetField(25,34));
92 0 : adfGeoTransform[1] = 50;
93 0 : adfGeoTransform[2] = 0;
94 0 : adfGeoTransform[3] = atoi(poRecord->GetField(35,44));
95 0 : adfGeoTransform[4] = 0;
96 0 : adfGeoTransform[5] = 50;
97 :
98 0 : nRasterDataType = 3; /* GDT_Int16 */
99 : }
100 :
101 : /* -------------------------------------------------------------------- */
102 : /* Parse if LANDFORM_PROFILE_DTM */
103 : /* -------------------------------------------------------------------- */
104 0 : else if( GetProductId() == NPC_LANDFORM_PROFILE_DTM )
105 : {
106 0 : nRasterXSize = atoi(poRecord->GetField(23,30));
107 0 : nRasterYSize = atoi(poRecord->GetField(31,38));
108 :
109 : // NOTE: unusual use of GeoTransform - the pixel origin is the
110 : // bottom left corner!
111 : adfGeoTransform[0] = atoi(poRecord->GetField(13,17))
112 0 : + GetXOrigin();
113 0 : adfGeoTransform[1] = atoi(poRecord->GetField(39,42));
114 0 : adfGeoTransform[2] = 0;
115 : adfGeoTransform[3] = atoi(poRecord->GetField(18,22))
116 0 : + GetYOrigin();
117 0 : adfGeoTransform[4] = 0;
118 0 : adfGeoTransform[5] = atoi(poRecord->GetField(43,46));
119 :
120 0 : nRasterDataType = 3; /* GDT_Int16 */
121 : }
122 :
123 : /* -------------------------------------------------------------------- */
124 : /* Initialize column offsets table. */
125 : /* -------------------------------------------------------------------- */
126 0 : delete poRecord;
127 :
128 0 : panColumnOffset = (long *) CPLCalloc(sizeof(long),nRasterXSize);
129 :
130 0 : GetFPPos( panColumnOffset+0, NULL );
131 :
132 : /* -------------------------------------------------------------------- */
133 : /* Create an OGRSFLayer for this file readers raster points. */
134 : /* -------------------------------------------------------------------- */
135 0 : if( poDS != NULL )
136 : {
137 0 : poRasterLayer = new OGRNTFRasterLayer( poDS, this );
138 0 : poDS->AddLayer( poRasterLayer );
139 : }
140 : }
141 :
142 : /************************************************************************/
143 : /* ReadRasterColumn() */
144 : /************************************************************************/
145 :
146 0 : CPLErr NTFFileReader::ReadRasterColumn( int iColumn, float *pafElev )
147 :
148 : {
149 : /* -------------------------------------------------------------------- */
150 : /* If we don't already have the scanline offset of the previous */
151 : /* line, force reading of previous records to establish it. */
152 : /* -------------------------------------------------------------------- */
153 0 : if( panColumnOffset[iColumn] == 0 )
154 : {
155 : int iPrev;
156 :
157 0 : for( iPrev = 0; iPrev < iColumn-1; iPrev++ )
158 : {
159 0 : if( panColumnOffset[iPrev+1] == 0 )
160 : {
161 : CPLErr eErr;
162 :
163 0 : eErr = ReadRasterColumn( iPrev, NULL );
164 0 : if( eErr != CE_None )
165 0 : return eErr;
166 : }
167 : }
168 : }
169 :
170 : /* -------------------------------------------------------------------- */
171 : /* If the dataset isn't open, open it now. */
172 : /* -------------------------------------------------------------------- */
173 0 : if( GetFP() == NULL )
174 0 : Open();
175 :
176 : /* -------------------------------------------------------------------- */
177 : /* Read requested record. */
178 : /* -------------------------------------------------------------------- */
179 : NTFRecord *poRecord;
180 :
181 0 : SetFPPos( panColumnOffset[iColumn], iColumn );
182 0 : poRecord = ReadRecord();
183 :
184 0 : if( iColumn < nRasterXSize-1 )
185 : {
186 0 : GetFPPos( panColumnOffset+iColumn+1, NULL );
187 : }
188 :
189 : /* -------------------------------------------------------------------- */
190 : /* Handle LANDRANGER DTM columns. */
191 : /* -------------------------------------------------------------------- */
192 0 : if( pafElev != NULL && GetProductId() == NPC_LANDRANGER_DTM )
193 : {
194 : double dfVScale, dfVOffset;
195 :
196 0 : dfVOffset = atoi(poRecord->GetField(56,65));
197 0 : dfVScale = atoi(poRecord->GetField(66,75)) * 0.001;
198 :
199 0 : for( int iPixel = 0; iPixel < nRasterXSize; iPixel++ )
200 : {
201 : pafElev[iPixel] = (float) (dfVOffset + dfVScale *
202 0 : atoi(poRecord->GetField(84+iPixel*4,87+iPixel*4)));
203 : }
204 : }
205 :
206 : /* -------------------------------------------------------------------- */
207 : /* Handle PROFILE */
208 : /* -------------------------------------------------------------------- */
209 0 : else if( pafElev != NULL && GetProductId() == NPC_LANDFORM_PROFILE_DTM )
210 : {
211 0 : for( int iPixel = 0; iPixel < nRasterXSize; iPixel++ )
212 : {
213 : pafElev[iPixel] = (float)
214 0 : (atoi(poRecord->GetField(19+iPixel*5,23+iPixel*5)) * GetZMult());
215 : }
216 : }
217 :
218 0 : delete poRecord;
219 :
220 0 : return CE_None;
221 : }
222 :
223 : /************************************************************************/
224 : /* ==================================================================== */
225 : /* OGRNTFRasterLayer */
226 : /* ==================================================================== */
227 : /************************************************************************/
228 :
229 : /************************************************************************/
230 : /* OGRNTFRasterLayer */
231 : /************************************************************************/
232 :
233 : OGRNTFRasterLayer::OGRNTFRasterLayer( OGRNTFDataSource *poDSIn,
234 0 : NTFFileReader * poReaderIn )
235 :
236 : {
237 : char szLayerName[128];
238 :
239 0 : sprintf( szLayerName, "DTM_%s", poReaderIn->GetTileName() );
240 0 : poFeatureDefn = new OGRFeatureDefn( szLayerName );
241 0 : poFeatureDefn->Reference();
242 0 : poFeatureDefn->SetGeomType( wkbPoint25D );
243 :
244 0 : OGRFieldDefn oHeight( "HEIGHT", OFTInteger );
245 0 : oHeight.SetWidth(5);
246 0 : poFeatureDefn->AddFieldDefn( &oHeight );
247 :
248 0 : poReader = poReaderIn;
249 0 : poDS = poDSIn;
250 0 : poFilterGeom = NULL;
251 :
252 : pafColumn = (float *) CPLCalloc(sizeof(float),
253 0 : poReader->GetRasterYSize());
254 0 : iColumnOffset = -1;
255 0 : iCurrentFC = 0;
256 :
257 : /* -------------------------------------------------------------------- */
258 : /* Check for DEM subsampling, and compute total feature count */
259 : /* accordingly. */
260 : /* -------------------------------------------------------------------- */
261 0 : if( poDS->GetOption( "DEM_SAMPLE" ) == NULL )
262 0 : nDEMSample = 1;
263 : else
264 0 : nDEMSample = MAX(1,atoi(poDS->GetOption("DEM_SAMPLE")));
265 :
266 : nFeatureCount = (poReader->GetRasterXSize() / nDEMSample)
267 0 : * (poReader->GetRasterYSize() / nDEMSample);
268 0 : }
269 :
270 : /************************************************************************/
271 : /* ~OGRNTFRasterLayer() */
272 : /************************************************************************/
273 :
274 0 : OGRNTFRasterLayer::~OGRNTFRasterLayer()
275 :
276 : {
277 0 : CPLFree( pafColumn );
278 0 : if( poFeatureDefn )
279 0 : poFeatureDefn->Release();
280 :
281 0 : if( poFilterGeom != NULL )
282 0 : delete poFilterGeom;
283 0 : }
284 :
285 : /************************************************************************/
286 : /* SetSpatialFilter() */
287 : /************************************************************************/
288 :
289 0 : void OGRNTFRasterLayer::SetSpatialFilter( OGRGeometry * poGeomIn )
290 :
291 : {
292 0 : if( poFilterGeom != NULL )
293 : {
294 0 : delete poFilterGeom;
295 0 : poFilterGeom = NULL;
296 : }
297 :
298 0 : if( poGeomIn != NULL )
299 0 : poFilterGeom = poGeomIn->clone();
300 0 : }
301 :
302 : /************************************************************************/
303 : /* ResetReading() */
304 : /************************************************************************/
305 :
306 0 : void OGRNTFRasterLayer::ResetReading()
307 :
308 : {
309 0 : iCurrentFC = 0;
310 0 : }
311 :
312 : /************************************************************************/
313 : /* GetNextFeature() */
314 : /************************************************************************/
315 :
316 0 : OGRFeature *OGRNTFRasterLayer::GetNextFeature()
317 :
318 : {
319 0 : if( iCurrentFC == 0 )
320 0 : iCurrentFC = 1;
321 : else
322 : {
323 : int iReqColumn, iReqRow;
324 :
325 0 : iReqColumn = (iCurrentFC - 1) / poReader->GetRasterYSize();
326 0 : iReqRow = iCurrentFC - iReqColumn * poReader->GetRasterXSize() - 1;
327 :
328 0 : if( iReqRow + nDEMSample > poReader->GetRasterYSize() )
329 : {
330 0 : iReqRow = 0;
331 0 : iReqColumn += nDEMSample;
332 : }
333 : else
334 : {
335 0 : iReqRow += nDEMSample;
336 : }
337 :
338 : iCurrentFC = iReqColumn * poReader->GetRasterYSize()
339 0 : + iReqRow + 1;
340 : }
341 :
342 0 : return GetFeature( (long) iCurrentFC );
343 : }
344 :
345 : /************************************************************************/
346 : /* GetFeature() */
347 : /************************************************************************/
348 :
349 0 : OGRFeature *OGRNTFRasterLayer::GetFeature( long nFeatureId )
350 :
351 : {
352 : int iReqColumn, iReqRow;
353 :
354 : /* -------------------------------------------------------------------- */
355 : /* Is this in the range of legal feature ids (pixels)? */
356 : /* -------------------------------------------------------------------- */
357 0 : if( nFeatureId < 1
358 : || nFeatureId > poReader->GetRasterXSize()*poReader->GetRasterYSize() )
359 : {
360 0 : return NULL;
361 : }
362 :
363 : /* -------------------------------------------------------------------- */
364 : /* Do we need to load a different column. */
365 : /* -------------------------------------------------------------------- */
366 0 : iReqColumn = (nFeatureId - 1) / poReader->GetRasterYSize();
367 0 : iReqRow = nFeatureId - iReqColumn * poReader->GetRasterXSize() - 1;
368 :
369 0 : if( iReqColumn != iColumnOffset )
370 : {
371 0 : iColumnOffset = iReqColumn;
372 0 : if( poReader->ReadRasterColumn( iReqColumn, pafColumn ) != CE_None )
373 0 : return NULL;
374 : }
375 :
376 : /* -------------------------------------------------------------------- */
377 : /* Create a corresponding feature. */
378 : /* -------------------------------------------------------------------- */
379 0 : OGRFeature *poFeature = new OGRFeature( poFeatureDefn );
380 0 : double *padfGeoTransform = poReader->GetGeoTransform();
381 :
382 0 : poFeature->SetFID( nFeatureId );
383 :
384 : // NOTE: unusual use of GeoTransform - the pixel origin is the
385 : // bottom left corner!
386 : poFeature->SetGeometryDirectly(
387 : new OGRPoint( padfGeoTransform[0] + padfGeoTransform[1] * iReqColumn,
388 : padfGeoTransform[3] + padfGeoTransform[5] * iReqRow,
389 0 : pafColumn[iReqRow] ) );
390 0 : poFeature->SetField( 0, (int) pafColumn[iReqRow] );
391 :
392 0 : return poFeature;
393 : }
394 :
395 : /************************************************************************/
396 : /* GetFeatureCount() */
397 : /* */
398 : /* If a spatial filter is in effect, we turn control over to */
399 : /* the generic counter. Otherwise we return the total count. */
400 : /* Eventually we should consider implementing a more efficient */
401 : /* way of counting features matching a spatial query. */
402 : /************************************************************************/
403 :
404 0 : int OGRNTFRasterLayer::GetFeatureCount( int bForce )
405 :
406 : {
407 0 : return nFeatureCount;
408 : }
409 :
410 : /************************************************************************/
411 : /* GetSpatialRef() */
412 : /************************************************************************/
413 :
414 0 : OGRSpatialReference *OGRNTFRasterLayer::GetSpatialRef()
415 :
416 : {
417 0 : return poDS->GetSpatialRef();
418 : }
419 :
420 : /************************************************************************/
421 : /* TestCapability() */
422 : /************************************************************************/
423 :
424 0 : int OGRNTFRasterLayer::TestCapability( const char * pszCap )
425 :
426 : {
427 0 : if( EQUAL(pszCap,OLCRandomRead) )
428 0 : return TRUE;
429 :
430 0 : else if( EQUAL(pszCap,OLCSequentialWrite)
431 : || EQUAL(pszCap,OLCRandomWrite) )
432 0 : return FALSE;
433 :
434 0 : else if( EQUAL(pszCap,OLCFastFeatureCount) )
435 0 : return TRUE;
436 :
437 0 : else if( EQUAL(pszCap,OLCFastSpatialFilter) )
438 0 : return FALSE;
439 :
440 : else
441 0 : return FALSE;
442 : }
443 :
444 :
|
