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