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