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