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 : /* ResetReading() */
274 : /************************************************************************/
275 :
276 0 : void OGRNTFRasterLayer::ResetReading()
277 :
278 : {
279 0 : iCurrentFC = 1;
280 0 : }
281 :
282 : /************************************************************************/
283 : /* GetNextFeature() */
284 : /************************************************************************/
285 :
286 0 : OGRFeature *OGRNTFRasterLayer::GetNextFeature()
287 :
288 : {
289 0 : if (iCurrentFC > static_cast<GIntBig>(poReader->GetRasterXSize()) *
290 0 : poReader->GetRasterYSize())
291 : {
292 0 : return nullptr;
293 : }
294 :
295 0 : OGRFeature *poFeature = GetFeature(iCurrentFC);
296 :
297 : int iReqColumn, iReqRow;
298 :
299 0 : iReqColumn =
300 0 : static_cast<int>((iCurrentFC - 1) / poReader->GetRasterYSize());
301 0 : iReqRow = static_cast<int>(
302 0 : iCurrentFC -
303 0 : static_cast<GIntBig>(iReqColumn) * poReader->GetRasterYSize() - 1);
304 :
305 0 : if (iReqRow + nDEMSample > poReader->GetRasterYSize())
306 : {
307 0 : iReqRow = 0;
308 0 : iReqColumn += nDEMSample;
309 : }
310 : else
311 : {
312 0 : iReqRow += nDEMSample;
313 : }
314 :
315 0 : iCurrentFC = static_cast<GIntBig>(iReqColumn) * poReader->GetRasterYSize() +
316 0 : iReqRow + 1;
317 :
318 0 : return poFeature;
319 : }
320 :
321 : /************************************************************************/
322 : /* GetFeature() */
323 : /************************************************************************/
324 :
325 0 : OGRFeature *OGRNTFRasterLayer::GetFeature(GIntBig nFeatureId)
326 :
327 : {
328 : int iReqColumn, iReqRow;
329 :
330 : /* -------------------------------------------------------------------- */
331 : /* Is this in the range of legal feature ids (pixels)? */
332 : /* -------------------------------------------------------------------- */
333 0 : if (nFeatureId < 1 ||
334 0 : nFeatureId > static_cast<GIntBig>(poReader->GetRasterXSize()) *
335 0 : poReader->GetRasterYSize())
336 : {
337 0 : return nullptr;
338 : }
339 :
340 : /* -------------------------------------------------------------------- */
341 : /* Do we need to load a different column. */
342 : /* -------------------------------------------------------------------- */
343 0 : iReqColumn =
344 0 : static_cast<int>((nFeatureId - 1) / poReader->GetRasterYSize());
345 0 : iReqRow = static_cast<int>(
346 0 : nFeatureId -
347 0 : static_cast<GIntBig>(iReqColumn) * poReader->GetRasterYSize() - 1);
348 :
349 0 : if (iReqColumn != iColumnOffset)
350 : {
351 0 : iColumnOffset = iReqColumn;
352 0 : if (poReader->ReadRasterColumn(iReqColumn, pafColumn) != CE_None)
353 0 : return nullptr;
354 : }
355 0 : if (iReqRow < 0 || iReqRow >= poReader->GetRasterYSize())
356 0 : return nullptr;
357 :
358 : /* -------------------------------------------------------------------- */
359 : /* Create a corresponding feature. */
360 : /* -------------------------------------------------------------------- */
361 0 : OGRFeature *poFeature = new OGRFeature(poFeatureDefn);
362 0 : double *padfGeoTransform = poReader->GetGeoTransform();
363 :
364 0 : poFeature->SetFID(nFeatureId);
365 :
366 : // NOTE: unusual use of GeoTransform - the pixel origin is the
367 : // bottom left corner!
368 0 : poFeature->SetGeometryDirectly(
369 0 : new OGRPoint(padfGeoTransform[0] + padfGeoTransform[1] * iReqColumn,
370 0 : padfGeoTransform[3] + padfGeoTransform[5] * iReqRow,
371 0 : pafColumn[iReqRow]));
372 0 : poFeature->SetField(0, pafColumn[iReqRow]);
373 :
374 0 : return poFeature;
375 : }
376 :
377 : /************************************************************************/
378 : /* GetFeatureCount() */
379 : /* */
380 : /* If a spatial filter is in effect, we turn control over to */
381 : /* the generic counter. Otherwise we return the total count. */
382 : /* Eventually we should consider implementing a more efficient */
383 : /* way of counting features matching a spatial query. */
384 : /************************************************************************/
385 :
386 0 : GIntBig OGRNTFRasterLayer::GetFeatureCount(CPL_UNUSED int bForce)
387 : {
388 0 : return nFeatureCount;
389 : }
390 :
391 : /************************************************************************/
392 : /* TestCapability() */
393 : /************************************************************************/
394 :
395 0 : int OGRNTFRasterLayer::TestCapability(const char *pszCap)
396 :
397 : {
398 0 : if (EQUAL(pszCap, OLCRandomRead))
399 0 : return TRUE;
400 :
401 0 : else if (EQUAL(pszCap, OLCFastFeatureCount))
402 0 : return TRUE;
403 :
404 0 : return FALSE;
405 : }
|
