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