Line data Source code
1 : /******************************************************************************
2 : *
3 : * Project: NDF Driver
4 : * Purpose: Implementation of NLAPS Data Format read support.
5 : * Author: Frank Warmerdam, warmerdam@pobox.com
6 : *
7 : ******************************************************************************
8 : * Copyright (c) 2005, Frank Warmerdam
9 : * Copyright (c) 2008-2011, Even Rouault <even dot rouault at spatialys.com>
10 : *
11 : * SPDX-License-Identifier: MIT
12 : ****************************************************************************/
13 :
14 : #include "cpl_string.h"
15 : #include "gdal_frmts.h"
16 : #include "ogr_spatialref.h"
17 : #include "rawdataset.h"
18 :
19 : /************************************************************************/
20 : /* ==================================================================== */
21 : /* NDFDataset */
22 : /* ==================================================================== */
23 : /************************************************************************/
24 :
25 : class NDFDataset final : public RawDataset
26 : {
27 : GDALGeoTransform m_gt{};
28 :
29 : OGRSpatialReference m_oSRS{};
30 : char **papszExtraFiles;
31 :
32 : char **papszHeader;
33 : const char *Get(const char *pszKey, const char *pszDefault);
34 :
35 : CPL_DISALLOW_COPY_ASSIGN(NDFDataset)
36 :
37 : CPLErr Close() override;
38 :
39 : public:
40 : NDFDataset();
41 : ~NDFDataset() override;
42 :
43 : CPLErr GetGeoTransform(GDALGeoTransform >) const override;
44 :
45 1 : const OGRSpatialReference *GetSpatialRef() const override
46 : {
47 1 : return m_oSRS.IsEmpty() ? nullptr : &m_oSRS;
48 : }
49 :
50 : char **GetFileList(void) override;
51 :
52 : static GDALDataset *Open(GDALOpenInfo *);
53 : static int Identify(GDALOpenInfo *);
54 : };
55 :
56 : /************************************************************************/
57 : /* NDFDataset() */
58 : /************************************************************************/
59 :
60 2 : NDFDataset::NDFDataset() : papszExtraFiles(nullptr), papszHeader(nullptr)
61 : {
62 2 : m_oSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
63 2 : }
64 :
65 : /************************************************************************/
66 : /* ~NDFDataset() */
67 : /************************************************************************/
68 :
69 4 : NDFDataset::~NDFDataset()
70 :
71 : {
72 2 : NDFDataset::Close();
73 4 : }
74 :
75 : /************************************************************************/
76 : /* Close() */
77 : /************************************************************************/
78 :
79 4 : CPLErr NDFDataset::Close()
80 : {
81 4 : CPLErr eErr = CE_None;
82 4 : if (nOpenFlags != OPEN_FLAGS_CLOSED)
83 : {
84 2 : if (NDFDataset::FlushCache(true) != CE_None)
85 0 : eErr = CE_Failure;
86 :
87 2 : CSLDestroy(papszHeader);
88 2 : CSLDestroy(papszExtraFiles);
89 :
90 2 : if (GDALPamDataset::Close() != CE_None)
91 0 : eErr = CE_Failure;
92 : }
93 4 : return eErr;
94 : }
95 :
96 : /************************************************************************/
97 : /* GetGeoTransform() */
98 : /************************************************************************/
99 :
100 1 : CPLErr NDFDataset::GetGeoTransform(GDALGeoTransform >) const
101 :
102 : {
103 1 : gt = m_gt;
104 1 : return CE_None;
105 : }
106 :
107 : /************************************************************************/
108 : /* Get() */
109 : /* */
110 : /* Fetch a value from the header by keyword. */
111 : /************************************************************************/
112 :
113 26 : const char *NDFDataset::Get(const char *pszKey, const char *pszDefault)
114 :
115 : {
116 26 : const char *pszResult = CSLFetchNameValue(papszHeader, pszKey);
117 :
118 26 : if (pszResult == nullptr)
119 0 : return pszDefault;
120 :
121 26 : return pszResult;
122 : }
123 :
124 : /************************************************************************/
125 : /* GetFileList() */
126 : /************************************************************************/
127 :
128 1 : char **NDFDataset::GetFileList()
129 :
130 : {
131 : // Main data file, etc.
132 1 : char **papszFileList = GDALPamDataset::GetFileList();
133 :
134 : // Header file.
135 1 : papszFileList = CSLInsertStrings(papszFileList, -1, papszExtraFiles);
136 :
137 1 : return papszFileList;
138 : }
139 :
140 : /************************************************************************/
141 : /* Identify() */
142 : /************************************************************************/
143 :
144 58238 : int NDFDataset::Identify(GDALOpenInfo *poOpenInfo)
145 :
146 : {
147 : /* -------------------------------------------------------------------- */
148 : /* The user must select the header file (i.e. .H1). */
149 : /* -------------------------------------------------------------------- */
150 61958 : return poOpenInfo->nHeaderBytes >= 50 &&
151 3720 : (STARTS_WITH_CI(
152 : reinterpret_cast<const char *>(poOpenInfo->pabyHeader),
153 3716 : "NDF_REVISION=2") ||
154 3716 : STARTS_WITH_CI(
155 : reinterpret_cast<const char *>(poOpenInfo->pabyHeader),
156 58238 : "NDF_REVISION=0"));
157 : }
158 :
159 : /************************************************************************/
160 : /* Open() */
161 : /************************************************************************/
162 :
163 2 : GDALDataset *NDFDataset::Open(GDALOpenInfo *poOpenInfo)
164 :
165 : {
166 : /* -------------------------------------------------------------------- */
167 : /* The user must select the header file (i.e. .H1). */
168 : /* -------------------------------------------------------------------- */
169 2 : if (!Identify(poOpenInfo) || poOpenInfo->fpL == nullptr)
170 0 : return nullptr;
171 :
172 : /* -------------------------------------------------------------------- */
173 : /* Confirm the requested access is supported. */
174 : /* -------------------------------------------------------------------- */
175 2 : if (poOpenInfo->eAccess == GA_Update)
176 : {
177 0 : ReportUpdateNotSupportedByDriver("NDF");
178 0 : return nullptr;
179 : }
180 : /* -------------------------------------------------------------------- */
181 : /* Read and process the header into a local name/value */
182 : /* stringlist. We just take off the trailing semicolon. The */
183 : /* keyword is already separated from the value by an equal */
184 : /* sign. */
185 : /* -------------------------------------------------------------------- */
186 :
187 2 : const char *pszLine = nullptr;
188 2 : const int nHeaderMax = 1000;
189 2 : int nHeaderLines = 0;
190 : char **papszHeader =
191 2 : static_cast<char **>(CPLMalloc(sizeof(char *) * (nHeaderMax + 1)));
192 :
193 210 : while (nHeaderLines < nHeaderMax &&
194 212 : (pszLine = CPLReadLineL(poOpenInfo->fpL)) != nullptr &&
195 106 : !EQUAL(pszLine, "END_OF_HDR;"))
196 : {
197 104 : if (strstr(pszLine, "=") == nullptr)
198 0 : break;
199 :
200 104 : char *pszFixed = CPLStrdup(pszLine);
201 104 : if (pszFixed[strlen(pszFixed) - 1] == ';')
202 104 : pszFixed[strlen(pszFixed) - 1] = '\0';
203 :
204 104 : papszHeader[nHeaderLines++] = pszFixed;
205 104 : papszHeader[nHeaderLines] = nullptr;
206 : }
207 2 : CPL_IGNORE_RET_VAL(VSIFCloseL(poOpenInfo->fpL));
208 2 : poOpenInfo->fpL = nullptr;
209 :
210 2 : if (CSLFetchNameValue(papszHeader, "PIXELS_PER_LINE") == nullptr ||
211 2 : CSLFetchNameValue(papszHeader, "LINES_PER_DATA_FILE") == nullptr ||
212 6 : CSLFetchNameValue(papszHeader, "BITS_PER_PIXEL") == nullptr ||
213 2 : CSLFetchNameValue(papszHeader, "PIXEL_FORMAT") == nullptr)
214 : {
215 0 : CPLError(CE_Failure, CPLE_AppDefined,
216 : "Dataset appears to be NDF but is missing a required field.");
217 0 : CSLDestroy(papszHeader);
218 0 : return nullptr;
219 : }
220 :
221 4 : if (!EQUAL(CSLFetchNameValue(papszHeader, "PIXEL_FORMAT"), "BYTE") ||
222 2 : !EQUAL(CSLFetchNameValue(papszHeader, "BITS_PER_PIXEL"), "8"))
223 : {
224 0 : CPLError(CE_Failure, CPLE_AppDefined,
225 : "Currently NDF driver supports only 8bit BYTE format.");
226 0 : CSLDestroy(papszHeader);
227 0 : return nullptr;
228 : }
229 :
230 : /* -------------------------------------------------------------------- */
231 : /* Confirm the requested access is supported. */
232 : /* -------------------------------------------------------------------- */
233 2 : if (poOpenInfo->eAccess == GA_Update)
234 : {
235 0 : CSLDestroy(papszHeader);
236 0 : ReportUpdateNotSupportedByDriver("NDF");
237 0 : return nullptr;
238 : }
239 :
240 : /* -------------------------------------------------------------------- */
241 : /* Create a corresponding GDALDataset. */
242 : /* -------------------------------------------------------------------- */
243 4 : auto poDS = std::make_unique<NDFDataset>();
244 2 : poDS->papszHeader = papszHeader;
245 :
246 2 : poDS->nRasterXSize = atoi(poDS->Get("PIXELS_PER_LINE", ""));
247 2 : poDS->nRasterYSize = atoi(poDS->Get("LINES_PER_DATA_FILE", ""));
248 :
249 : /* -------------------------------------------------------------------- */
250 : /* Create a raw raster band for each file. */
251 : /* -------------------------------------------------------------------- */
252 : const char *pszBand =
253 2 : CSLFetchNameValue(papszHeader, "NUMBER_OF_BANDS_IN_VOLUME");
254 2 : if (pszBand == nullptr)
255 : {
256 0 : CPLError(CE_Failure, CPLE_AppDefined, "Cannot find band count");
257 0 : return nullptr;
258 : }
259 2 : const int nBands = atoi(pszBand);
260 :
261 4 : if (!GDALCheckDatasetDimensions(poDS->nRasterXSize, poDS->nRasterYSize) ||
262 2 : !GDALCheckBandCount(nBands, FALSE))
263 : {
264 0 : return nullptr;
265 : }
266 :
267 4 : for (int iBand = 0; iBand < nBands; iBand++)
268 : {
269 : char szKey[100];
270 2 : snprintf(szKey, sizeof(szKey), "BAND%d_FILENAME", iBand + 1);
271 2 : CPLString osFilename = poDS->Get(szKey, "");
272 :
273 : // NDF1 file do not include the band filenames.
274 2 : if (osFilename.empty())
275 : {
276 : char szBandExtension[15];
277 0 : snprintf(szBandExtension, sizeof(szBandExtension), "I%d",
278 : iBand + 1);
279 : osFilename =
280 0 : CPLResetExtensionSafe(poOpenInfo->pszFilename, szBandExtension);
281 : }
282 : else
283 : {
284 2 : if (CPLHasPathTraversal(osFilename.c_str()))
285 : {
286 0 : CPLError(CE_Failure, CPLE_AppDefined,
287 : "Path traversal detected in %s", osFilename.c_str());
288 0 : return nullptr;
289 : }
290 2 : CPLString osBasePath = CPLGetPathSafe(poOpenInfo->pszFilename);
291 2 : osFilename = CPLFormFilenameSafe(osBasePath, osFilename, nullptr);
292 : }
293 :
294 2 : VSILFILE *fpRaw = VSIFOpenL(osFilename, "rb");
295 2 : if (fpRaw == nullptr)
296 : {
297 0 : CPLError(CE_Failure, CPLE_AppDefined,
298 : "Failed to open band file: %s", osFilename.c_str());
299 0 : return nullptr;
300 : }
301 2 : poDS->papszExtraFiles = CSLAddString(poDS->papszExtraFiles, osFilename);
302 :
303 : auto poBand = RawRasterBand::Create(
304 2 : poDS.get(), iBand + 1, fpRaw, 0, 1, poDS->nRasterXSize, GDT_Byte,
305 : RawRasterBand::ByteOrder::ORDER_LITTLE_ENDIAN,
306 2 : RawRasterBand::OwnFP::YES);
307 2 : if (!poBand)
308 0 : return nullptr;
309 :
310 2 : snprintf(szKey, sizeof(szKey), "BAND%d_NAME", iBand + 1);
311 2 : poBand->SetDescription(poDS->Get(szKey, ""));
312 :
313 2 : snprintf(szKey, sizeof(szKey), "BAND%d_WAVELENGTHS", iBand + 1);
314 2 : poBand->SetMetadataItem("WAVELENGTHS", poDS->Get(szKey, ""));
315 :
316 2 : snprintf(szKey, sizeof(szKey), "BAND%d_RADIOMETRIC_GAINS/BIAS",
317 : iBand + 1);
318 2 : poBand->SetMetadataItem("RADIOMETRIC_GAINS_BIAS", poDS->Get(szKey, ""));
319 :
320 2 : poDS->SetBand(iBand + 1, std::move(poBand));
321 : }
322 :
323 : /* -------------------------------------------------------------------- */
324 : /* Fetch and parse USGS projection parameters. */
325 : /* -------------------------------------------------------------------- */
326 2 : double adfUSGSParams[15] = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
327 : const CPLStringList aosParamTokens(CSLTokenizeStringComplex(
328 4 : poDS->Get("USGS_PROJECTION_NUMBER", ""), ",", FALSE, TRUE));
329 :
330 2 : if (aosParamTokens.size() >= 15)
331 : {
332 0 : for (int i = 0; i < 15; i++)
333 0 : adfUSGSParams[i] = CPLAtof(aosParamTokens[i]);
334 : }
335 :
336 : /* -------------------------------------------------------------------- */
337 : /* Minimal georef support ... should add full USGS style */
338 : /* support at some point. */
339 : /* -------------------------------------------------------------------- */
340 2 : const int nUSGSProjection = atoi(poDS->Get("USGS_PROJECTION_NUMBER", ""));
341 2 : const int nZone = atoi(poDS->Get("USGS_MAP_ZONE", "0"));
342 :
343 4 : OGRSpatialReference oSRS;
344 2 : oSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
345 2 : oSRS.importFromUSGS(nUSGSProjection, nZone, adfUSGSParams, 12);
346 :
347 4 : const CPLString osDatum = poDS->Get("HORIZONTAL_DATUM", "");
348 2 : if (EQUAL(osDatum, "WGS84") || EQUAL(osDatum, "NAD83") ||
349 0 : EQUAL(osDatum, "NAD27"))
350 : {
351 2 : oSRS.SetWellKnownGeogCS(osDatum);
352 : }
353 0 : else if (STARTS_WITH_CI(osDatum, "NAD27"))
354 : {
355 0 : oSRS.SetWellKnownGeogCS("NAD27");
356 : }
357 : else
358 : {
359 0 : CPLError(CE_Warning, CPLE_AppDefined,
360 : "Unrecognized datum name in NLAPS/NDF file:%s, "
361 : "assuming WGS84.",
362 : osDatum.c_str());
363 0 : oSRS.SetWellKnownGeogCS("WGS84");
364 : }
365 :
366 2 : if (oSRS.GetRoot() != nullptr)
367 : {
368 2 : poDS->m_oSRS = std::move(oSRS);
369 : }
370 :
371 : /* -------------------------------------------------------------------- */
372 : /* Get geotransform. */
373 : /* -------------------------------------------------------------------- */
374 : char **papszUL =
375 2 : CSLTokenizeString2(poDS->Get("UPPER_LEFT_CORNER", ""), ",", 0);
376 : char **papszUR =
377 2 : CSLTokenizeString2(poDS->Get("UPPER_RIGHT_CORNER", ""), ",", 0);
378 : char **papszLL =
379 2 : CSLTokenizeString2(poDS->Get("LOWER_LEFT_CORNER", ""), ",", 0);
380 :
381 4 : if (CSLCount(papszUL) == 4 && CSLCount(papszUR) == 4 &&
382 2 : CSLCount(papszLL) == 4)
383 : {
384 2 : poDS->m_gt[0] = CPLAtof(papszUL[2]);
385 4 : poDS->m_gt[1] = (CPLAtof(papszUR[2]) - CPLAtof(papszUL[2])) /
386 2 : (poDS->nRasterXSize - 1);
387 4 : poDS->m_gt[2] = (CPLAtof(papszUR[3]) - CPLAtof(papszUL[3])) /
388 2 : (poDS->nRasterXSize - 1);
389 :
390 2 : poDS->m_gt[3] = CPLAtof(papszUL[3]);
391 4 : poDS->m_gt[4] = (CPLAtof(papszLL[2]) - CPLAtof(papszUL[2])) /
392 2 : (poDS->nRasterYSize - 1);
393 4 : poDS->m_gt[5] = (CPLAtof(papszLL[3]) - CPLAtof(papszUL[3])) /
394 2 : (poDS->nRasterYSize - 1);
395 :
396 : // Move origin up-left half a pixel.
397 2 : poDS->m_gt[0] -= poDS->m_gt[1] * 0.5;
398 2 : poDS->m_gt[0] -= poDS->m_gt[4] * 0.5;
399 2 : poDS->m_gt[3] -= poDS->m_gt[2] * 0.5;
400 2 : poDS->m_gt[3] -= poDS->m_gt[5] * 0.5;
401 : }
402 :
403 2 : CSLDestroy(papszUL);
404 2 : CSLDestroy(papszLL);
405 2 : CSLDestroy(papszUR);
406 :
407 : /* -------------------------------------------------------------------- */
408 : /* Initialize any PAM information. */
409 : /* -------------------------------------------------------------------- */
410 2 : poDS->SetDescription(poOpenInfo->pszFilename);
411 2 : poDS->TryLoadXML();
412 :
413 : /* -------------------------------------------------------------------- */
414 : /* Check for overviews. */
415 : /* -------------------------------------------------------------------- */
416 2 : poDS->oOvManager.Initialize(poDS.get(), poOpenInfo->pszFilename);
417 :
418 2 : return poDS.release();
419 : }
420 :
421 : /************************************************************************/
422 : /* GDALRegister_NDF() */
423 : /************************************************************************/
424 :
425 2024 : void GDALRegister_NDF()
426 :
427 : {
428 2024 : if (GDALGetDriverByName("NDF") != nullptr)
429 283 : return;
430 :
431 1741 : GDALDriver *poDriver = new GDALDriver();
432 :
433 1741 : poDriver->SetDescription("NDF");
434 1741 : poDriver->SetMetadataItem(GDAL_DCAP_RASTER, "YES");
435 1741 : poDriver->SetMetadataItem(GDAL_DMD_LONGNAME, "NLAPS Data Format");
436 1741 : poDriver->SetMetadataItem(GDAL_DMD_HELPTOPIC, "drivers/raster/ndf.html");
437 1741 : poDriver->SetMetadataItem(GDAL_DCAP_VIRTUALIO, "YES");
438 :
439 1741 : poDriver->pfnIdentify = NDFDataset::Identify;
440 1741 : poDriver->pfnOpen = NDFDataset::Open;
441 :
442 1741 : GetGDALDriverManager()->RegisterDriver(poDriver);
443 : }
|