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