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 57718 : int NDFDataset::Identify(GDALOpenInfo *poOpenInfo)
145 :
146 : {
147 : /* -------------------------------------------------------------------- */
148 : /* The user must select the header file (i.e. .H1). */
149 : /* -------------------------------------------------------------------- */
150 61363 : return poOpenInfo->nHeaderBytes >= 50 &&
151 3645 : (STARTS_WITH_CI(
152 : reinterpret_cast<const char *>(poOpenInfo->pabyHeader),
153 3641 : "NDF_REVISION=2") ||
154 3641 : STARTS_WITH_CI(
155 : reinterpret_cast<const char *>(poOpenInfo->pabyHeader),
156 57718 : "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 : CPLString osBasePath = CPLGetPathSafe(poOpenInfo->pszFilename);
285 2 : osFilename = CPLFormFilenameSafe(osBasePath, osFilename, nullptr);
286 : }
287 :
288 2 : VSILFILE *fpRaw = VSIFOpenL(osFilename, "rb");
289 2 : if (fpRaw == nullptr)
290 : {
291 0 : CPLError(CE_Failure, CPLE_AppDefined,
292 : "Failed to open band file: %s", osFilename.c_str());
293 0 : return nullptr;
294 : }
295 2 : poDS->papszExtraFiles = CSLAddString(poDS->papszExtraFiles, osFilename);
296 :
297 : auto poBand = RawRasterBand::Create(
298 2 : poDS.get(), iBand + 1, fpRaw, 0, 1, poDS->nRasterXSize, GDT_Byte,
299 : RawRasterBand::ByteOrder::ORDER_LITTLE_ENDIAN,
300 2 : RawRasterBand::OwnFP::YES);
301 2 : if (!poBand)
302 0 : return nullptr;
303 :
304 2 : snprintf(szKey, sizeof(szKey), "BAND%d_NAME", iBand + 1);
305 2 : poBand->SetDescription(poDS->Get(szKey, ""));
306 :
307 2 : snprintf(szKey, sizeof(szKey), "BAND%d_WAVELENGTHS", iBand + 1);
308 2 : poBand->SetMetadataItem("WAVELENGTHS", poDS->Get(szKey, ""));
309 :
310 2 : snprintf(szKey, sizeof(szKey), "BAND%d_RADIOMETRIC_GAINS/BIAS",
311 : iBand + 1);
312 2 : poBand->SetMetadataItem("RADIOMETRIC_GAINS_BIAS", poDS->Get(szKey, ""));
313 :
314 2 : poDS->SetBand(iBand + 1, std::move(poBand));
315 : }
316 :
317 : /* -------------------------------------------------------------------- */
318 : /* Fetch and parse USGS projection parameters. */
319 : /* -------------------------------------------------------------------- */
320 2 : double adfUSGSParams[15] = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
321 : const CPLStringList aosParamTokens(CSLTokenizeStringComplex(
322 4 : poDS->Get("USGS_PROJECTION_NUMBER", ""), ",", FALSE, TRUE));
323 :
324 2 : if (aosParamTokens.size() >= 15)
325 : {
326 0 : for (int i = 0; i < 15; i++)
327 0 : adfUSGSParams[i] = CPLAtof(aosParamTokens[i]);
328 : }
329 :
330 : /* -------------------------------------------------------------------- */
331 : /* Minimal georef support ... should add full USGS style */
332 : /* support at some point. */
333 : /* -------------------------------------------------------------------- */
334 2 : const int nUSGSProjection = atoi(poDS->Get("USGS_PROJECTION_NUMBER", ""));
335 2 : const int nZone = atoi(poDS->Get("USGS_MAP_ZONE", "0"));
336 :
337 4 : OGRSpatialReference oSRS;
338 2 : oSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
339 2 : oSRS.importFromUSGS(nUSGSProjection, nZone, adfUSGSParams, 12);
340 :
341 4 : const CPLString osDatum = poDS->Get("HORIZONTAL_DATUM", "");
342 2 : if (EQUAL(osDatum, "WGS84") || EQUAL(osDatum, "NAD83") ||
343 0 : EQUAL(osDatum, "NAD27"))
344 : {
345 2 : oSRS.SetWellKnownGeogCS(osDatum);
346 : }
347 0 : else if (STARTS_WITH_CI(osDatum, "NAD27"))
348 : {
349 0 : oSRS.SetWellKnownGeogCS("NAD27");
350 : }
351 : else
352 : {
353 0 : CPLError(CE_Warning, CPLE_AppDefined,
354 : "Unrecognized datum name in NLAPS/NDF file:%s, "
355 : "assuming WGS84.",
356 : osDatum.c_str());
357 0 : oSRS.SetWellKnownGeogCS("WGS84");
358 : }
359 :
360 2 : if (oSRS.GetRoot() != nullptr)
361 : {
362 2 : poDS->m_oSRS = std::move(oSRS);
363 : }
364 :
365 : /* -------------------------------------------------------------------- */
366 : /* Get geotransform. */
367 : /* -------------------------------------------------------------------- */
368 : char **papszUL =
369 2 : CSLTokenizeString2(poDS->Get("UPPER_LEFT_CORNER", ""), ",", 0);
370 : char **papszUR =
371 2 : CSLTokenizeString2(poDS->Get("UPPER_RIGHT_CORNER", ""), ",", 0);
372 : char **papszLL =
373 2 : CSLTokenizeString2(poDS->Get("LOWER_LEFT_CORNER", ""), ",", 0);
374 :
375 4 : if (CSLCount(papszUL) == 4 && CSLCount(papszUR) == 4 &&
376 2 : CSLCount(papszLL) == 4)
377 : {
378 2 : poDS->m_gt[0] = CPLAtof(papszUL[2]);
379 4 : poDS->m_gt[1] = (CPLAtof(papszUR[2]) - CPLAtof(papszUL[2])) /
380 2 : (poDS->nRasterXSize - 1);
381 4 : poDS->m_gt[2] = (CPLAtof(papszUR[3]) - CPLAtof(papszUL[3])) /
382 2 : (poDS->nRasterXSize - 1);
383 :
384 2 : poDS->m_gt[3] = CPLAtof(papszUL[3]);
385 4 : poDS->m_gt[4] = (CPLAtof(papszLL[2]) - CPLAtof(papszUL[2])) /
386 2 : (poDS->nRasterYSize - 1);
387 4 : poDS->m_gt[5] = (CPLAtof(papszLL[3]) - CPLAtof(papszUL[3])) /
388 2 : (poDS->nRasterYSize - 1);
389 :
390 : // Move origin up-left half a pixel.
391 2 : poDS->m_gt[0] -= poDS->m_gt[1] * 0.5;
392 2 : poDS->m_gt[0] -= poDS->m_gt[4] * 0.5;
393 2 : poDS->m_gt[3] -= poDS->m_gt[2] * 0.5;
394 2 : poDS->m_gt[3] -= poDS->m_gt[5] * 0.5;
395 : }
396 :
397 2 : CSLDestroy(papszUL);
398 2 : CSLDestroy(papszLL);
399 2 : CSLDestroy(papszUR);
400 :
401 : /* -------------------------------------------------------------------- */
402 : /* Initialize any PAM information. */
403 : /* -------------------------------------------------------------------- */
404 2 : poDS->SetDescription(poOpenInfo->pszFilename);
405 2 : poDS->TryLoadXML();
406 :
407 : /* -------------------------------------------------------------------- */
408 : /* Check for overviews. */
409 : /* -------------------------------------------------------------------- */
410 2 : poDS->oOvManager.Initialize(poDS.get(), poOpenInfo->pszFilename);
411 :
412 2 : return poDS.release();
413 : }
414 :
415 : /************************************************************************/
416 : /* GDALRegister_NDF() */
417 : /************************************************************************/
418 :
419 1935 : void GDALRegister_NDF()
420 :
421 : {
422 1935 : if (GDALGetDriverByName("NDF") != nullptr)
423 282 : return;
424 :
425 1653 : GDALDriver *poDriver = new GDALDriver();
426 :
427 1653 : poDriver->SetDescription("NDF");
428 1653 : poDriver->SetMetadataItem(GDAL_DCAP_RASTER, "YES");
429 1653 : poDriver->SetMetadataItem(GDAL_DMD_LONGNAME, "NLAPS Data Format");
430 1653 : poDriver->SetMetadataItem(GDAL_DMD_HELPTOPIC, "drivers/raster/ndf.html");
431 1653 : poDriver->SetMetadataItem(GDAL_DCAP_VIRTUALIO, "YES");
432 :
433 1653 : poDriver->pfnIdentify = NDFDataset::Identify;
434 1653 : poDriver->pfnOpen = NDFDataset::Open;
435 :
436 1653 : GetGDALDriverManager()->RegisterDriver(poDriver);
437 : }
|
