1 : /******************************************************************************
2 : * $Id: ntf_raster.cpp 21977 2011-03-18 19:53:18Z 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 21977 2011-03-18 19:53:18Z 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 0 : 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 0 : 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 0 : 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", OFTReal );
245 0 : poFeatureDefn->AddFieldDefn( &oHeight );
246 :
247 0 : poReader = poReaderIn;
248 0 : poDS = poDSIn;
249 0 : poFilterGeom = NULL;
250 :
251 : pafColumn = (float *) CPLCalloc(sizeof(float),
252 0 : poReader->GetRasterYSize());
253 0 : iColumnOffset = -1;
254 0 : iCurrentFC = 0;
255 :
256 : /* -------------------------------------------------------------------- */
257 : /* Check for DEM subsampling, and compute total feature count */
258 : /* accordingly. */
259 : /* -------------------------------------------------------------------- */
260 0 : if( poDS->GetOption( "DEM_SAMPLE" ) == NULL )
261 0 : nDEMSample = 1;
262 : else
263 0 : nDEMSample = MAX(1,atoi(poDS->GetOption("DEM_SAMPLE")));
264 :
265 : nFeatureCount = (poReader->GetRasterXSize() / nDEMSample)
266 0 : * (poReader->GetRasterYSize() / nDEMSample);
267 0 : }
268 :
269 : /************************************************************************/
270 : /* ~OGRNTFRasterLayer() */
271 : /************************************************************************/
272 :
273 0 : OGRNTFRasterLayer::~OGRNTFRasterLayer()
274 :
275 : {
276 0 : CPLFree( pafColumn );
277 0 : if( poFeatureDefn )
278 0 : poFeatureDefn->Release();
279 :
280 0 : if( poFilterGeom != NULL )
281 0 : delete poFilterGeom;
282 0 : }
283 :
284 : /************************************************************************/
285 : /* SetSpatialFilter() */
286 : /************************************************************************/
287 :
288 0 : void OGRNTFRasterLayer::SetSpatialFilter( OGRGeometry * poGeomIn )
289 :
290 : {
291 0 : if( poFilterGeom != NULL )
292 : {
293 0 : delete poFilterGeom;
294 0 : poFilterGeom = NULL;
295 : }
296 :
297 0 : if( poGeomIn != NULL )
298 0 : poFilterGeom = poGeomIn->clone();
299 0 : }
300 :
301 : /************************************************************************/
302 : /* ResetReading() */
303 : /************************************************************************/
304 :
305 0 : void OGRNTFRasterLayer::ResetReading()
306 :
307 : {
308 0 : iCurrentFC = 0;
309 0 : }
310 :
311 : /************************************************************************/
312 : /* GetNextFeature() */
313 : /************************************************************************/
314 :
315 0 : OGRFeature *OGRNTFRasterLayer::GetNextFeature()
316 :
317 : {
318 0 : if( iCurrentFC == 0 )
319 0 : iCurrentFC = 1;
320 : else
321 : {
322 : int iReqColumn, iReqRow;
323 :
324 0 : iReqColumn = (iCurrentFC - 1) / poReader->GetRasterYSize();
325 0 : iReqRow = iCurrentFC - iReqColumn * poReader->GetRasterXSize() - 1;
326 :
327 0 : if( iReqRow + nDEMSample > poReader->GetRasterYSize() )
328 : {
329 0 : iReqRow = 0;
330 0 : iReqColumn += nDEMSample;
331 : }
332 : else
333 : {
334 0 : iReqRow += nDEMSample;
335 : }
336 :
337 : iCurrentFC = iReqColumn * poReader->GetRasterYSize()
338 0 : + iReqRow + 1;
339 : }
340 :
341 0 : return GetFeature( (long) iCurrentFC );
342 : }
343 :
344 : /************************************************************************/
345 : /* GetFeature() */
346 : /************************************************************************/
347 :
348 0 : OGRFeature *OGRNTFRasterLayer::GetFeature( long nFeatureId )
349 :
350 : {
351 : int iReqColumn, iReqRow;
352 :
353 : /* -------------------------------------------------------------------- */
354 : /* Is this in the range of legal feature ids (pixels)? */
355 : /* -------------------------------------------------------------------- */
356 0 : if( nFeatureId < 1
357 : || nFeatureId > poReader->GetRasterXSize()*poReader->GetRasterYSize() )
358 : {
359 0 : return NULL;
360 : }
361 :
362 : /* -------------------------------------------------------------------- */
363 : /* Do we need to load a different column. */
364 : /* -------------------------------------------------------------------- */
365 0 : iReqColumn = (nFeatureId - 1) / poReader->GetRasterYSize();
366 0 : iReqRow = nFeatureId - iReqColumn * poReader->GetRasterXSize() - 1;
367 :
368 0 : if( iReqColumn != iColumnOffset )
369 : {
370 0 : iColumnOffset = iReqColumn;
371 0 : if( poReader->ReadRasterColumn( iReqColumn, pafColumn ) != CE_None )
372 0 : return NULL;
373 : }
374 :
375 : /* -------------------------------------------------------------------- */
376 : /* Create a corresponding feature. */
377 : /* -------------------------------------------------------------------- */
378 0 : OGRFeature *poFeature = new OGRFeature( poFeatureDefn );
379 0 : double *padfGeoTransform = poReader->GetGeoTransform();
380 :
381 0 : poFeature->SetFID( nFeatureId );
382 :
383 : // NOTE: unusual use of GeoTransform - the pixel origin is the
384 : // bottom left corner!
385 : poFeature->SetGeometryDirectly(
386 0 : new OGRPoint( padfGeoTransform[0] + padfGeoTransform[1] * iReqColumn,
387 0 : padfGeoTransform[3] + padfGeoTransform[5] * iReqRow,
388 0 : pafColumn[iReqRow] ) );
389 0 : poFeature->SetField( 0, pafColumn[iReqRow] );
390 :
391 0 : return poFeature;
392 : }
393 :
394 : /************************************************************************/
395 : /* GetFeatureCount() */
396 : /* */
397 : /* If a spatial filter is in effect, we turn control over to */
398 : /* the generic counter. Otherwise we return the total count. */
399 : /* Eventually we should consider implementing a more efficient */
400 : /* way of counting features matching a spatial query. */
401 : /************************************************************************/
402 :
403 0 : int OGRNTFRasterLayer::GetFeatureCount( int bForce )
404 :
405 : {
406 0 : return nFeatureCount;
407 : }
408 :
409 : /************************************************************************/
410 : /* GetSpatialRef() */
411 : /************************************************************************/
412 :
413 0 : OGRSpatialReference *OGRNTFRasterLayer::GetSpatialRef()
414 :
415 : {
416 0 : return poDS->GetSpatialRef();
417 : }
418 :
419 : /************************************************************************/
420 : /* TestCapability() */
421 : /************************************************************************/
422 :
423 0 : int OGRNTFRasterLayer::TestCapability( const char * pszCap )
424 :
425 : {
426 0 : if( EQUAL(pszCap,OLCRandomRead) )
427 0 : return TRUE;
428 :
429 0 : else if( EQUAL(pszCap,OLCSequentialWrite)
430 : || EQUAL(pszCap,OLCRandomWrite) )
431 0 : return FALSE;
432 :
433 0 : else if( EQUAL(pszCap,OLCFastFeatureCount) )
434 0 : return TRUE;
435 :
436 0 : else if( EQUAL(pszCap,OLCFastSpatialFilter) )
437 0 : return FALSE;
438 :
439 : else
440 0 : return FALSE;
441 : }
442 :
443 :
|
