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