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