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