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