Line data Source code
1 : /******************************************************************************
2 : *
3 : * Project: ERMapper .ers Driver
4 : * Purpose: Implementation of .ers driver.
5 : * Author: Frank Warmerdam, warmerdam@pobox.com
6 : *
7 : ******************************************************************************
8 : * Copyright (c) 2007, Frank Warmerdam <warmerdam@pobox.com>
9 : * Copyright (c) 2008-2013, Even Rouault <even dot rouault at spatialys.com>
10 : *
11 : * SPDX-License-Identifier: MIT
12 : ****************************************************************************/
13 :
14 : #include "cpl_string.h"
15 : #include "ershdrnode.h"
16 : #include "gdal_frmts.h"
17 : #include "gdal_proxy.h"
18 : #include "ogr_spatialref.h"
19 : #include "rawdataset.h"
20 :
21 : #include <limits>
22 :
23 : /************************************************************************/
24 : /* ==================================================================== */
25 : /* ERSDataset */
26 : /* ==================================================================== */
27 : /************************************************************************/
28 :
29 : class ERSRasterBand;
30 :
31 : class ERSDataset final : public RawDataset
32 : {
33 : friend class ERSRasterBand;
34 :
35 : VSILFILE *fpImage; // Image data file.
36 : GDALDataset *poDepFile;
37 :
38 : int bGotTransform;
39 : double adfGeoTransform[6];
40 : OGRSpatialReference m_oSRS{};
41 :
42 : CPLString osRawFilename;
43 :
44 : int bHDRDirty;
45 : ERSHdrNode *poHeader;
46 :
47 : const char *Find(const char *, const char *);
48 :
49 : int nGCPCount;
50 : GDAL_GCP *pasGCPList;
51 : OGRSpatialReference m_oGCPSRS{};
52 :
53 : void ReadGCPs();
54 :
55 : int bHasNoDataValue;
56 : double dfNoDataValue;
57 :
58 : CPLString osProj, osProjForced;
59 : CPLString osDatum, osDatumForced;
60 : CPLString osUnits, osUnitsForced;
61 : void WriteProjectionInfo(const char *pszProj, const char *pszDatum,
62 : const char *pszUnits);
63 :
64 : CPLStringList oERSMetadataList;
65 :
66 : protected:
67 : int CloseDependentDatasets() override;
68 :
69 : CPLErr Close() override;
70 :
71 : public:
72 : ERSDataset();
73 : ~ERSDataset() override;
74 :
75 : CPLErr FlushCache(bool bAtClosing) override;
76 : CPLErr GetGeoTransform(double *padfTransform) override;
77 : CPLErr SetGeoTransform(double *padfTransform) override;
78 : const OGRSpatialReference *GetSpatialRef() const override;
79 : CPLErr SetSpatialRef(const OGRSpatialReference *poSRS) override;
80 : char **GetFileList(void) override;
81 :
82 : int GetGCPCount() override;
83 : const OGRSpatialReference *GetGCPSpatialRef() const override;
84 : const GDAL_GCP *GetGCPs() override;
85 : CPLErr SetGCPs(int nGCPCountIn, const GDAL_GCP *pasGCPListIn,
86 : const OGRSpatialReference *poSRS) override;
87 :
88 : char **GetMetadataDomainList() override;
89 : const char *GetMetadataItem(const char *pszName,
90 : const char *pszDomain = "") override;
91 : char **GetMetadata(const char *pszDomain = "") override;
92 :
93 : static GDALDataset *Open(GDALOpenInfo *);
94 : static int Identify(GDALOpenInfo *);
95 : static GDALDataset *Create(const char *pszFilename, int nXSize, int nYSize,
96 : int nBandsIn, GDALDataType eType,
97 : char **papszParamList);
98 : };
99 :
100 : /************************************************************************/
101 : /* ERSDataset() */
102 : /************************************************************************/
103 :
104 62 : ERSDataset::ERSDataset()
105 : : fpImage(nullptr), poDepFile(nullptr), bGotTransform(FALSE),
106 : bHDRDirty(FALSE), poHeader(nullptr), nGCPCount(0), pasGCPList(nullptr),
107 62 : bHasNoDataValue(FALSE), dfNoDataValue(0.0)
108 : {
109 62 : m_oSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
110 62 : m_oGCPSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
111 62 : adfGeoTransform[0] = 0.0;
112 62 : adfGeoTransform[1] = 1.0;
113 62 : adfGeoTransform[2] = 0.0;
114 62 : adfGeoTransform[3] = 0.0;
115 62 : adfGeoTransform[4] = 0.0;
116 62 : adfGeoTransform[5] = 1.0;
117 62 : }
118 :
119 : /************************************************************************/
120 : /* ~ERSDataset() */
121 : /************************************************************************/
122 :
123 124 : ERSDataset::~ERSDataset()
124 :
125 : {
126 62 : ERSDataset::Close();
127 124 : }
128 :
129 : /************************************************************************/
130 : /* Close() */
131 : /************************************************************************/
132 :
133 123 : CPLErr ERSDataset::Close()
134 : {
135 123 : CPLErr eErr = CE_None;
136 123 : if (nOpenFlags != OPEN_FLAGS_CLOSED)
137 : {
138 62 : if (ERSDataset::FlushCache(true) != CE_None)
139 0 : eErr = CE_Failure;
140 :
141 62 : if (fpImage != nullptr)
142 : {
143 60 : VSIFCloseL(fpImage);
144 : }
145 :
146 62 : ERSDataset::CloseDependentDatasets();
147 :
148 62 : if (nGCPCount > 0)
149 : {
150 3 : GDALDeinitGCPs(nGCPCount, pasGCPList);
151 3 : CPLFree(pasGCPList);
152 : }
153 :
154 62 : if (poHeader != nullptr)
155 62 : delete poHeader;
156 :
157 62 : if (GDALPamDataset::Close() != CE_None)
158 0 : eErr = CE_Failure;
159 : }
160 123 : return eErr;
161 : }
162 :
163 : /************************************************************************/
164 : /* CloseDependentDatasets() */
165 : /************************************************************************/
166 :
167 62 : int ERSDataset::CloseDependentDatasets()
168 : {
169 62 : int bHasDroppedRef = GDALPamDataset::CloseDependentDatasets();
170 :
171 62 : if (poDepFile != nullptr)
172 : {
173 1 : bHasDroppedRef = TRUE;
174 :
175 2 : for (int iBand = 0; iBand < nBands; iBand++)
176 : {
177 1 : delete papoBands[iBand];
178 1 : papoBands[iBand] = nullptr;
179 : }
180 1 : nBands = 0;
181 :
182 1 : GDALClose((GDALDatasetH)poDepFile);
183 1 : poDepFile = nullptr;
184 : }
185 :
186 62 : return bHasDroppedRef;
187 : }
188 :
189 : /************************************************************************/
190 : /* FlushCache() */
191 : /************************************************************************/
192 :
193 63 : CPLErr ERSDataset::FlushCache(bool bAtClosing)
194 :
195 : {
196 63 : CPLErr eErr = CE_None;
197 63 : if (bHDRDirty)
198 : {
199 37 : VSILFILE *fpERS = VSIFOpenL(GetDescription(), "w");
200 37 : if (fpERS == nullptr)
201 : {
202 0 : eErr = CE_Failure;
203 0 : CPLError(CE_Failure, CPLE_OpenFailed,
204 0 : "Unable to rewrite %s header.", GetDescription());
205 : }
206 : else
207 : {
208 37 : if (VSIFPrintfL(fpERS, "DatasetHeader Begin\n") <= 0)
209 0 : eErr = CE_Failure;
210 37 : poHeader->WriteSelf(fpERS, 1);
211 37 : if (VSIFPrintfL(fpERS, "DatasetHeader End\n") <= 0)
212 0 : eErr = CE_Failure;
213 37 : if (VSIFCloseL(fpERS) != 0)
214 0 : eErr = CE_Failure;
215 : }
216 : }
217 :
218 63 : if (RawDataset::FlushCache(bAtClosing) != CE_None)
219 0 : eErr = CE_Failure;
220 63 : return eErr;
221 : }
222 :
223 : /************************************************************************/
224 : /* GetMetadataDomainList() */
225 : /************************************************************************/
226 :
227 0 : char **ERSDataset::GetMetadataDomainList()
228 : {
229 0 : return BuildMetadataDomainList(GDALPamDataset::GetMetadataDomainList(),
230 0 : TRUE, "ERS", nullptr);
231 : }
232 :
233 : /************************************************************************/
234 : /* GetMetadataItem() */
235 : /************************************************************************/
236 :
237 60 : const char *ERSDataset::GetMetadataItem(const char *pszName,
238 : const char *pszDomain)
239 : {
240 60 : if (pszDomain != nullptr && EQUAL(pszDomain, "ERS") && pszName != nullptr)
241 : {
242 15 : if (EQUAL(pszName, "PROJ"))
243 5 : return osProj.size() ? osProj.c_str() : nullptr;
244 10 : if (EQUAL(pszName, "DATUM"))
245 5 : return osDatum.size() ? osDatum.c_str() : nullptr;
246 5 : if (EQUAL(pszName, "UNITS"))
247 5 : return osUnits.size() ? osUnits.c_str() : nullptr;
248 : }
249 45 : return GDALPamDataset::GetMetadataItem(pszName, pszDomain);
250 : }
251 :
252 : /************************************************************************/
253 : /* GetMetadata() */
254 : /************************************************************************/
255 :
256 17 : char **ERSDataset::GetMetadata(const char *pszDomain)
257 :
258 : {
259 17 : if (pszDomain != nullptr && EQUAL(pszDomain, "ERS"))
260 : {
261 1 : oERSMetadataList.Clear();
262 1 : if (!osProj.empty())
263 : oERSMetadataList.AddString(
264 1 : CPLSPrintf("%s=%s", "PROJ", osProj.c_str()));
265 1 : if (!osDatum.empty())
266 : oERSMetadataList.AddString(
267 1 : CPLSPrintf("%s=%s", "DATUM", osDatum.c_str()));
268 1 : if (!osUnits.empty())
269 : oERSMetadataList.AddString(
270 1 : CPLSPrintf("%s=%s", "UNITS", osUnits.c_str()));
271 1 : return oERSMetadataList.List();
272 : }
273 :
274 16 : return GDALPamDataset::GetMetadata(pszDomain);
275 : }
276 :
277 : /************************************************************************/
278 : /* GetGCPCount() */
279 : /************************************************************************/
280 :
281 3 : int ERSDataset::GetGCPCount()
282 :
283 : {
284 3 : return nGCPCount;
285 : }
286 :
287 : /************************************************************************/
288 : /* GetGCPSpatialRef() */
289 : /************************************************************************/
290 :
291 1 : const OGRSpatialReference *ERSDataset::GetGCPSpatialRef() const
292 :
293 : {
294 1 : return m_oGCPSRS.IsEmpty() ? nullptr : &m_oGCPSRS;
295 : }
296 :
297 : /************************************************************************/
298 : /* GetGCPs() */
299 : /************************************************************************/
300 :
301 1 : const GDAL_GCP *ERSDataset::GetGCPs()
302 :
303 : {
304 1 : return pasGCPList;
305 : }
306 :
307 : /************************************************************************/
308 : /* SetGCPs() */
309 : /************************************************************************/
310 :
311 1 : CPLErr ERSDataset::SetGCPs(int nGCPCountIn, const GDAL_GCP *pasGCPListIn,
312 : const OGRSpatialReference *poSRS)
313 :
314 : {
315 : /* -------------------------------------------------------------------- */
316 : /* Clean old gcps. */
317 : /* -------------------------------------------------------------------- */
318 1 : m_oGCPSRS.Clear();
319 :
320 1 : if (nGCPCount > 0)
321 : {
322 0 : GDALDeinitGCPs(nGCPCount, pasGCPList);
323 0 : CPLFree(pasGCPList);
324 :
325 0 : pasGCPList = nullptr;
326 0 : nGCPCount = 0;
327 : }
328 :
329 : /* -------------------------------------------------------------------- */
330 : /* Copy new ones. */
331 : /* -------------------------------------------------------------------- */
332 1 : nGCPCount = nGCPCountIn;
333 1 : pasGCPList = GDALDuplicateGCPs(nGCPCount, pasGCPListIn);
334 1 : if (poSRS)
335 1 : m_oGCPSRS = *poSRS;
336 :
337 : /* -------------------------------------------------------------------- */
338 : /* Setup the header contents corresponding to these GCPs. */
339 : /* -------------------------------------------------------------------- */
340 1 : bHDRDirty = TRUE;
341 :
342 1 : poHeader->Set("RasterInfo.WarpControl.WarpType", "Polynomial");
343 1 : if (nGCPCount > 6)
344 0 : poHeader->Set("RasterInfo.WarpControl.WarpOrder", "2");
345 : else
346 1 : poHeader->Set("RasterInfo.WarpControl.WarpOrder", "1");
347 1 : poHeader->Set("RasterInfo.WarpControl.WarpSampling", "Nearest");
348 :
349 : /* -------------------------------------------------------------------- */
350 : /* Translate the projection. */
351 : /* -------------------------------------------------------------------- */
352 : char szERSProj[32], szERSDatum[32], szERSUnits[32];
353 :
354 1 : m_oGCPSRS.exportToERM(szERSProj, szERSDatum, szERSUnits);
355 :
356 : /* Write the above computed values, unless they have been overridden by */
357 : /* the creation options PROJ, DATUM or UNITS */
358 :
359 1 : poHeader->Set("RasterInfo.WarpControl.CoordinateSpace.Datum",
360 2 : CPLString().Printf("\"%s\"", (osDatum.size())
361 0 : ? osDatum.c_str()
362 1 : : szERSDatum));
363 1 : poHeader->Set("RasterInfo.WarpControl.CoordinateSpace.Projection",
364 2 : CPLString().Printf("\"%s\"", (osProj.size()) ? osProj.c_str()
365 1 : : szERSProj));
366 1 : poHeader->Set("RasterInfo.WarpControl.CoordinateSpace.CoordinateType",
367 2 : CPLString().Printf("EN"));
368 1 : poHeader->Set("RasterInfo.WarpControl.CoordinateSpace.Units",
369 2 : CPLString().Printf("\"%s\"", (osUnits.size())
370 0 : ? osUnits.c_str()
371 1 : : szERSUnits));
372 1 : poHeader->Set("RasterInfo.WarpControl.CoordinateSpace.Rotation", "0:0:0.0");
373 :
374 : /* -------------------------------------------------------------------- */
375 : /* Translate the GCPs. */
376 : /* -------------------------------------------------------------------- */
377 1 : CPLString osControlPoints = "{\n";
378 :
379 5 : for (int iGCP = 0; iGCP < nGCPCount; iGCP++)
380 : {
381 8 : CPLString osLine;
382 :
383 8 : CPLString osId = pasGCPList[iGCP].pszId;
384 4 : if (osId.empty())
385 4 : osId.Printf("%d", iGCP + 1);
386 :
387 : osLine.Printf(
388 : "\t\t\t\t\"%s\"\tYes\tYes\t%.6f\t%.6f\t%.15g\t%.15g\t%.15g\n",
389 4 : osId.c_str(), pasGCPList[iGCP].dfGCPPixel,
390 4 : pasGCPList[iGCP].dfGCPLine, pasGCPList[iGCP].dfGCPX,
391 4 : pasGCPList[iGCP].dfGCPY, pasGCPList[iGCP].dfGCPZ);
392 4 : osControlPoints += osLine;
393 : }
394 1 : osControlPoints += "\t\t}";
395 :
396 1 : poHeader->Set("RasterInfo.WarpControl.ControlPoints", osControlPoints);
397 :
398 2 : return CE_None;
399 : }
400 :
401 : /************************************************************************/
402 : /* GetSpatialRef() */
403 : /************************************************************************/
404 :
405 3 : const OGRSpatialReference *ERSDataset::GetSpatialRef() const
406 :
407 : {
408 : // try xml first
409 3 : const auto poSRS = GDALPamDataset::GetSpatialRef();
410 3 : if (poSRS)
411 0 : return poSRS;
412 :
413 3 : return m_oSRS.IsEmpty() ? nullptr : &m_oSRS;
414 : }
415 :
416 : /************************************************************************/
417 : /* SetSpatialRef() */
418 : /************************************************************************/
419 :
420 34 : CPLErr ERSDataset::SetSpatialRef(const OGRSpatialReference *poSRS)
421 :
422 : {
423 34 : if (poSRS == nullptr && m_oSRS.IsEmpty())
424 0 : return CE_None;
425 34 : if (poSRS != nullptr && poSRS->IsSame(&m_oSRS))
426 0 : return CE_None;
427 :
428 34 : m_oSRS.Clear();
429 34 : if (poSRS)
430 34 : m_oSRS = *poSRS;
431 :
432 : char szERSProj[32], szERSDatum[32], szERSUnits[32];
433 :
434 34 : m_oSRS.exportToERM(szERSProj, szERSDatum, szERSUnits);
435 :
436 : /* Write the above computed values, unless they have been overridden by */
437 : /* the creation options PROJ, DATUM or UNITS */
438 34 : if (!osProjForced.empty())
439 1 : osProj = osProjForced;
440 : else
441 33 : osProj = szERSProj;
442 34 : if (!osDatumForced.empty())
443 1 : osDatum = osDatumForced;
444 : else
445 33 : osDatum = szERSDatum;
446 34 : if (!osUnitsForced.empty())
447 1 : osUnits = osUnitsForced;
448 : else
449 33 : osUnits = szERSUnits;
450 :
451 34 : WriteProjectionInfo(osProj, osDatum, osUnits);
452 :
453 34 : return CE_None;
454 : }
455 :
456 : /************************************************************************/
457 : /* WriteProjectionInfo() */
458 : /************************************************************************/
459 :
460 36 : void ERSDataset::WriteProjectionInfo(const char *pszProj, const char *pszDatum,
461 : const char *pszUnits)
462 : {
463 36 : bHDRDirty = TRUE;
464 36 : poHeader->Set("CoordinateSpace.Datum",
465 72 : CPLString().Printf("\"%s\"", pszDatum));
466 36 : poHeader->Set("CoordinateSpace.Projection",
467 72 : CPLString().Printf("\"%s\"", pszProj));
468 36 : poHeader->Set("CoordinateSpace.CoordinateType", CPLString().Printf("EN"));
469 36 : poHeader->Set("CoordinateSpace.Units",
470 72 : CPLString().Printf("\"%s\"", pszUnits));
471 36 : poHeader->Set("CoordinateSpace.Rotation", "0:0:0.0");
472 :
473 : /* -------------------------------------------------------------------- */
474 : /* It seems that CoordinateSpace needs to come before */
475 : /* RasterInfo. Try moving it up manually. */
476 : /* -------------------------------------------------------------------- */
477 36 : int iRasterInfo = -1;
478 36 : int iCoordSpace = -1;
479 :
480 250 : for (int i = 0; i < poHeader->nItemCount; i++)
481 : {
482 250 : if (EQUAL(poHeader->papszItemName[i], "RasterInfo"))
483 34 : iRasterInfo = i;
484 :
485 250 : if (EQUAL(poHeader->papszItemName[i], "CoordinateSpace"))
486 : {
487 36 : iCoordSpace = i;
488 36 : break;
489 : }
490 : }
491 :
492 36 : if (iCoordSpace > iRasterInfo && iRasterInfo != -1)
493 : {
494 68 : for (int i = iCoordSpace; i > 0 && i != iRasterInfo; i--)
495 : {
496 :
497 34 : ERSHdrNode *poTemp = poHeader->papoItemChild[i];
498 34 : poHeader->papoItemChild[i] = poHeader->papoItemChild[i - 1];
499 34 : poHeader->papoItemChild[i - 1] = poTemp;
500 :
501 34 : char *pszTemp = poHeader->papszItemName[i];
502 34 : poHeader->papszItemName[i] = poHeader->papszItemName[i - 1];
503 34 : poHeader->papszItemName[i - 1] = pszTemp;
504 :
505 34 : pszTemp = poHeader->papszItemValue[i];
506 34 : poHeader->papszItemValue[i] = poHeader->papszItemValue[i - 1];
507 34 : poHeader->papszItemValue[i - 1] = pszTemp;
508 : }
509 : }
510 36 : }
511 :
512 : /************************************************************************/
513 : /* GetGeoTransform() */
514 : /************************************************************************/
515 :
516 22 : CPLErr ERSDataset::GetGeoTransform(double *padfTransform)
517 :
518 : {
519 22 : if (bGotTransform)
520 : {
521 22 : memcpy(padfTransform, adfGeoTransform, sizeof(double) * 6);
522 22 : return CE_None;
523 : }
524 :
525 0 : return GDALPamDataset::GetGeoTransform(padfTransform);
526 : }
527 :
528 : /************************************************************************/
529 : /* SetGeoTransform() */
530 : /************************************************************************/
531 :
532 32 : CPLErr ERSDataset::SetGeoTransform(double *padfTransform)
533 :
534 : {
535 32 : if (memcmp(padfTransform, adfGeoTransform, sizeof(double) * 6) == 0)
536 0 : return CE_None;
537 :
538 32 : if (adfGeoTransform[2] != 0 || adfGeoTransform[4] != 0)
539 : {
540 0 : CPLError(CE_Failure, CPLE_AppDefined,
541 : "Rotated and skewed geotransforms not currently supported for "
542 : "ERS driver.");
543 0 : return CE_Failure;
544 : }
545 :
546 32 : bGotTransform = TRUE;
547 32 : memcpy(adfGeoTransform, padfTransform, sizeof(double) * 6);
548 :
549 32 : bHDRDirty = TRUE;
550 :
551 32 : poHeader->Set("RasterInfo.CellInfo.Xdimension",
552 64 : CPLString().Printf("%.15g", fabs(adfGeoTransform[1])));
553 32 : poHeader->Set("RasterInfo.CellInfo.Ydimension",
554 64 : CPLString().Printf("%.15g", fabs(adfGeoTransform[5])));
555 32 : poHeader->Set("RasterInfo.RegistrationCoord.Eastings",
556 64 : CPLString().Printf("%.15g", adfGeoTransform[0]));
557 32 : poHeader->Set("RasterInfo.RegistrationCoord.Northings",
558 64 : CPLString().Printf("%.15g", adfGeoTransform[3]));
559 :
560 64 : if (CPLAtof(poHeader->Find("RasterInfo.RegistrationCellX", "0")) != 0.0 ||
561 32 : CPLAtof(poHeader->Find("RasterInfo.RegistrationCellY", "0")) != 0.0)
562 : {
563 : // Reset RegistrationCellX/Y to 0 if the header gets rewritten (#5493)
564 0 : poHeader->Set("RasterInfo.RegistrationCellX", "0");
565 0 : poHeader->Set("RasterInfo.RegistrationCellY", "0");
566 : }
567 :
568 32 : return CE_None;
569 : }
570 :
571 : /************************************************************************/
572 : /* ERSDMS2Dec() */
573 : /* */
574 : /* Convert ERS DMS format to decimal degrees. Input is like */
575 : /* "-180:00:00". */
576 : /************************************************************************/
577 :
578 10 : static double ERSDMS2Dec(const char *pszDMS)
579 :
580 : {
581 10 : char **papszTokens = CSLTokenizeStringComplex(pszDMS, ":", FALSE, FALSE);
582 :
583 10 : if (CSLCount(papszTokens) != 3)
584 : {
585 0 : CSLDestroy(papszTokens);
586 0 : return CPLAtof(pszDMS);
587 : }
588 :
589 10 : double dfResult = fabs(CPLAtof(papszTokens[0])) +
590 10 : CPLAtof(papszTokens[1]) / 60.0 +
591 10 : CPLAtof(papszTokens[2]) / 3600.0;
592 :
593 10 : if (CPLAtof(papszTokens[0]) < 0)
594 8 : dfResult *= -1;
595 :
596 10 : CSLDestroy(papszTokens);
597 10 : return dfResult;
598 : }
599 :
600 : /************************************************************************/
601 : /* GetFileList() */
602 : /************************************************************************/
603 :
604 10 : char **ERSDataset::GetFileList()
605 :
606 : {
607 : static thread_local int nRecLevel = 0;
608 10 : if (nRecLevel > 0)
609 0 : return nullptr;
610 :
611 : // Main data file, etc.
612 10 : char **papszFileList = GDALPamDataset::GetFileList();
613 :
614 : // Add raw data file if we have one.
615 10 : if (!osRawFilename.empty())
616 10 : papszFileList = CSLAddString(papszFileList, osRawFilename);
617 :
618 : // If we have a dependent file, merge its list of files in.
619 10 : if (poDepFile)
620 : {
621 0 : nRecLevel++;
622 0 : char **papszDepFiles = poDepFile->GetFileList();
623 0 : nRecLevel--;
624 0 : papszFileList = CSLInsertStrings(papszFileList, -1, papszDepFiles);
625 0 : CSLDestroy(papszDepFiles);
626 : }
627 :
628 10 : return papszFileList;
629 : }
630 :
631 : /************************************************************************/
632 : /* ReadGCPs() */
633 : /* */
634 : /* Read the GCPs from the header. */
635 : /************************************************************************/
636 :
637 2 : void ERSDataset::ReadGCPs()
638 :
639 : {
640 : const char *pszCP =
641 2 : poHeader->Find("RasterInfo.WarpControl.ControlPoints", nullptr);
642 :
643 2 : if (pszCP == nullptr)
644 0 : return;
645 :
646 : /* -------------------------------------------------------------------- */
647 : /* Parse the control points. They will look something like: */
648 : /* */
649 : /* "1035" Yes No 2344.650885 3546.419458 483270.73 3620906.21 3.105 */
650 : /* -------------------------------------------------------------------- */
651 2 : char **papszTokens = CSLTokenizeStringComplex(pszCP, "{ \t}", TRUE, FALSE);
652 2 : int nItemCount = CSLCount(papszTokens);
653 :
654 : /* -------------------------------------------------------------------- */
655 : /* Work out if we have elevation values or not. */
656 : /* -------------------------------------------------------------------- */
657 : int nItemsPerLine;
658 :
659 2 : if (nItemCount == 7)
660 0 : nItemsPerLine = 7;
661 2 : else if (nItemCount == 8)
662 0 : nItemsPerLine = 8;
663 2 : else if (nItemCount < 14)
664 : {
665 0 : CPLDebug("ERS", "Invalid item count for ControlPoints");
666 0 : CSLDestroy(papszTokens);
667 0 : return;
668 : }
669 2 : else if (EQUAL(papszTokens[8], "Yes") || EQUAL(papszTokens[8], "No"))
670 0 : nItemsPerLine = 7;
671 2 : else if (EQUAL(papszTokens[9], "Yes") || EQUAL(papszTokens[9], "No"))
672 2 : nItemsPerLine = 8;
673 : else
674 : {
675 0 : CPLDebug("ERS", "Invalid format for ControlPoints");
676 0 : CSLDestroy(papszTokens);
677 0 : return;
678 : }
679 :
680 : /* -------------------------------------------------------------------- */
681 : /* Setup GCPs. */
682 : /* -------------------------------------------------------------------- */
683 2 : CPLAssert(nGCPCount == 0);
684 :
685 2 : nGCPCount = nItemCount / nItemsPerLine;
686 2 : pasGCPList = (GDAL_GCP *)CPLCalloc(nGCPCount, sizeof(GDAL_GCP));
687 2 : GDALInitGCPs(nGCPCount, pasGCPList);
688 :
689 10 : for (int iGCP = 0; iGCP < nGCPCount; iGCP++)
690 : {
691 8 : GDAL_GCP *psGCP = pasGCPList + iGCP;
692 :
693 8 : CPLFree(psGCP->pszId);
694 8 : psGCP->pszId = CPLStrdup(papszTokens[iGCP * nItemsPerLine + 0]);
695 8 : psGCP->dfGCPPixel = CPLAtof(papszTokens[iGCP * nItemsPerLine + 3]);
696 8 : psGCP->dfGCPLine = CPLAtof(papszTokens[iGCP * nItemsPerLine + 4]);
697 8 : psGCP->dfGCPX = CPLAtof(papszTokens[iGCP * nItemsPerLine + 5]);
698 8 : psGCP->dfGCPY = CPLAtof(papszTokens[iGCP * nItemsPerLine + 6]);
699 8 : if (nItemsPerLine == 8)
700 8 : psGCP->dfGCPZ = CPLAtof(papszTokens[iGCP * nItemsPerLine + 7]);
701 : }
702 :
703 2 : CSLDestroy(papszTokens);
704 :
705 : /* -------------------------------------------------------------------- */
706 : /* Parse the GCP projection. */
707 : /* -------------------------------------------------------------------- */
708 : osProj =
709 2 : poHeader->Find("RasterInfo.WarpControl.CoordinateSpace.Projection", "");
710 : osDatum =
711 2 : poHeader->Find("RasterInfo.WarpControl.CoordinateSpace.Datum", "");
712 : osUnits =
713 2 : poHeader->Find("RasterInfo.WarpControl.CoordinateSpace.Units", "");
714 :
715 2 : m_oGCPSRS.importFromERM(!osProj.empty() ? osProj.c_str() : "RAW",
716 2 : !osDatum.empty() ? osDatum.c_str() : "WGS84",
717 2 : !osUnits.empty() ? osUnits.c_str() : "METERS");
718 : }
719 :
720 : /************************************************************************/
721 : /* ==================================================================== */
722 : /* ERSRasterBand */
723 : /* ==================================================================== */
724 : /************************************************************************/
725 :
726 : class ERSRasterBand final : public RawRasterBand
727 : {
728 : public:
729 : ERSRasterBand(GDALDataset *poDS, int nBand, VSILFILE *fpRaw,
730 : vsi_l_offset nImgOffset, int nPixelOffset, int nLineOffset,
731 : GDALDataType eDataType, int bNativeOrder);
732 :
733 : double GetNoDataValue(int *pbSuccess = nullptr) override;
734 : CPLErr SetNoDataValue(double) override;
735 : };
736 :
737 : /************************************************************************/
738 : /* ERSRasterBand() */
739 : /************************************************************************/
740 :
741 102 : ERSRasterBand::ERSRasterBand(GDALDataset *poDSIn, int nBandIn,
742 : VSILFILE *fpRawIn, vsi_l_offset nImgOffsetIn,
743 : int nPixelOffsetIn, int nLineOffsetIn,
744 102 : GDALDataType eDataTypeIn, int bNativeOrderIn)
745 : : RawRasterBand(poDSIn, nBandIn, fpRawIn, nImgOffsetIn, nPixelOffsetIn,
746 : nLineOffsetIn, eDataTypeIn, bNativeOrderIn,
747 102 : RawRasterBand::OwnFP::NO)
748 : {
749 102 : }
750 :
751 : /************************************************************************/
752 : /* GetNoDataValue() */
753 : /************************************************************************/
754 :
755 15 : double ERSRasterBand::GetNoDataValue(int *pbSuccess)
756 : {
757 15 : ERSDataset *poGDS = cpl::down_cast<ERSDataset *>(poDS);
758 15 : if (poGDS->bHasNoDataValue)
759 : {
760 1 : if (pbSuccess)
761 1 : *pbSuccess = TRUE;
762 1 : return poGDS->dfNoDataValue;
763 : }
764 :
765 14 : return RawRasterBand::GetNoDataValue(pbSuccess);
766 : }
767 :
768 : /************************************************************************/
769 : /* SetNoDataValue() */
770 : /************************************************************************/
771 :
772 1 : CPLErr ERSRasterBand::SetNoDataValue(double dfNoDataValue)
773 : {
774 1 : ERSDataset *poGDS = cpl::down_cast<ERSDataset *>(poDS);
775 1 : if (!poGDS->bHasNoDataValue || poGDS->dfNoDataValue != dfNoDataValue)
776 : {
777 1 : poGDS->bHasNoDataValue = TRUE;
778 1 : poGDS->dfNoDataValue = dfNoDataValue;
779 :
780 1 : poGDS->bHDRDirty = TRUE;
781 1 : poGDS->poHeader->Set("RasterInfo.NullCellValue",
782 2 : CPLString().Printf("%.16g", dfNoDataValue));
783 : }
784 1 : return CE_None;
785 : }
786 :
787 : /************************************************************************/
788 : /* Identify() */
789 : /************************************************************************/
790 :
791 53931 : int ERSDataset::Identify(GDALOpenInfo *poOpenInfo)
792 :
793 : {
794 : /* -------------------------------------------------------------------- */
795 : /* We assume the user selects the .ers file. */
796 : /* -------------------------------------------------------------------- */
797 53931 : CPLString osHeader((const char *)poOpenInfo->pabyHeader,
798 107856 : poOpenInfo->nHeaderBytes);
799 :
800 53928 : if (osHeader.ifind("Algorithm Begin") != std::string::npos)
801 : {
802 0 : CPLError(CE_Failure, CPLE_OpenFailed,
803 : "%s appears to be an algorithm ERS file, which is not "
804 : "currently supported.",
805 : poOpenInfo->pszFilename);
806 0 : return FALSE;
807 : }
808 :
809 53924 : if (osHeader.ifind("DatasetHeader ") != std::string::npos)
810 89 : return TRUE;
811 :
812 53836 : return FALSE;
813 : }
814 :
815 : /************************************************************************/
816 : /* ERSProxyRasterBand */
817 : /************************************************************************/
818 :
819 : namespace
820 : {
821 :
822 63 : static int &GetRecLevel()
823 : {
824 : static thread_local int nRecLevel = 0;
825 63 : return nRecLevel;
826 : }
827 :
828 : class ERSProxyRasterBand final : public GDALProxyRasterBand
829 : {
830 : public:
831 1 : explicit ERSProxyRasterBand(GDALRasterBand *poUnderlyingBand)
832 1 : : m_poUnderlyingBand(poUnderlyingBand)
833 : {
834 1 : poUnderlyingBand->GetBlockSize(&nBlockXSize, &nBlockYSize);
835 1 : eDataType = poUnderlyingBand->GetRasterDataType();
836 1 : }
837 :
838 : int GetOverviewCount() override;
839 :
840 : protected:
841 1 : GDALRasterBand *RefUnderlyingRasterBand(bool /*bForceOpen*/) const override
842 : {
843 1 : return m_poUnderlyingBand;
844 : }
845 :
846 : private:
847 : GDALRasterBand *m_poUnderlyingBand;
848 : };
849 :
850 0 : int ERSProxyRasterBand::GetOverviewCount()
851 : {
852 0 : int &nRecLevel = GetRecLevel();
853 0 : nRecLevel++;
854 0 : int nRet = GDALProxyRasterBand::GetOverviewCount();
855 0 : nRecLevel--;
856 0 : return nRet;
857 : }
858 :
859 : } // namespace
860 :
861 : /************************************************************************/
862 : /* Open() */
863 : /************************************************************************/
864 :
865 63 : GDALDataset *ERSDataset::Open(GDALOpenInfo *poOpenInfo)
866 :
867 : {
868 63 : if (!Identify(poOpenInfo) || poOpenInfo->fpL == nullptr)
869 0 : return nullptr;
870 :
871 63 : int &nRecLevel = GetRecLevel();
872 : // cppcheck-suppress knownConditionTrueFalse
873 63 : if (nRecLevel)
874 : {
875 1 : CPLError(CE_Failure, CPLE_AppDefined,
876 : "Attempt at recursively opening ERS dataset");
877 1 : return nullptr;
878 : }
879 :
880 : /* -------------------------------------------------------------------- */
881 : /* Ingest the file as a tree of header nodes. */
882 : /* -------------------------------------------------------------------- */
883 62 : ERSHdrNode *poHeader = new ERSHdrNode();
884 :
885 62 : if (!poHeader->ParseHeader(poOpenInfo->fpL))
886 : {
887 0 : delete poHeader;
888 0 : VSIFCloseL(poOpenInfo->fpL);
889 0 : poOpenInfo->fpL = nullptr;
890 0 : return nullptr;
891 : }
892 :
893 62 : VSIFCloseL(poOpenInfo->fpL);
894 62 : poOpenInfo->fpL = nullptr;
895 :
896 : /* -------------------------------------------------------------------- */
897 : /* Do we have the minimum required information from this header? */
898 : /* -------------------------------------------------------------------- */
899 62 : if (poHeader->Find("RasterInfo.NrOfLines") == nullptr ||
900 124 : poHeader->Find("RasterInfo.NrOfCellsPerLine") == nullptr ||
901 62 : poHeader->Find("RasterInfo.NrOfBands") == nullptr)
902 : {
903 0 : if (poHeader->FindNode("Algorithm") != nullptr)
904 : {
905 0 : CPLError(CE_Failure, CPLE_OpenFailed,
906 : "%s appears to be an algorithm ERS file, which is not "
907 : "currently supported.",
908 : poOpenInfo->pszFilename);
909 : }
910 0 : delete poHeader;
911 0 : return nullptr;
912 : }
913 :
914 : /* -------------------------------------------------------------------- */
915 : /* Create a corresponding GDALDataset. */
916 : /* -------------------------------------------------------------------- */
917 124 : auto poDS = std::make_unique<ERSDataset>();
918 62 : poDS->poHeader = poHeader;
919 62 : poDS->eAccess = poOpenInfo->eAccess;
920 :
921 : /* -------------------------------------------------------------------- */
922 : /* Capture some information from the file that is of interest. */
923 : /* -------------------------------------------------------------------- */
924 62 : int nBands = atoi(poHeader->Find("RasterInfo.NrOfBands"));
925 62 : poDS->nRasterXSize = atoi(poHeader->Find("RasterInfo.NrOfCellsPerLine"));
926 62 : poDS->nRasterYSize = atoi(poHeader->Find("RasterInfo.NrOfLines"));
927 :
928 124 : if (!GDALCheckDatasetDimensions(poDS->nRasterXSize, poDS->nRasterYSize) ||
929 62 : !GDALCheckBandCount(nBands, FALSE))
930 : {
931 0 : return nullptr;
932 : }
933 :
934 : /* -------------------------------------------------------------------- */
935 : /* Get the HeaderOffset if it exists in the header */
936 : /* -------------------------------------------------------------------- */
937 62 : GIntBig nHeaderOffset = 0;
938 62 : const char *pszHeaderOffset = poHeader->Find("HeaderOffset");
939 62 : if (pszHeaderOffset != nullptr)
940 : {
941 2 : nHeaderOffset = CPLAtoGIntBig(pszHeaderOffset);
942 2 : if (nHeaderOffset < 0)
943 : {
944 0 : CPLError(CE_Failure, CPLE_AppDefined,
945 : "Illegal value for HeaderOffset: %s", pszHeaderOffset);
946 0 : return nullptr;
947 : }
948 : }
949 :
950 : /* -------------------------------------------------------------------- */
951 : /* Establish the data type. */
952 : /* -------------------------------------------------------------------- */
953 : CPLString osCellType =
954 124 : poHeader->Find("RasterInfo.CellType", "Unsigned8BitInteger");
955 : GDALDataType eType;
956 62 : if (EQUAL(osCellType, "Unsigned8BitInteger"))
957 29 : eType = GDT_Byte;
958 33 : else if (EQUAL(osCellType, "Signed8BitInteger"))
959 6 : eType = GDT_Int8;
960 27 : else if (EQUAL(osCellType, "Unsigned16BitInteger"))
961 3 : eType = GDT_UInt16;
962 24 : else if (EQUAL(osCellType, "Signed16BitInteger"))
963 6 : eType = GDT_Int16;
964 18 : else if (EQUAL(osCellType, "Unsigned32BitInteger"))
965 3 : eType = GDT_UInt32;
966 15 : else if (EQUAL(osCellType, "Signed32BitInteger"))
967 3 : eType = GDT_Int32;
968 12 : else if (EQUAL(osCellType, "IEEE4ByteReal"))
969 9 : eType = GDT_Float32;
970 3 : else if (EQUAL(osCellType, "IEEE8ByteReal"))
971 3 : eType = GDT_Float64;
972 : else
973 : {
974 0 : CPLDebug("ERS", "Unknown CellType '%s'", osCellType.c_str());
975 0 : eType = GDT_Byte;
976 : }
977 :
978 : /* -------------------------------------------------------------------- */
979 : /* Pick up the word order. */
980 : /* -------------------------------------------------------------------- */
981 : const int bNative =
982 : #ifdef CPL_LSB
983 62 : EQUAL(poHeader->Find("ByteOrder", "LSBFirst"), "LSBFirst")
984 : #else
985 : EQUAL(poHeader->Find("ByteOrder", "MSBFirst"), "MSBFirst")
986 : #endif
987 : ;
988 : /* -------------------------------------------------------------------- */
989 : /* Figure out the name of the target file. */
990 : /* -------------------------------------------------------------------- */
991 124 : CPLString osPath = CPLGetPath(poOpenInfo->pszFilename);
992 124 : CPLString osDataFile = poHeader->Find("DataFile", "");
993 :
994 62 : if (osDataFile.length() == 0) // just strip off extension.
995 : {
996 61 : osDataFile = CPLGetFilename(poOpenInfo->pszFilename);
997 61 : osDataFile = osDataFile.substr(0, osDataFile.find_last_of('.'));
998 : }
999 :
1000 124 : CPLString osDataFilePath = CPLFormFilename(osPath, osDataFile, nullptr);
1001 :
1002 : /* -------------------------------------------------------------------- */
1003 : /* DataSetType = Translated files are links to things like ecw */
1004 : /* files. */
1005 : /* -------------------------------------------------------------------- */
1006 62 : if (EQUAL(poHeader->Find("DataSetType", ""), "Translated"))
1007 : {
1008 2 : nRecLevel++;
1009 2 : poDS->poDepFile = GDALDataset::FromHandle(
1010 : GDALOpen(osDataFilePath, poOpenInfo->eAccess));
1011 2 : nRecLevel--;
1012 :
1013 2 : if (poDS->poDepFile != nullptr &&
1014 1 : poDS->poDepFile->GetRasterXSize() == poDS->GetRasterXSize() &&
1015 4 : poDS->poDepFile->GetRasterYSize() == poDS->GetRasterYSize() &&
1016 1 : poDS->poDepFile->GetRasterCount() >= nBands)
1017 : {
1018 2 : for (int iBand = 0; iBand < nBands; iBand++)
1019 : {
1020 : // Assume pixel interleaved.
1021 2 : poDS->SetBand(iBand + 1,
1022 : new ERSProxyRasterBand(
1023 1 : poDS->poDepFile->GetRasterBand(iBand + 1)));
1024 : }
1025 : }
1026 : else
1027 : {
1028 1 : delete poDS->poDepFile;
1029 1 : poDS->poDepFile = nullptr;
1030 : }
1031 : }
1032 :
1033 : /* ==================================================================== */
1034 : /* While ERStorage indicates a raw file. */
1035 : /* ==================================================================== */
1036 60 : else if (EQUAL(poHeader->Find("DataSetType", ""), "ERStorage"))
1037 : {
1038 : // Open data file.
1039 60 : if (poOpenInfo->eAccess == GA_Update)
1040 38 : poDS->fpImage = VSIFOpenL(osDataFilePath, "r+");
1041 : else
1042 22 : poDS->fpImage = VSIFOpenL(osDataFilePath, "r");
1043 :
1044 60 : poDS->osRawFilename = std::move(osDataFilePath);
1045 :
1046 60 : if (poDS->fpImage != nullptr && nBands > 0)
1047 : {
1048 60 : int iWordSize = GDALGetDataTypeSizeBytes(eType);
1049 :
1050 60 : const auto knIntMax = std::numeric_limits<int>::max();
1051 120 : if (nBands > knIntMax / iWordSize ||
1052 60 : poDS->nRasterXSize > knIntMax / (nBands * iWordSize))
1053 : {
1054 0 : CPLError(CE_Failure, CPLE_AppDefined,
1055 : "int overflow: too large nBands and/or nRasterXSize");
1056 0 : return nullptr;
1057 : }
1058 :
1059 60 : if (!RAWDatasetCheckMemoryUsage(
1060 60 : poDS->nRasterXSize, poDS->nRasterYSize, nBands, iWordSize,
1061 60 : iWordSize, iWordSize * nBands * poDS->nRasterXSize,
1062 : nHeaderOffset,
1063 60 : static_cast<vsi_l_offset>(iWordSize) * poDS->nRasterXSize,
1064 60 : poDS->fpImage))
1065 : {
1066 0 : return nullptr;
1067 : }
1068 60 : if (nHeaderOffset > std::numeric_limits<GIntBig>::max() -
1069 60 : static_cast<GIntBig>(nBands - 1) *
1070 60 : iWordSize * poDS->nRasterXSize)
1071 : {
1072 0 : CPLError(CE_Failure, CPLE_AppDefined,
1073 : "int overflow: too large nHeaderOffset");
1074 0 : return nullptr;
1075 : }
1076 :
1077 162 : for (int iBand = 0; iBand < nBands; iBand++)
1078 : {
1079 : // Assume pixel interleaved.
1080 : auto poBand = std::make_unique<ERSRasterBand>(
1081 102 : poDS.get(), iBand + 1, poDS->fpImage,
1082 0 : nHeaderOffset + static_cast<vsi_l_offset>(iWordSize) *
1083 102 : iBand * poDS->nRasterXSize,
1084 102 : iWordSize, iWordSize * nBands * poDS->nRasterXSize, eType,
1085 102 : bNative);
1086 102 : if (!poBand->IsValid())
1087 0 : return nullptr;
1088 102 : poDS->SetBand(iBand + 1, std::move(poBand));
1089 : }
1090 : }
1091 : }
1092 :
1093 : /* -------------------------------------------------------------------- */
1094 : /* Otherwise we have an error! */
1095 : /* -------------------------------------------------------------------- */
1096 62 : if (poDS->nBands == 0)
1097 : {
1098 1 : return nullptr;
1099 : }
1100 :
1101 : /* -------------------------------------------------------------------- */
1102 : /* Look for band descriptions. */
1103 : /* -------------------------------------------------------------------- */
1104 61 : ERSHdrNode *poRI = poHeader->FindNode("RasterInfo");
1105 :
1106 357 : for (int iChild = 0, iBand = 0;
1107 357 : poRI != nullptr && iChild < poRI->nItemCount && iBand < poDS->nBands;
1108 : iChild++)
1109 : {
1110 296 : if (poRI->papoItemChild[iChild] != nullptr &&
1111 33 : EQUAL(poRI->papszItemName[iChild], "BandId"))
1112 : {
1113 : const char *pszValue =
1114 12 : poRI->papoItemChild[iChild]->Find("Value", nullptr);
1115 :
1116 12 : iBand++;
1117 12 : if (pszValue)
1118 : {
1119 12 : CPLPushErrorHandler(CPLQuietErrorHandler);
1120 12 : poDS->GetRasterBand(iBand)->SetDescription(pszValue);
1121 12 : CPLPopErrorHandler();
1122 : }
1123 :
1124 12 : pszValue = poRI->papoItemChild[iChild]->Find("Units", nullptr);
1125 12 : if (pszValue)
1126 : {
1127 4 : CPLPushErrorHandler(CPLQuietErrorHandler);
1128 4 : poDS->GetRasterBand(iBand)->SetUnitType(pszValue);
1129 4 : CPLPopErrorHandler();
1130 : }
1131 : }
1132 : }
1133 :
1134 : /* -------------------------------------------------------------------- */
1135 : /* Look for projection. */
1136 : /* -------------------------------------------------------------------- */
1137 61 : poDS->osProj = poHeader->Find("CoordinateSpace.Projection", "");
1138 61 : poDS->osDatum = poHeader->Find("CoordinateSpace.Datum", "");
1139 61 : poDS->osUnits = poHeader->Find("CoordinateSpace.Units", "");
1140 :
1141 219 : poDS->m_oSRS.importFromERM(
1142 61 : !poDS->osProj.empty() ? poDS->osProj.c_str() : "RAW",
1143 61 : !poDS->osDatum.empty() ? poDS->osDatum.c_str() : "WGS84",
1144 61 : !poDS->osUnits.empty() ? poDS->osUnits.c_str() : "METERS");
1145 :
1146 : /* -------------------------------------------------------------------- */
1147 : /* Look for the geotransform. */
1148 : /* -------------------------------------------------------------------- */
1149 61 : if (poHeader->Find("RasterInfo.RegistrationCoord.Eastings", nullptr))
1150 : {
1151 6 : poDS->bGotTransform = TRUE;
1152 6 : poDS->adfGeoTransform[0] = CPLAtof(
1153 : poHeader->Find("RasterInfo.RegistrationCoord.Eastings", ""));
1154 12 : poDS->adfGeoTransform[1] =
1155 6 : CPLAtof(poHeader->Find("RasterInfo.CellInfo.Xdimension", "1.0"));
1156 6 : poDS->adfGeoTransform[2] = 0.0;
1157 6 : poDS->adfGeoTransform[3] = CPLAtof(
1158 : poHeader->Find("RasterInfo.RegistrationCoord.Northings", ""));
1159 6 : poDS->adfGeoTransform[4] = 0.0;
1160 6 : poDS->adfGeoTransform[5] =
1161 6 : -CPLAtof(poHeader->Find("RasterInfo.CellInfo.Ydimension", "1.0"));
1162 : }
1163 60 : else if (poHeader->Find("RasterInfo.RegistrationCoord.Latitude", nullptr) &&
1164 5 : poHeader->Find("RasterInfo.CellInfo.Xdimension", nullptr))
1165 : {
1166 5 : poDS->bGotTransform = TRUE;
1167 5 : poDS->adfGeoTransform[0] = ERSDMS2Dec(
1168 : poHeader->Find("RasterInfo.RegistrationCoord.Longitude", ""));
1169 10 : poDS->adfGeoTransform[1] =
1170 5 : CPLAtof(poHeader->Find("RasterInfo.CellInfo.Xdimension", ""));
1171 5 : poDS->adfGeoTransform[2] = 0.0;
1172 5 : poDS->adfGeoTransform[3] = ERSDMS2Dec(
1173 : poHeader->Find("RasterInfo.RegistrationCoord.Latitude", ""));
1174 5 : poDS->adfGeoTransform[4] = 0.0;
1175 5 : poDS->adfGeoTransform[5] =
1176 5 : -CPLAtof(poHeader->Find("RasterInfo.CellInfo.Ydimension", ""));
1177 : }
1178 :
1179 : /* -------------------------------------------------------------------- */
1180 : /* Adjust if we have a registration cell. */
1181 : /* -------------------------------------------------------------------- */
1182 :
1183 : /* http://geospatial.intergraph.com/Libraries/Tech_Docs/ERDAS_ER_Mapper_Customization_Guide.sflb.ashx
1184 : */
1185 : /* Page 27 : */
1186 : /* RegistrationCellX and RegistrationCellY : The image X and Y
1187 : coordinates of the cell which corresponds to the Registration
1188 : Coordinate. Note that the RegistrationCellX and
1189 : RegistrationCellY can be fractional values. If
1190 : RegistrationCellX and RegistrationCellY are not specified,
1191 : they are assumed to be (0,0), which is the top left corner of the
1192 : image.
1193 : */
1194 : double dfCellX =
1195 61 : CPLAtof(poHeader->Find("RasterInfo.RegistrationCellX", "0"));
1196 : double dfCellY =
1197 61 : CPLAtof(poHeader->Find("RasterInfo.RegistrationCellY", "0"));
1198 :
1199 61 : if (poDS->bGotTransform)
1200 : {
1201 11 : poDS->adfGeoTransform[0] -= dfCellX * poDS->adfGeoTransform[1] +
1202 11 : dfCellY * poDS->adfGeoTransform[2];
1203 22 : poDS->adfGeoTransform[3] -= dfCellX * poDS->adfGeoTransform[4] +
1204 11 : dfCellY * poDS->adfGeoTransform[5];
1205 : }
1206 :
1207 : /* -------------------------------------------------------------------- */
1208 : /* Check for null values. */
1209 : /* -------------------------------------------------------------------- */
1210 61 : if (poHeader->Find("RasterInfo.NullCellValue", nullptr))
1211 : {
1212 8 : poDS->bHasNoDataValue = TRUE;
1213 16 : poDS->dfNoDataValue =
1214 8 : CPLAtofM(poHeader->Find("RasterInfo.NullCellValue"));
1215 :
1216 8 : if (poDS->poDepFile != nullptr)
1217 : {
1218 0 : CPLPushErrorHandler(CPLQuietErrorHandler);
1219 :
1220 0 : for (int iBand = 1; iBand <= poDS->nBands; iBand++)
1221 0 : poDS->GetRasterBand(iBand)->SetNoDataValue(poDS->dfNoDataValue);
1222 :
1223 0 : CPLPopErrorHandler();
1224 : }
1225 : }
1226 :
1227 : /* -------------------------------------------------------------------- */
1228 : /* Do we have an "All" region? */
1229 : /* -------------------------------------------------------------------- */
1230 61 : ERSHdrNode *poAll = nullptr;
1231 :
1232 366 : for (int iChild = 0; poRI != nullptr && iChild < poRI->nItemCount; iChild++)
1233 : {
1234 305 : if (poRI->papoItemChild[iChild] != nullptr &&
1235 39 : EQUAL(poRI->papszItemName[iChild], "RegionInfo"))
1236 : {
1237 3 : if (EQUAL(poRI->papoItemChild[iChild]->Find("RegionName", ""),
1238 : "All"))
1239 3 : poAll = poRI->papoItemChild[iChild];
1240 : }
1241 : }
1242 :
1243 : /* -------------------------------------------------------------------- */
1244 : /* Do we have statistics? */
1245 : /* -------------------------------------------------------------------- */
1246 61 : if (poAll && poAll->FindNode("Stats"))
1247 : {
1248 3 : CPLPushErrorHandler(CPLQuietErrorHandler);
1249 :
1250 6 : for (int iBand = 1; iBand <= poDS->nBands; iBand++)
1251 : {
1252 : const char *pszValue =
1253 3 : poAll->FindElem("Stats.MinimumValue", iBand - 1);
1254 :
1255 3 : if (pszValue)
1256 3 : poDS->GetRasterBand(iBand)->SetMetadataItem(
1257 3 : "STATISTICS_MINIMUM", pszValue);
1258 :
1259 3 : pszValue = poAll->FindElem("Stats.MaximumValue", iBand - 1);
1260 :
1261 3 : if (pszValue)
1262 3 : poDS->GetRasterBand(iBand)->SetMetadataItem(
1263 3 : "STATISTICS_MAXIMUM", pszValue);
1264 :
1265 3 : pszValue = poAll->FindElem("Stats.MeanValue", iBand - 1);
1266 :
1267 3 : if (pszValue)
1268 3 : poDS->GetRasterBand(iBand)->SetMetadataItem("STATISTICS_MEAN",
1269 3 : pszValue);
1270 :
1271 3 : pszValue = poAll->FindElem("Stats.MedianValue", iBand - 1);
1272 :
1273 3 : if (pszValue)
1274 3 : poDS->GetRasterBand(iBand)->SetMetadataItem("STATISTICS_MEDIAN",
1275 3 : pszValue);
1276 : }
1277 :
1278 3 : CPLPopErrorHandler();
1279 : }
1280 :
1281 : /* -------------------------------------------------------------------- */
1282 : /* Do we have GCPs. */
1283 : /* -------------------------------------------------------------------- */
1284 61 : if (poHeader->FindNode("RasterInfo.WarpControl"))
1285 2 : poDS->ReadGCPs();
1286 :
1287 : /* -------------------------------------------------------------------- */
1288 : /* Initialize any PAM information. */
1289 : /* -------------------------------------------------------------------- */
1290 61 : poDS->SetDescription(poOpenInfo->pszFilename);
1291 61 : poDS->TryLoadXML();
1292 :
1293 : // if no SR in xml, try aux
1294 61 : const OGRSpatialReference *poSRS = poDS->GDALPamDataset::GetSpatialRef();
1295 61 : if (poSRS == nullptr)
1296 : {
1297 : // try aux
1298 : auto poAuxDS = std::unique_ptr<GDALDataset>(GDALFindAssociatedAuxFile(
1299 122 : poOpenInfo->pszFilename, GA_ReadOnly, poDS.get()));
1300 61 : if (poAuxDS)
1301 : {
1302 0 : poSRS = poAuxDS->GetSpatialRef();
1303 0 : if (poSRS)
1304 : {
1305 0 : poDS->m_oSRS = *poSRS;
1306 : }
1307 : }
1308 : }
1309 : /* -------------------------------------------------------------------- */
1310 : /* Check for overviews. */
1311 : /* -------------------------------------------------------------------- */
1312 61 : poDS->oOvManager.Initialize(poDS.get(), poOpenInfo->pszFilename);
1313 :
1314 61 : return poDS.release();
1315 : }
1316 :
1317 : /************************************************************************/
1318 : /* Create() */
1319 : /************************************************************************/
1320 :
1321 67 : GDALDataset *ERSDataset::Create(const char *pszFilename, int nXSize, int nYSize,
1322 : int nBandsIn, GDALDataType eType,
1323 : char **papszOptions)
1324 :
1325 : {
1326 : /* -------------------------------------------------------------------- */
1327 : /* Verify settings. */
1328 : /* -------------------------------------------------------------------- */
1329 67 : if (nBandsIn <= 0)
1330 : {
1331 1 : CPLError(CE_Failure, CPLE_NotSupported,
1332 : "ERS driver does not support %d bands.\n", nBandsIn);
1333 1 : return nullptr;
1334 : }
1335 :
1336 66 : if (eType != GDT_Byte && eType != GDT_Int8 && eType != GDT_Int16 &&
1337 29 : eType != GDT_UInt16 && eType != GDT_Int32 && eType != GDT_UInt32 &&
1338 19 : eType != GDT_Float32 && eType != GDT_Float64)
1339 : {
1340 16 : CPLError(
1341 : CE_Failure, CPLE_AppDefined,
1342 : "The ERS driver does not supporting creating files of types %s.",
1343 : GDALGetDataTypeName(eType));
1344 16 : return nullptr;
1345 : }
1346 :
1347 : /* -------------------------------------------------------------------- */
1348 : /* Work out the name we want to use for the .ers and binary */
1349 : /* data files. */
1350 : /* -------------------------------------------------------------------- */
1351 100 : CPLString osBinFile, osErsFile;
1352 :
1353 50 : if (EQUAL(CPLGetExtension(pszFilename), "ers"))
1354 : {
1355 7 : osErsFile = pszFilename;
1356 7 : osBinFile = osErsFile.substr(0, osErsFile.length() - 4);
1357 : }
1358 : else
1359 : {
1360 43 : osBinFile = pszFilename;
1361 43 : osErsFile = osBinFile + ".ers";
1362 : }
1363 :
1364 : /* -------------------------------------------------------------------- */
1365 : /* Work out some values we will write. */
1366 : /* -------------------------------------------------------------------- */
1367 50 : const char *pszCellType = "Unsigned8BitInteger";
1368 :
1369 50 : if (eType == GDT_Byte)
1370 28 : pszCellType = "Unsigned8BitInteger";
1371 22 : else if (eType == GDT_Int8)
1372 3 : pszCellType = "Signed8BitInteger";
1373 19 : else if (eType == GDT_Int16)
1374 3 : pszCellType = "Signed16BitInteger";
1375 16 : else if (eType == GDT_UInt16)
1376 3 : pszCellType = "Unsigned16BitInteger";
1377 13 : else if (eType == GDT_Int32)
1378 3 : pszCellType = "Signed32BitInteger";
1379 10 : else if (eType == GDT_UInt32)
1380 3 : pszCellType = "Unsigned32BitInteger";
1381 7 : else if (eType == GDT_Float32)
1382 4 : pszCellType = "IEEE4ByteReal";
1383 3 : else if (eType == GDT_Float64)
1384 3 : pszCellType = "IEEE8ByteReal";
1385 : else
1386 : {
1387 0 : CPLAssert(false);
1388 : }
1389 :
1390 : /* -------------------------------------------------------------------- */
1391 : /* Handling for signed eight bit data. */
1392 : /* -------------------------------------------------------------------- */
1393 50 : const char *pszPixelType = CSLFetchNameValue(papszOptions, "PIXELTYPE");
1394 50 : if (pszPixelType && EQUAL(pszPixelType, "SIGNEDBYTE") && eType == GDT_Byte)
1395 0 : pszCellType = "Signed8BitInteger";
1396 :
1397 : /* -------------------------------------------------------------------- */
1398 : /* Write binary file. */
1399 : /* -------------------------------------------------------------------- */
1400 50 : VSILFILE *fpBin = VSIFOpenL(osBinFile, "w");
1401 :
1402 50 : if (fpBin == nullptr)
1403 : {
1404 3 : CPLError(CE_Failure, CPLE_FileIO, "Failed to create %s:\n%s",
1405 3 : osBinFile.c_str(), VSIStrerror(errno));
1406 3 : return nullptr;
1407 : }
1408 :
1409 : GUIntBig nSize =
1410 47 : nXSize * (GUIntBig)nYSize * nBandsIn * (GDALGetDataTypeSize(eType) / 8);
1411 47 : GByte byZero = 0;
1412 94 : if (VSIFSeekL(fpBin, nSize - 1, SEEK_SET) != 0 ||
1413 47 : VSIFWriteL(&byZero, 1, 1, fpBin) != 1)
1414 : {
1415 10 : CPLError(CE_Failure, CPLE_FileIO, "Failed to write %s:\n%s",
1416 10 : osBinFile.c_str(), VSIStrerror(errno));
1417 10 : VSIFCloseL(fpBin);
1418 10 : return nullptr;
1419 : }
1420 37 : VSIFCloseL(fpBin);
1421 :
1422 : /* -------------------------------------------------------------------- */
1423 : /* Try writing header file. */
1424 : /* -------------------------------------------------------------------- */
1425 37 : VSILFILE *fpERS = VSIFOpenL(osErsFile, "w");
1426 :
1427 37 : if (fpERS == nullptr)
1428 : {
1429 0 : CPLError(CE_Failure, CPLE_FileIO, "Failed to create %s:\n%s",
1430 0 : osErsFile.c_str(), VSIStrerror(errno));
1431 0 : return nullptr;
1432 : }
1433 :
1434 37 : VSIFPrintfL(fpERS, "DatasetHeader Begin\n");
1435 37 : VSIFPrintfL(fpERS, "\tVersion\t\t = \"6.0\"\n");
1436 37 : VSIFPrintfL(fpERS, "\tName\t\t= \"%s\"\n", CPLGetFilename(osErsFile));
1437 :
1438 : // Last updated requires timezone info which we don't necessarily get
1439 : // get from VSICTime() so perhaps it is better to omit this.
1440 : // VSIFPrintfL( fpERS, "\tLastUpdated\t= %s",
1441 : // VSICTime( VSITime( NULL ) ) );
1442 :
1443 37 : VSIFPrintfL(fpERS, "\tDataSetType\t= ERStorage\n");
1444 37 : VSIFPrintfL(fpERS, "\tDataType\t= Raster\n");
1445 37 : VSIFPrintfL(fpERS, "\tByteOrder\t= LSBFirst\n");
1446 37 : VSIFPrintfL(fpERS, "\tRasterInfo Begin\n");
1447 37 : VSIFPrintfL(fpERS, "\t\tCellType\t= %s\n", pszCellType);
1448 37 : VSIFPrintfL(fpERS, "\t\tNrOfLines\t= %d\n", nYSize);
1449 37 : VSIFPrintfL(fpERS, "\t\tNrOfCellsPerLine\t= %d\n", nXSize);
1450 37 : VSIFPrintfL(fpERS, "\t\tNrOfBands\t= %d\n", nBandsIn);
1451 37 : VSIFPrintfL(fpERS, "\tRasterInfo End\n");
1452 37 : if (VSIFPrintfL(fpERS, "DatasetHeader End\n") < 17)
1453 : {
1454 0 : CPLError(CE_Failure, CPLE_FileIO, "Failed to write %s:\n%s",
1455 0 : osErsFile.c_str(), VSIStrerror(errno));
1456 0 : return nullptr;
1457 : }
1458 :
1459 37 : VSIFCloseL(fpERS);
1460 :
1461 : /* -------------------------------------------------------------------- */
1462 : /* Reopen. */
1463 : /* -------------------------------------------------------------------- */
1464 74 : GDALOpenInfo oOpenInfo(osErsFile, GA_Update);
1465 37 : ERSDataset *poDS = (ERSDataset *)Open(&oOpenInfo);
1466 37 : if (poDS == nullptr)
1467 0 : return nullptr;
1468 :
1469 : /* -------------------------------------------------------------------- */
1470 : /* Fetch DATUM, PROJ and UNITS creation option */
1471 : /* -------------------------------------------------------------------- */
1472 37 : const char *pszDatum = CSLFetchNameValue(papszOptions, "DATUM");
1473 37 : if (pszDatum)
1474 : {
1475 2 : poDS->osDatumForced = pszDatum;
1476 2 : poDS->osDatum = pszDatum;
1477 : }
1478 37 : const char *pszProj = CSLFetchNameValue(papszOptions, "PROJ");
1479 37 : if (pszProj)
1480 : {
1481 2 : poDS->osProjForced = pszProj;
1482 2 : poDS->osProj = pszProj;
1483 : }
1484 37 : const char *pszUnits = CSLFetchNameValue(papszOptions, "UNITS");
1485 37 : if (pszUnits)
1486 : {
1487 2 : poDS->osUnitsForced = pszUnits;
1488 2 : poDS->osUnits = pszUnits;
1489 : }
1490 :
1491 37 : if (pszDatum || pszProj || pszUnits)
1492 : {
1493 2 : poDS->WriteProjectionInfo(pszProj ? pszProj : "RAW",
1494 : pszDatum ? pszDatum : "RAW",
1495 : pszUnits ? pszUnits : "METERS");
1496 : }
1497 :
1498 37 : return poDS;
1499 : }
1500 :
1501 : /************************************************************************/
1502 : /* GDALRegister_ERS() */
1503 : /************************************************************************/
1504 :
1505 1595 : void GDALRegister_ERS()
1506 :
1507 : {
1508 1595 : if (GDALGetDriverByName("ERS") != nullptr)
1509 302 : return;
1510 :
1511 1293 : GDALDriver *poDriver = new GDALDriver();
1512 :
1513 1293 : poDriver->SetDescription("ERS");
1514 1293 : poDriver->SetMetadataItem(GDAL_DCAP_RASTER, "YES");
1515 1293 : poDriver->SetMetadataItem(GDAL_DMD_LONGNAME, "ERMapper .ers Labelled");
1516 1293 : poDriver->SetMetadataItem(GDAL_DMD_HELPTOPIC, "drivers/raster/ers.html");
1517 1293 : poDriver->SetMetadataItem(GDAL_DMD_EXTENSION, "ers");
1518 1293 : poDriver->SetMetadataItem(GDAL_DMD_CREATIONDATATYPES,
1519 : "Byte Int8 Int16 UInt16 Int32 UInt32 "
1520 1293 : "Float32 Float64");
1521 :
1522 1293 : poDriver->SetMetadataItem(
1523 : GDAL_DMD_CREATIONOPTIONLIST,
1524 : "<CreationOptionList>"
1525 : " <Option name='PIXELTYPE' type='string' description='(deprecated, "
1526 : "use Int8 datatype) By setting this to SIGNEDBYTE, a new Byte file can "
1527 : "be forced to be written as signed byte'/>"
1528 : " <Option name='PROJ' type='string' description='ERS Projection "
1529 : "Name'/>"
1530 : " <Option name='DATUM' type='string' description='ERS Datum Name' />"
1531 : " <Option name='UNITS' type='string-select' description='ERS "
1532 : "Projection Units'>"
1533 : " <Value>METERS</Value>"
1534 : " <Value>FEET</Value>"
1535 : " </Option>"
1536 1293 : "</CreationOptionList>");
1537 :
1538 1293 : poDriver->SetMetadataItem(GDAL_DCAP_VIRTUALIO, "YES");
1539 :
1540 1293 : poDriver->pfnOpen = ERSDataset::Open;
1541 1293 : poDriver->pfnIdentify = ERSDataset::Identify;
1542 1293 : poDriver->pfnCreate = ERSDataset::Create;
1543 :
1544 1293 : GetGDALDriverManager()->RegisterDriver(poDriver);
1545 : }
|
