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