Line data Source code
1 : /******************************************************************************
2 : *
3 : * Project: ACE2 Driver
4 : * Purpose: Implementation of ACE2 elevation format read support.
5 : * http://tethys.eaprs.cse.dmu.ac.uk/ACE2/shared/documentation
6 : * Author: Even Rouault, <even dot rouault at spatialys.com>
7 : *
8 : ******************************************************************************
9 : * Copyright (c) 2011-2012, Even Rouault <even dot rouault at spatialys.com>
10 : *
11 : * SPDX-License-Identifier: MIT
12 : ****************************************************************************/
13 :
14 : #include "gdal_frmts.h"
15 : #include "gdal_priv.h"
16 : #include "gdal_pam.h"
17 : #include "ogr_spatialref.h"
18 : #include "rawdataset.h"
19 :
20 : static const char *const apszCategorySource[] = {
21 : "Pure SRTM (above 60deg N pure GLOBE data, below 60S pure ACE [original] "
22 : "data)",
23 : "SRTM voids filled by interpolation and/or altimeter data",
24 : "SRTM data warped using the ERS-1 Geodetic Mission",
25 : "SRTM data warped using EnviSat & ERS-2 data",
26 : "Mean lake level data derived from Altimetry",
27 : "GLOBE/ACE data warped using combined altimetry (only above 60deg N)",
28 : "Pure altimetry data (derived from ERS-1 Geodetic Mission, ERS-2 and "
29 : "EnviSat data using Delaunay Triangulation",
30 : nullptr};
31 :
32 : static const char *const apszCategoryQuality[] = {
33 : "Generic - use base datasets",
34 : "Accuracy of greater than +/- 16m",
35 : "Accuracy between +/- 16m - +/- 10m",
36 : "Accuracy between +/-10m - +/-5m",
37 : "Accuracy between +/-5m - +/-1m",
38 : "Accuracy between +/-1m",
39 : nullptr};
40 :
41 : static const char *const apszCategoryConfidence[] = {
42 : "No confidence could be derived due to lack of data",
43 : "Heights generated by interpolation",
44 : "Low confidence",
45 : "Low confidence",
46 : "Low confidence",
47 : "Medium confidence",
48 : "Medium confidence",
49 : "Medium confidence",
50 : "Medium confidence",
51 : "Medium confidence",
52 : "Medium confidence",
53 : "Medium confidence",
54 : "Medium confidence",
55 : "High confidence",
56 : "High confidence",
57 : "High confidence",
58 : "High confidence",
59 : "Inland water confidence",
60 : "Inland water confidence",
61 : "Inland water confidence",
62 : "Inland water confidence",
63 : "Inland water confidence",
64 : nullptr};
65 :
66 : /************************************************************************/
67 : /* ==================================================================== */
68 : /* ACE2Dataset */
69 : /* ==================================================================== */
70 : /************************************************************************/
71 :
72 : class ACE2Dataset final : public GDALPamDataset
73 : {
74 : friend class ACE2RasterBand;
75 :
76 : OGRSpatialReference m_oSRS{};
77 : GDALGeoTransform m_gt{};
78 :
79 : public:
80 : ACE2Dataset();
81 :
82 1 : const OGRSpatialReference *GetSpatialRef() const override
83 : {
84 1 : return &m_oSRS;
85 : }
86 :
87 : CPLErr GetGeoTransform(GDALGeoTransform >) const override;
88 :
89 : static GDALDataset *Open(GDALOpenInfo *);
90 : static int Identify(GDALOpenInfo *);
91 : };
92 :
93 : /************************************************************************/
94 : /* ==================================================================== */
95 : /* ACE2RasterBand */
96 : /* ==================================================================== */
97 : /************************************************************************/
98 :
99 : class ACE2RasterBand final : public RawRasterBand
100 : {
101 : public:
102 : ACE2RasterBand(VSILFILE *fpRaw, GDALDataType eDataType, int nXSize,
103 : int nYSize);
104 :
105 : const char *GetUnitType() override;
106 : char **GetCategoryNames() override;
107 : };
108 :
109 : /************************************************************************/
110 : /* ==================================================================== */
111 : /* ACE2Dataset */
112 : /* ==================================================================== */
113 : /************************************************************************/
114 :
115 : /************************************************************************/
116 : /* ACE2Dataset() */
117 : /************************************************************************/
118 :
119 2 : ACE2Dataset::ACE2Dataset()
120 : {
121 2 : m_oSRS.SetFromUserInput(SRS_WKT_WGS84_LAT_LONG);
122 2 : m_oSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
123 2 : }
124 :
125 : /************************************************************************/
126 : /* GetGeoTransform() */
127 : /************************************************************************/
128 :
129 1 : CPLErr ACE2Dataset::GetGeoTransform(GDALGeoTransform >) const
130 :
131 : {
132 1 : gt = m_gt;
133 1 : return CE_None;
134 : }
135 :
136 : /************************************************************************/
137 : /* ACE2RasterBand() */
138 : /************************************************************************/
139 :
140 2 : ACE2RasterBand::ACE2RasterBand(VSILFILE *fpRawIn, GDALDataType eDataTypeIn,
141 2 : int nXSize, int nYSize)
142 : : RawRasterBand(fpRawIn, 0, GDALGetDataTypeSizeBytes(eDataTypeIn),
143 4 : nXSize * GDALGetDataTypeSizeBytes(eDataTypeIn), eDataTypeIn,
144 2 : CPL_IS_LSB, nXSize, nYSize, RawRasterBand::OwnFP::YES)
145 : {
146 2 : }
147 :
148 : /************************************************************************/
149 : /* GetUnitType() */
150 : /************************************************************************/
151 :
152 0 : const char *ACE2RasterBand::GetUnitType()
153 : {
154 0 : if (eDataType == GDT_Float32)
155 0 : return "m";
156 :
157 0 : return "";
158 : }
159 :
160 : /************************************************************************/
161 : /* GetCategoryNames() */
162 : /************************************************************************/
163 :
164 0 : char **ACE2RasterBand::GetCategoryNames()
165 : {
166 0 : if (eDataType != GDT_Int16)
167 0 : return nullptr;
168 :
169 0 : const char *pszName = poDS->GetDescription();
170 :
171 0 : if (strstr(pszName, "_SOURCE_"))
172 0 : return const_cast<char **>(apszCategorySource);
173 0 : if (strstr(pszName, "_QUALITY_"))
174 0 : return const_cast<char **>(apszCategoryQuality);
175 0 : if (strstr(pszName, "_CONF_"))
176 0 : return const_cast<char **>(apszCategoryConfidence);
177 :
178 0 : return nullptr;
179 : }
180 :
181 : /************************************************************************/
182 : /* Identify() */
183 : /************************************************************************/
184 :
185 58968 : int ACE2Dataset::Identify(GDALOpenInfo *poOpenInfo)
186 :
187 : {
188 117932 : if (!(poOpenInfo->IsExtensionEqualToCI("ACE2") ||
189 58964 : strstr(poOpenInfo->pszFilename, ".ACE2.gz") ||
190 58964 : strstr(poOpenInfo->pszFilename, ".ace2.gz")))
191 58964 : return FALSE;
192 :
193 4 : return TRUE;
194 : }
195 :
196 : /************************************************************************/
197 : /* Open() */
198 : /************************************************************************/
199 :
200 2 : GDALDataset *ACE2Dataset::Open(GDALOpenInfo *poOpenInfo)
201 :
202 : {
203 2 : if (!Identify(poOpenInfo))
204 0 : return nullptr;
205 :
206 4 : const std::string osBasename = CPLGetBasenameSafe(poOpenInfo->pszFilename);
207 2 : const char *pszBasename = osBasename.c_str();
208 :
209 2 : if (strlen(pszBasename) < 7)
210 0 : return nullptr;
211 :
212 : /* Determine southwest coordinates from filename */
213 :
214 : /* e.g. 30S120W_5M.ACE2 */
215 2 : char pszLatLonValueString[4] = {'\0'};
216 2 : memset(pszLatLonValueString, 0, 4);
217 : // cppcheck-suppress redundantCopy
218 2 : strncpy(pszLatLonValueString, &pszBasename[0], 2);
219 2 : int southWestLat = atoi(pszLatLonValueString);
220 2 : memset(pszLatLonValueString, 0, 4);
221 : // cppcheck-suppress redundantCopy
222 2 : strncpy(pszLatLonValueString, &pszBasename[3], 3);
223 2 : int southWestLon = atoi(pszLatLonValueString);
224 :
225 2 : if (pszBasename[2] == 'N' || pszBasename[2] == 'n')
226 : /*southWestLat = southWestLat*/;
227 0 : else if (pszBasename[2] == 'S' || pszBasename[2] == 's')
228 0 : southWestLat = southWestLat * -1;
229 : else
230 0 : return nullptr;
231 :
232 2 : if (pszBasename[6] == 'E' || pszBasename[6] == 'e')
233 : /*southWestLon = southWestLon*/;
234 0 : else if (pszBasename[6] == 'W' || pszBasename[6] == 'w')
235 0 : southWestLon = southWestLon * -1;
236 : else
237 0 : return nullptr;
238 :
239 2 : GDALDataType eDT = GDT_Unknown;
240 2 : if (strstr(pszBasename, "_CONF_") || strstr(pszBasename, "_QUALITY_") ||
241 2 : strstr(pszBasename, "_SOURCE_"))
242 0 : eDT = GDT_Int16;
243 : else
244 2 : eDT = GDT_Float32;
245 2 : const int nWordSize = GDALGetDataTypeSizeBytes(eDT);
246 :
247 : VSIStatBufL sStat;
248 2 : if (strstr(pszBasename, "_5M"))
249 2 : sStat.st_size = 180 * 180 * nWordSize;
250 0 : else if (strstr(pszBasename, "_30S"))
251 0 : sStat.st_size = 1800 * 1800 * nWordSize;
252 0 : else if (strstr(pszBasename, "_9S"))
253 0 : sStat.st_size = 6000 * 6000 * nWordSize;
254 0 : else if (strstr(pszBasename, "_3S"))
255 0 : sStat.st_size = 18000 * 18000 * nWordSize;
256 : /* Check file size otherwise */
257 0 : else if (VSIStatL(poOpenInfo->pszFilename, &sStat) != 0)
258 : {
259 0 : return nullptr;
260 : }
261 :
262 2 : int nXSize = 0;
263 2 : int nYSize = 0;
264 :
265 2 : double dfPixelSize = 0;
266 2 : if (sStat.st_size == 180 * 180 * nWordSize)
267 : {
268 : /* 5 minute */
269 2 : nXSize = 180;
270 2 : nYSize = 180;
271 2 : dfPixelSize = 5.0 / 60;
272 : }
273 0 : else if (sStat.st_size == 1800 * 1800 * nWordSize)
274 : {
275 : /* 30 s */
276 0 : nXSize = 1800;
277 0 : nYSize = 1800;
278 0 : dfPixelSize = 30.0 / 3600;
279 : }
280 0 : else if (sStat.st_size == 6000 * 6000 * nWordSize)
281 : {
282 : /* 9 s */
283 0 : nXSize = 6000;
284 0 : nYSize = 6000;
285 0 : dfPixelSize = 9.0 / 3600;
286 : }
287 0 : else if (sStat.st_size == 18000 * 18000 * nWordSize)
288 : {
289 : /* 3 s */
290 0 : nXSize = 18000;
291 0 : nYSize = 18000;
292 0 : dfPixelSize = 3.0 / 3600;
293 : }
294 : else
295 0 : return nullptr;
296 :
297 : /* -------------------------------------------------------------------- */
298 : /* Open file. */
299 : /* -------------------------------------------------------------------- */
300 :
301 4 : CPLString osFilename = poOpenInfo->pszFilename;
302 2 : if ((strstr(poOpenInfo->pszFilename, ".ACE2.gz") ||
303 2 : strstr(poOpenInfo->pszFilename, ".ace2.gz")) &&
304 0 : !STARTS_WITH(poOpenInfo->pszFilename, "/vsigzip/"))
305 0 : osFilename = "/vsigzip/" + osFilename;
306 :
307 2 : VSILFILE *fpImage = VSIFOpenL(osFilename, "rb+");
308 2 : if (fpImage == nullptr)
309 0 : return nullptr;
310 :
311 : /* -------------------------------------------------------------------- */
312 : /* Create the dataset. */
313 : /* -------------------------------------------------------------------- */
314 4 : auto poDS = std::make_unique<ACE2Dataset>();
315 :
316 2 : poDS->nRasterXSize = nXSize;
317 2 : poDS->nRasterYSize = nYSize;
318 :
319 2 : poDS->m_gt[0] = southWestLon;
320 2 : poDS->m_gt[1] = dfPixelSize;
321 2 : poDS->m_gt[2] = 0.0;
322 2 : poDS->m_gt[3] = southWestLat + nYSize * dfPixelSize;
323 2 : poDS->m_gt[4] = 0.0;
324 2 : poDS->m_gt[5] = -dfPixelSize;
325 :
326 : /* -------------------------------------------------------------------- */
327 : /* Create band information objects */
328 : /* -------------------------------------------------------------------- */
329 : auto poBand =
330 4 : std::make_unique<ACE2RasterBand>(fpImage, eDT, nXSize, nYSize);
331 2 : if (!poBand->IsValid())
332 0 : return nullptr;
333 2 : poDS->SetBand(1, std::move(poBand));
334 :
335 : /* -------------------------------------------------------------------- */
336 : /* Initialize any PAM information. */
337 : /* -------------------------------------------------------------------- */
338 2 : poDS->SetDescription(poOpenInfo->pszFilename);
339 2 : poDS->TryLoadXML();
340 :
341 : /* -------------------------------------------------------------------- */
342 : /* Check for overviews. */
343 : /* -------------------------------------------------------------------- */
344 2 : poDS->oOvManager.Initialize(poDS.get(), poOpenInfo->pszFilename);
345 :
346 2 : return poDS.release();
347 : }
348 :
349 : /************************************************************************/
350 : /* GDALRegister_ACE2() */
351 : /************************************************************************/
352 :
353 2033 : void GDALRegister_ACE2()
354 :
355 : {
356 2033 : if (GDALGetDriverByName("ACE2") != nullptr)
357 283 : return;
358 :
359 1750 : GDALDriver *poDriver = new GDALDriver();
360 :
361 1750 : poDriver->SetDescription("ACE2");
362 1750 : poDriver->SetMetadataItem(GDAL_DCAP_RASTER, "YES");
363 1750 : poDriver->SetMetadataItem(GDAL_DMD_LONGNAME, "ACE2");
364 1750 : poDriver->SetMetadataItem(GDAL_DMD_HELPTOPIC, "drivers/raster/ace2.html");
365 1750 : poDriver->SetMetadataItem(GDAL_DMD_EXTENSION, "ACE2");
366 1750 : poDriver->SetMetadataItem(GDAL_DCAP_VIRTUALIO, "YES");
367 :
368 1750 : poDriver->pfnOpen = ACE2Dataset::Open;
369 1750 : poDriver->pfnIdentify = ACE2Dataset::Identify;
370 :
371 1750 : GetGDALDriverManager()->RegisterDriver(poDriver);
372 : }
|
