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 59489 : int ERSDataset::Identify(GDALOpenInfo *poOpenInfo)
786 :
787 : {
788 : /* -------------------------------------------------------------------- */
789 : /* We assume the user selects the .ers file. */
790 : /* -------------------------------------------------------------------- */
791 59489 : CPLString osHeader(reinterpret_cast<const char *>(poOpenInfo->pabyHeader),
792 118978 : poOpenInfo->nHeaderBytes);
793 :
794 59488 : 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 59488 : if (osHeader.ifind("DatasetHeader ") != std::string::npos)
804 89 : return TRUE;
805 :
806 59400 : 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.length() == 0) // just strip off extension.
991 : {
992 61 : osDataFile = CPLGetFilename(poOpenInfo->pszFilename);
993 61 : osDataFile = osDataFile.substr(0, osDataFile.find_last_of('.'));
994 : }
995 :
996 124 : CPLString osDataFilePath = CPLFormFilenameSafe(osPath, osDataFile, nullptr);
997 :
998 : /* -------------------------------------------------------------------- */
999 : /* DataSetType = Translated files are links to things like ecw */
1000 : /* files. */
1001 : /* -------------------------------------------------------------------- */
1002 62 : if (EQUAL(poHeader->Find("DataSetType", ""), "Translated"))
1003 : {
1004 2 : nRecLevel++;
1005 2 : poDS->poDepFile = GDALDataset::FromHandle(
1006 : GDALOpen(osDataFilePath, poOpenInfo->eAccess));
1007 2 : nRecLevel--;
1008 :
1009 2 : if (poDS->poDepFile != nullptr &&
1010 1 : poDS->poDepFile->GetRasterXSize() == poDS->GetRasterXSize() &&
1011 4 : poDS->poDepFile->GetRasterYSize() == poDS->GetRasterYSize() &&
1012 1 : poDS->poDepFile->GetRasterCount() >= nBands)
1013 : {
1014 2 : for (int iBand = 0; iBand < nBands; iBand++)
1015 : {
1016 : // Assume pixel interleaved.
1017 2 : poDS->SetBand(iBand + 1,
1018 : new ERSProxyRasterBand(
1019 1 : poDS->poDepFile->GetRasterBand(iBand + 1)));
1020 : }
1021 : }
1022 : else
1023 : {
1024 1 : delete poDS->poDepFile;
1025 1 : poDS->poDepFile = nullptr;
1026 : }
1027 : }
1028 :
1029 : /* ==================================================================== */
1030 : /* While ERStorage indicates a raw file. */
1031 : /* ==================================================================== */
1032 60 : else if (EQUAL(poHeader->Find("DataSetType", ""), "ERStorage"))
1033 : {
1034 : // Open data file.
1035 60 : if (poOpenInfo->eAccess == GA_Update)
1036 38 : poDS->fpImage = VSIFOpenL(osDataFilePath, "r+");
1037 : else
1038 22 : poDS->fpImage = VSIFOpenL(osDataFilePath, "r");
1039 :
1040 60 : poDS->osRawFilename = std::move(osDataFilePath);
1041 :
1042 60 : if (poDS->fpImage != nullptr && nBands > 0)
1043 : {
1044 60 : int iWordSize = GDALGetDataTypeSizeBytes(eType);
1045 :
1046 60 : const auto knIntMax = std::numeric_limits<int>::max();
1047 120 : if (nBands > knIntMax / iWordSize ||
1048 60 : poDS->nRasterXSize > knIntMax / (nBands * iWordSize))
1049 : {
1050 0 : CPLError(CE_Failure, CPLE_AppDefined,
1051 : "int overflow: too large nBands and/or nRasterXSize");
1052 0 : return nullptr;
1053 : }
1054 :
1055 60 : if (!RAWDatasetCheckMemoryUsage(
1056 60 : poDS->nRasterXSize, poDS->nRasterYSize, nBands, iWordSize,
1057 60 : iWordSize, iWordSize * nBands * poDS->nRasterXSize,
1058 : nHeaderOffset,
1059 60 : static_cast<vsi_l_offset>(iWordSize) * poDS->nRasterXSize,
1060 60 : poDS->fpImage))
1061 : {
1062 0 : return nullptr;
1063 : }
1064 60 : if (nHeaderOffset > std::numeric_limits<GIntBig>::max() -
1065 60 : static_cast<GIntBig>(nBands - 1) *
1066 60 : iWordSize * poDS->nRasterXSize)
1067 : {
1068 0 : CPLError(CE_Failure, CPLE_AppDefined,
1069 : "int overflow: too large nHeaderOffset");
1070 0 : return nullptr;
1071 : }
1072 :
1073 162 : for (int iBand = 0; iBand < nBands; iBand++)
1074 : {
1075 : // Assume pixel interleaved.
1076 : auto poBand = std::make_unique<ERSRasterBand>(
1077 102 : poDS.get(), iBand + 1, poDS->fpImage,
1078 0 : nHeaderOffset + static_cast<vsi_l_offset>(iWordSize) *
1079 102 : iBand * poDS->nRasterXSize,
1080 102 : iWordSize, iWordSize * nBands * poDS->nRasterXSize, eType,
1081 102 : bNative);
1082 102 : if (!poBand->IsValid())
1083 0 : return nullptr;
1084 102 : poDS->SetBand(iBand + 1, std::move(poBand));
1085 : }
1086 : }
1087 : }
1088 :
1089 : /* -------------------------------------------------------------------- */
1090 : /* Otherwise we have an error! */
1091 : /* -------------------------------------------------------------------- */
1092 62 : if (poDS->nBands == 0)
1093 : {
1094 1 : return nullptr;
1095 : }
1096 :
1097 : /* -------------------------------------------------------------------- */
1098 : /* Look for band descriptions. */
1099 : /* -------------------------------------------------------------------- */
1100 61 : ERSHdrNode *poRI = poHeader->FindNode("RasterInfo");
1101 :
1102 357 : for (int iChild = 0, iBand = 0;
1103 357 : poRI != nullptr && iChild < poRI->nItemCount && iBand < poDS->nBands;
1104 : iChild++)
1105 : {
1106 296 : if (poRI->papoItemChild[iChild] != nullptr &&
1107 33 : EQUAL(poRI->papszItemName[iChild], "BandId"))
1108 : {
1109 : const char *pszValue =
1110 12 : poRI->papoItemChild[iChild]->Find("Value", nullptr);
1111 :
1112 12 : iBand++;
1113 12 : if (pszValue)
1114 : {
1115 12 : CPLPushErrorHandler(CPLQuietErrorHandler);
1116 12 : poDS->GetRasterBand(iBand)->SetDescription(pszValue);
1117 12 : CPLPopErrorHandler();
1118 : }
1119 :
1120 12 : pszValue = poRI->papoItemChild[iChild]->Find("Units", nullptr);
1121 12 : if (pszValue)
1122 : {
1123 4 : CPLPushErrorHandler(CPLQuietErrorHandler);
1124 4 : poDS->GetRasterBand(iBand)->SetUnitType(pszValue);
1125 4 : CPLPopErrorHandler();
1126 : }
1127 : }
1128 : }
1129 :
1130 : /* -------------------------------------------------------------------- */
1131 : /* Look for projection. */
1132 : /* -------------------------------------------------------------------- */
1133 61 : poDS->osProj = poHeader->Find("CoordinateSpace.Projection", "");
1134 61 : poDS->osDatum = poHeader->Find("CoordinateSpace.Datum", "");
1135 61 : poDS->osUnits = poHeader->Find("CoordinateSpace.Units", "");
1136 :
1137 219 : poDS->m_oSRS.importFromERM(
1138 61 : !poDS->osProj.empty() ? poDS->osProj.c_str() : "RAW",
1139 61 : !poDS->osDatum.empty() ? poDS->osDatum.c_str() : "WGS84",
1140 61 : !poDS->osUnits.empty() ? poDS->osUnits.c_str() : "METERS");
1141 :
1142 : /* -------------------------------------------------------------------- */
1143 : /* Look for the geotransform. */
1144 : /* -------------------------------------------------------------------- */
1145 61 : if (poHeader->Find("RasterInfo.RegistrationCoord.Eastings", nullptr))
1146 : {
1147 6 : poDS->bGotTransform = TRUE;
1148 6 : poDS->m_gt[0] = CPLAtof(
1149 : poHeader->Find("RasterInfo.RegistrationCoord.Eastings", ""));
1150 6 : poDS->m_gt[1] =
1151 6 : CPLAtof(poHeader->Find("RasterInfo.CellInfo.Xdimension", "1.0"));
1152 6 : poDS->m_gt[2] = 0.0;
1153 6 : poDS->m_gt[3] = CPLAtof(
1154 : poHeader->Find("RasterInfo.RegistrationCoord.Northings", ""));
1155 6 : poDS->m_gt[4] = 0.0;
1156 6 : poDS->m_gt[5] =
1157 6 : -CPLAtof(poHeader->Find("RasterInfo.CellInfo.Ydimension", "1.0"));
1158 : }
1159 60 : else if (poHeader->Find("RasterInfo.RegistrationCoord.Latitude", nullptr) &&
1160 5 : poHeader->Find("RasterInfo.CellInfo.Xdimension", nullptr))
1161 : {
1162 5 : poDS->bGotTransform = TRUE;
1163 5 : poDS->m_gt[0] = ERSDMS2Dec(
1164 : poHeader->Find("RasterInfo.RegistrationCoord.Longitude", ""));
1165 5 : poDS->m_gt[1] =
1166 5 : CPLAtof(poHeader->Find("RasterInfo.CellInfo.Xdimension", ""));
1167 5 : poDS->m_gt[2] = 0.0;
1168 5 : poDS->m_gt[3] = ERSDMS2Dec(
1169 : poHeader->Find("RasterInfo.RegistrationCoord.Latitude", ""));
1170 5 : poDS->m_gt[4] = 0.0;
1171 5 : poDS->m_gt[5] =
1172 5 : -CPLAtof(poHeader->Find("RasterInfo.CellInfo.Ydimension", ""));
1173 : }
1174 :
1175 : /* -------------------------------------------------------------------- */
1176 : /* Adjust if we have a registration cell. */
1177 : /* -------------------------------------------------------------------- */
1178 :
1179 : /* http://geospatial.intergraph.com/Libraries/Tech_Docs/ERDAS_ER_Mapper_Customization_Guide.sflb.ashx
1180 : */
1181 : /* Page 27 : */
1182 : /* RegistrationCellX and RegistrationCellY : The image X and Y
1183 : coordinates of the cell which corresponds to the Registration
1184 : Coordinate. Note that the RegistrationCellX and
1185 : RegistrationCellY can be fractional values. If
1186 : RegistrationCellX and RegistrationCellY are not specified,
1187 : they are assumed to be (0,0), which is the top left corner of the
1188 : image.
1189 : */
1190 : double dfCellX =
1191 61 : CPLAtof(poHeader->Find("RasterInfo.RegistrationCellX", "0"));
1192 : double dfCellY =
1193 61 : CPLAtof(poHeader->Find("RasterInfo.RegistrationCellY", "0"));
1194 :
1195 61 : if (poDS->bGotTransform)
1196 : {
1197 11 : poDS->m_gt[0] -= dfCellX * poDS->m_gt[1] + dfCellY * poDS->m_gt[2];
1198 11 : poDS->m_gt[3] -= dfCellX * poDS->m_gt[4] + dfCellY * poDS->m_gt[5];
1199 : }
1200 :
1201 : /* -------------------------------------------------------------------- */
1202 : /* Check for null values. */
1203 : /* -------------------------------------------------------------------- */
1204 61 : if (poHeader->Find("RasterInfo.NullCellValue", nullptr))
1205 : {
1206 8 : poDS->bHasNoDataValue = TRUE;
1207 16 : poDS->dfNoDataValue =
1208 8 : CPLAtofM(poHeader->Find("RasterInfo.NullCellValue"));
1209 :
1210 8 : if (poDS->poDepFile != nullptr)
1211 : {
1212 0 : CPLPushErrorHandler(CPLQuietErrorHandler);
1213 :
1214 0 : for (int iBand = 1; iBand <= poDS->nBands; iBand++)
1215 0 : poDS->GetRasterBand(iBand)->SetNoDataValue(poDS->dfNoDataValue);
1216 :
1217 0 : CPLPopErrorHandler();
1218 : }
1219 : }
1220 :
1221 : /* -------------------------------------------------------------------- */
1222 : /* Do we have an "All" region? */
1223 : /* -------------------------------------------------------------------- */
1224 61 : ERSHdrNode *poAll = nullptr;
1225 :
1226 366 : for (int iChild = 0; poRI != nullptr && iChild < poRI->nItemCount; iChild++)
1227 : {
1228 305 : if (poRI->papoItemChild[iChild] != nullptr &&
1229 39 : EQUAL(poRI->papszItemName[iChild], "RegionInfo"))
1230 : {
1231 3 : if (EQUAL(poRI->papoItemChild[iChild]->Find("RegionName", ""),
1232 : "All"))
1233 3 : poAll = poRI->papoItemChild[iChild];
1234 : }
1235 : }
1236 :
1237 : /* -------------------------------------------------------------------- */
1238 : /* Do we have statistics? */
1239 : /* -------------------------------------------------------------------- */
1240 61 : if (poAll && poAll->FindNode("Stats"))
1241 : {
1242 3 : CPLPushErrorHandler(CPLQuietErrorHandler);
1243 :
1244 6 : for (int iBand = 1; iBand <= poDS->nBands; iBand++)
1245 : {
1246 : const char *pszValue =
1247 3 : poAll->FindElem("Stats.MinimumValue", iBand - 1);
1248 :
1249 3 : if (pszValue)
1250 3 : poDS->GetRasterBand(iBand)->SetMetadataItem(
1251 3 : "STATISTICS_MINIMUM", pszValue);
1252 :
1253 3 : pszValue = poAll->FindElem("Stats.MaximumValue", iBand - 1);
1254 :
1255 3 : if (pszValue)
1256 3 : poDS->GetRasterBand(iBand)->SetMetadataItem(
1257 3 : "STATISTICS_MAXIMUM", pszValue);
1258 :
1259 3 : pszValue = poAll->FindElem("Stats.MeanValue", iBand - 1);
1260 :
1261 3 : if (pszValue)
1262 3 : poDS->GetRasterBand(iBand)->SetMetadataItem("STATISTICS_MEAN",
1263 3 : pszValue);
1264 :
1265 3 : pszValue = poAll->FindElem("Stats.MedianValue", iBand - 1);
1266 :
1267 3 : if (pszValue)
1268 3 : poDS->GetRasterBand(iBand)->SetMetadataItem("STATISTICS_MEDIAN",
1269 3 : pszValue);
1270 : }
1271 :
1272 3 : CPLPopErrorHandler();
1273 : }
1274 :
1275 : /* -------------------------------------------------------------------- */
1276 : /* Do we have GCPs. */
1277 : /* -------------------------------------------------------------------- */
1278 61 : if (poHeader->FindNode("RasterInfo.WarpControl"))
1279 2 : poDS->ReadGCPs();
1280 :
1281 : /* -------------------------------------------------------------------- */
1282 : /* Initialize any PAM information. */
1283 : /* -------------------------------------------------------------------- */
1284 61 : poDS->SetDescription(poOpenInfo->pszFilename);
1285 61 : poDS->TryLoadXML();
1286 :
1287 : // if no SR in xml, try aux
1288 61 : const OGRSpatialReference *poSRS = poDS->GDALPamDataset::GetSpatialRef();
1289 61 : if (poSRS == nullptr)
1290 : {
1291 : // try aux
1292 : auto poAuxDS = std::unique_ptr<GDALDataset>(GDALFindAssociatedAuxFile(
1293 122 : poOpenInfo->pszFilename, GA_ReadOnly, poDS.get()));
1294 61 : if (poAuxDS)
1295 : {
1296 0 : poSRS = poAuxDS->GetSpatialRef();
1297 0 : if (poSRS)
1298 : {
1299 0 : poDS->m_oSRS = *poSRS;
1300 : }
1301 : }
1302 : }
1303 : /* -------------------------------------------------------------------- */
1304 : /* Check for overviews. */
1305 : /* -------------------------------------------------------------------- */
1306 61 : poDS->oOvManager.Initialize(poDS.get(), poOpenInfo->pszFilename);
1307 :
1308 61 : return poDS.release();
1309 : }
1310 :
1311 : /************************************************************************/
1312 : /* Create() */
1313 : /************************************************************************/
1314 :
1315 67 : GDALDataset *ERSDataset::Create(const char *pszFilename, int nXSize, int nYSize,
1316 : int nBandsIn, GDALDataType eType,
1317 : char **papszOptions)
1318 :
1319 : {
1320 : /* -------------------------------------------------------------------- */
1321 : /* Verify settings. */
1322 : /* -------------------------------------------------------------------- */
1323 67 : if (nBandsIn <= 0)
1324 : {
1325 1 : CPLError(CE_Failure, CPLE_NotSupported,
1326 : "ERS driver does not support %d bands.\n", nBandsIn);
1327 1 : return nullptr;
1328 : }
1329 :
1330 66 : if (eType != GDT_Byte && eType != GDT_Int8 && eType != GDT_Int16 &&
1331 29 : eType != GDT_UInt16 && eType != GDT_Int32 && eType != GDT_UInt32 &&
1332 19 : eType != GDT_Float32 && eType != GDT_Float64)
1333 : {
1334 16 : CPLError(
1335 : CE_Failure, CPLE_AppDefined,
1336 : "The ERS driver does not supporting creating files of types %s.",
1337 : GDALGetDataTypeName(eType));
1338 16 : return nullptr;
1339 : }
1340 :
1341 : /* -------------------------------------------------------------------- */
1342 : /* Work out the name we want to use for the .ers and binary */
1343 : /* data files. */
1344 : /* -------------------------------------------------------------------- */
1345 100 : CPLString osBinFile, osErsFile;
1346 :
1347 50 : if (EQUAL(CPLGetExtensionSafe(pszFilename).c_str(), "ers"))
1348 : {
1349 7 : osErsFile = pszFilename;
1350 7 : osBinFile = osErsFile.substr(0, osErsFile.length() - 4);
1351 : }
1352 : else
1353 : {
1354 43 : osBinFile = pszFilename;
1355 43 : osErsFile = osBinFile + ".ers";
1356 : }
1357 :
1358 : /* -------------------------------------------------------------------- */
1359 : /* Work out some values we will write. */
1360 : /* -------------------------------------------------------------------- */
1361 50 : const char *pszCellType = "Unsigned8BitInteger";
1362 50 : CPL_IGNORE_RET_VAL(pszCellType); // Make CSA happy
1363 :
1364 50 : if (eType == GDT_Byte)
1365 28 : pszCellType = "Unsigned8BitInteger";
1366 22 : else if (eType == GDT_Int8)
1367 3 : pszCellType = "Signed8BitInteger";
1368 19 : else if (eType == GDT_Int16)
1369 3 : pszCellType = "Signed16BitInteger";
1370 16 : else if (eType == GDT_UInt16)
1371 3 : pszCellType = "Unsigned16BitInteger";
1372 13 : else if (eType == GDT_Int32)
1373 3 : pszCellType = "Signed32BitInteger";
1374 10 : else if (eType == GDT_UInt32)
1375 3 : pszCellType = "Unsigned32BitInteger";
1376 7 : else if (eType == GDT_Float32)
1377 4 : pszCellType = "IEEE4ByteReal";
1378 3 : else if (eType == GDT_Float64)
1379 3 : pszCellType = "IEEE8ByteReal";
1380 : else
1381 : {
1382 0 : CPLAssert(false);
1383 : }
1384 :
1385 : /* -------------------------------------------------------------------- */
1386 : /* Handling for signed eight bit data. */
1387 : /* -------------------------------------------------------------------- */
1388 50 : const char *pszPixelType = CSLFetchNameValue(papszOptions, "PIXELTYPE");
1389 50 : if (pszPixelType && EQUAL(pszPixelType, "SIGNEDBYTE") && eType == GDT_Byte)
1390 0 : pszCellType = "Signed8BitInteger";
1391 :
1392 : /* -------------------------------------------------------------------- */
1393 : /* Write binary file. */
1394 : /* -------------------------------------------------------------------- */
1395 50 : VSILFILE *fpBin = VSIFOpenL(osBinFile, "w");
1396 :
1397 50 : if (fpBin == nullptr)
1398 : {
1399 3 : CPLError(CE_Failure, CPLE_FileIO, "Failed to create %s:\n%s",
1400 3 : osBinFile.c_str(), VSIStrerror(errno));
1401 3 : return nullptr;
1402 : }
1403 :
1404 47 : const GUIntBig nSize = static_cast<GUIntBig>(nXSize) * nYSize * nBandsIn *
1405 47 : GDALGetDataTypeSizeBytes(eType);
1406 47 : GByte byZero = 0;
1407 94 : if (VSIFSeekL(fpBin, nSize - 1, SEEK_SET) != 0 ||
1408 47 : VSIFWriteL(&byZero, 1, 1, fpBin) != 1)
1409 : {
1410 10 : CPLError(CE_Failure, CPLE_FileIO, "Failed to write %s:\n%s",
1411 10 : osBinFile.c_str(), VSIStrerror(errno));
1412 10 : VSIFCloseL(fpBin);
1413 10 : return nullptr;
1414 : }
1415 37 : VSIFCloseL(fpBin);
1416 :
1417 : /* -------------------------------------------------------------------- */
1418 : /* Try writing header file. */
1419 : /* -------------------------------------------------------------------- */
1420 37 : VSILFILE *fpERS = VSIFOpenL(osErsFile, "w");
1421 :
1422 37 : if (fpERS == nullptr)
1423 : {
1424 0 : CPLError(CE_Failure, CPLE_FileIO, "Failed to create %s:\n%s",
1425 0 : osErsFile.c_str(), VSIStrerror(errno));
1426 0 : return nullptr;
1427 : }
1428 :
1429 37 : VSIFPrintfL(fpERS, "DatasetHeader Begin\n");
1430 37 : VSIFPrintfL(fpERS, "\tVersion\t\t = \"6.0\"\n");
1431 37 : VSIFPrintfL(fpERS, "\tName\t\t= \"%s\"\n", CPLGetFilename(osErsFile));
1432 :
1433 : // Last updated requires timezone info which we don't necessarily get
1434 : // get from VSICTime() so perhaps it is better to omit this.
1435 : // VSIFPrintfL( fpERS, "\tLastUpdated\t= %s",
1436 : // VSICTime( VSITime( NULL ) ) );
1437 :
1438 37 : VSIFPrintfL(fpERS, "\tDataSetType\t= ERStorage\n");
1439 37 : VSIFPrintfL(fpERS, "\tDataType\t= Raster\n");
1440 37 : VSIFPrintfL(fpERS, "\tByteOrder\t= LSBFirst\n");
1441 37 : VSIFPrintfL(fpERS, "\tRasterInfo Begin\n");
1442 37 : VSIFPrintfL(fpERS, "\t\tCellType\t= %s\n", pszCellType);
1443 37 : VSIFPrintfL(fpERS, "\t\tNrOfLines\t= %d\n", nYSize);
1444 37 : VSIFPrintfL(fpERS, "\t\tNrOfCellsPerLine\t= %d\n", nXSize);
1445 37 : VSIFPrintfL(fpERS, "\t\tNrOfBands\t= %d\n", nBandsIn);
1446 37 : VSIFPrintfL(fpERS, "\tRasterInfo End\n");
1447 37 : if (VSIFPrintfL(fpERS, "DatasetHeader End\n") < 17)
1448 : {
1449 0 : CPLError(CE_Failure, CPLE_FileIO, "Failed to write %s:\n%s",
1450 0 : osErsFile.c_str(), VSIStrerror(errno));
1451 0 : return nullptr;
1452 : }
1453 :
1454 37 : VSIFCloseL(fpERS);
1455 :
1456 : /* -------------------------------------------------------------------- */
1457 : /* Reopen. */
1458 : /* -------------------------------------------------------------------- */
1459 74 : GDALOpenInfo oOpenInfo(osErsFile, GA_Update);
1460 37 : ERSDataset *poDS = cpl::down_cast<ERSDataset *>(Open(&oOpenInfo));
1461 37 : if (poDS == nullptr)
1462 0 : return nullptr;
1463 :
1464 : /* -------------------------------------------------------------------- */
1465 : /* Fetch DATUM, PROJ and UNITS creation option */
1466 : /* -------------------------------------------------------------------- */
1467 37 : const char *pszDatum = CSLFetchNameValue(papszOptions, "DATUM");
1468 37 : if (pszDatum)
1469 : {
1470 2 : poDS->osDatumForced = pszDatum;
1471 2 : poDS->osDatum = pszDatum;
1472 : }
1473 37 : const char *pszProj = CSLFetchNameValue(papszOptions, "PROJ");
1474 37 : if (pszProj)
1475 : {
1476 2 : poDS->osProjForced = pszProj;
1477 2 : poDS->osProj = pszProj;
1478 : }
1479 37 : const char *pszUnits = CSLFetchNameValue(papszOptions, "UNITS");
1480 37 : if (pszUnits)
1481 : {
1482 2 : poDS->osUnitsForced = pszUnits;
1483 2 : poDS->osUnits = pszUnits;
1484 : }
1485 :
1486 37 : if (pszDatum || pszProj || pszUnits)
1487 : {
1488 2 : poDS->WriteProjectionInfo(pszProj ? pszProj : "RAW",
1489 : pszDatum ? pszDatum : "RAW",
1490 : pszUnits ? pszUnits : "METERS");
1491 : }
1492 :
1493 37 : return poDS;
1494 : }
1495 :
1496 : /************************************************************************/
1497 : /* GDALRegister_ERS() */
1498 : /************************************************************************/
1499 :
1500 1964 : void GDALRegister_ERS()
1501 :
1502 : {
1503 1964 : if (GDALGetDriverByName("ERS") != nullptr)
1504 283 : return;
1505 :
1506 1681 : GDALDriver *poDriver = new GDALDriver();
1507 :
1508 1681 : poDriver->SetDescription("ERS");
1509 1681 : poDriver->SetMetadataItem(GDAL_DCAP_RASTER, "YES");
1510 1681 : poDriver->SetMetadataItem(GDAL_DMD_LONGNAME, "ERMapper .ers Labelled");
1511 1681 : poDriver->SetMetadataItem(GDAL_DMD_HELPTOPIC, "drivers/raster/ers.html");
1512 1681 : poDriver->SetMetadataItem(GDAL_DMD_EXTENSION, "ers");
1513 1681 : poDriver->SetMetadataItem(GDAL_DMD_CREATIONDATATYPES,
1514 : "Byte Int8 Int16 UInt16 Int32 UInt32 "
1515 1681 : "Float32 Float64");
1516 :
1517 1681 : poDriver->SetMetadataItem(
1518 : GDAL_DMD_CREATIONOPTIONLIST,
1519 : "<CreationOptionList>"
1520 : " <Option name='PIXELTYPE' type='string' description='(deprecated, "
1521 : "use Int8 datatype) By setting this to SIGNEDBYTE, a new Byte file can "
1522 : "be forced to be written as signed byte'/>"
1523 : " <Option name='PROJ' type='string' description='ERS Projection "
1524 : "Name'/>"
1525 : " <Option name='DATUM' type='string' description='ERS Datum Name' />"
1526 : " <Option name='UNITS' type='string-select' description='ERS "
1527 : "Projection Units'>"
1528 : " <Value>METERS</Value>"
1529 : " <Value>FEET</Value>"
1530 : " </Option>"
1531 1681 : "</CreationOptionList>");
1532 :
1533 1681 : poDriver->SetMetadataItem(GDAL_DCAP_VIRTUALIO, "YES");
1534 :
1535 1681 : poDriver->pfnOpen = ERSDataset::Open;
1536 1681 : poDriver->pfnIdentify = ERSDataset::Identify;
1537 1681 : poDriver->pfnCreate = ERSDataset::Create;
1538 :
1539 1681 : GetGDALDriverManager()->RegisterDriver(poDriver);
1540 : }
|