1 : /******************************************************************************
2 : * $Id: ace2dataset.cpp 23772 2012-01-21 08:31:07Z rouault $
3 : *
4 : * Project: ACE2 Driver
5 : * Purpose: Implementation of ACE2 elevation format read support.
6 : * http://tethys.eaprs.cse.dmu.ac.uk/ACE2/shared/documentation
7 : * Author: Even Rouault, <even dot rouault at mines dash paris dot org>
8 : *
9 : ******************************************************************************
10 : * Copyright (c) 2011, Even Rouault
11 : *
12 : * Permission is hereby granted, free of charge, to any person obtaining a
13 : * copy of this software and associated documentation files (the "Software"),
14 : * to deal in the Software without restriction, including without limitation
15 : * the rights to use, copy, modify, merge, publish, distribute, sublicense,
16 : * and/or sell copies of the Software, and to permit persons to whom the
17 : * Software is furnished to do so, subject to the following conditions:
18 : *
19 : * The above copyright notice and this permission notice shall be included
20 : * in all copies or substantial portions of the Software.
21 : *
22 : * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
23 : * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
24 : * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
25 : * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
26 : * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
27 : * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
28 : * DEALINGS IN THE SOFTWARE.
29 : ****************************************************************************/
30 :
31 : #include "rawdataset.h"
32 : #include "ogr_spatialref.h"
33 :
34 : CPL_CVSID("$Id: ace2dataset.cpp 23772 2012-01-21 08:31:07Z rouault $");
35 :
36 : CPL_C_START
37 : void GDALRegister_ACE2(void);
38 : CPL_C_END
39 :
40 : static const char * const apszCategorySource[] =
41 : {
42 : "Pure SRTM (above 60deg N pure GLOBE data, below 60S pure ACE [original] data)",
43 : "SRTM voids filled by interpolation and/or altimeter data",
44 : "SRTM data warped using the ERS-1 Geodetic Mission",
45 : "SRTM data warped using EnviSat & ERS-2 data",
46 : "Mean lake level data derived from Altimetry",
47 : "GLOBE/ACE data warped using combined altimetry (only above 60deg N)",
48 : "Pure altimetry data (derived from ERS-1 Geodetic Mission, ERS-2 and EnviSat data using Delaunay Triangulation",
49 : NULL
50 : };
51 :
52 : static const char * const apszCategoryQuality[] =
53 : {
54 : "Generic - use base datasets",
55 : "Accuracy of greater than +/- 16m",
56 : "Accuracy between +/- 16m - +/- 10m",
57 : "Accuracy between +/-10m - +/-5m",
58 : "Accuracy between +/-5m - +/-1m",
59 : "Accuracy between +/-1m",
60 : NULL
61 : };
62 :
63 : static const char * const apszCategoryConfidence[] =
64 : {
65 : "No confidence could be derived due to lack of data",
66 : "Heights generated by interpolation",
67 : "Low confidence",
68 : "Low confidence",
69 : "Low confidence",
70 : "Medium confidence",
71 : "Medium confidence",
72 : "Medium confidence",
73 : "Medium confidence",
74 : "Medium confidence",
75 : "Medium confidence",
76 : "Medium confidence",
77 : "Medium confidence",
78 : "High confidence",
79 : "High confidence",
80 : "High confidence",
81 : "High confidence",
82 : "Inland water confidence",
83 : "Inland water confidence",
84 : "Inland water confidence",
85 : "Inland water confidence",
86 : "Inland water confidence",
87 : NULL
88 : };
89 :
90 : /************************************************************************/
91 : /* ==================================================================== */
92 : /* ACE2Dataset */
93 : /* ==================================================================== */
94 : /************************************************************************/
95 :
96 : class ACE2Dataset : public GDALPamDataset
97 1 : {
98 : friend class ACE2RasterBand;
99 :
100 : double adfGeoTransform[6];
101 :
102 : public:
103 :
104 : ACE2Dataset();
105 :
106 : virtual const char *GetProjectionRef(void);
107 : virtual CPLErr GetGeoTransform( double * );
108 :
109 : static GDALDataset *Open( GDALOpenInfo * );
110 : static int Identify( GDALOpenInfo * );
111 : };
112 :
113 : /************************************************************************/
114 : /* ==================================================================== */
115 : /* ACE2RasterBand */
116 : /* ==================================================================== */
117 : /************************************************************************/
118 :
119 : class ACE2RasterBand : public RawRasterBand
120 1 : {
121 : public:
122 : ACE2RasterBand(VSILFILE* fpRaw,
123 : GDALDataType eDataType,
124 : int nXSize, int nYSize);
125 :
126 : virtual const char *GetUnitType();
127 : virtual char **GetCategoryNames();
128 : };
129 :
130 : /************************************************************************/
131 : /* ==================================================================== */
132 : /* ACE2Dataset */
133 : /* ==================================================================== */
134 : /************************************************************************/
135 :
136 : /************************************************************************/
137 : /* ACE2Dataset() */
138 : /************************************************************************/
139 :
140 1 : ACE2Dataset::ACE2Dataset()
141 : {
142 1 : adfGeoTransform[0] = 0.0;
143 1 : adfGeoTransform[1] = 1.0;
144 1 : adfGeoTransform[2] = 0.0;
145 1 : adfGeoTransform[3] = 0.0;
146 1 : adfGeoTransform[4] = 0.0;
147 1 : adfGeoTransform[5] = 1.0;
148 1 : }
149 :
150 : /************************************************************************/
151 : /* GetGeoTransform() */
152 : /************************************************************************/
153 :
154 1 : CPLErr ACE2Dataset::GetGeoTransform( double * padfTransform )
155 :
156 : {
157 1 : memcpy( padfTransform, adfGeoTransform, sizeof(double)*6 );
158 1 : return CE_None;
159 : }
160 :
161 : /************************************************************************/
162 : /* GetProjectionRef() */
163 : /************************************************************************/
164 :
165 1 : const char *ACE2Dataset::GetProjectionRef()
166 :
167 : {
168 1 : return SRS_WKT_WGS84;
169 : }
170 :
171 : /************************************************************************/
172 : /* ACE2RasterBand() */
173 : /************************************************************************/
174 :
175 1 : ACE2RasterBand::ACE2RasterBand(VSILFILE* fpRaw,
176 : GDALDataType eDataType,
177 : int nXSize, int nYSize) :
178 : RawRasterBand( fpRaw, 0, GDALGetDataTypeSize(eDataType) / 8,
179 : nXSize * GDALGetDataTypeSize(eDataType) / 8, eDataType,
180 1 : CPL_IS_LSB, nXSize, nYSize, TRUE, TRUE)
181 : {
182 1 : }
183 :
184 : /************************************************************************/
185 : /* GetUnitType() */
186 : /************************************************************************/
187 :
188 0 : const char *ACE2RasterBand::GetUnitType()
189 : {
190 0 : if (eDataType == GDT_Float32)
191 0 : return "m";
192 : else
193 0 : return "";
194 : }
195 :
196 : /************************************************************************/
197 : /* GetCategoryNames() */
198 : /************************************************************************/
199 :
200 0 : char **ACE2RasterBand::GetCategoryNames()
201 : {
202 0 : if (eDataType == GDT_Int16)
203 : {
204 0 : const char* pszName = poDS->GetDescription();
205 0 : if (strstr(pszName, "_SOURCE_"))
206 0 : return (char**) apszCategorySource;
207 0 : else if (strstr(pszName, "_QUALITY_"))
208 0 : return (char**) apszCategoryQuality;
209 0 : else if (strstr(pszName, "_CONF_"))
210 0 : return (char**) apszCategoryConfidence;
211 : }
212 :
213 0 : return NULL;
214 : }
215 :
216 : /************************************************************************/
217 : /* Identify() */
218 : /************************************************************************/
219 :
220 11923 : int ACE2Dataset::Identify( GDALOpenInfo * poOpenInfo )
221 :
222 : {
223 11923 : if (! (EQUAL(CPLGetExtension(poOpenInfo->pszFilename), "ACE2") ||
224 : strstr(poOpenInfo->pszFilename, ".ACE2.gz") ||
225 : strstr(poOpenInfo->pszFilename, ".ace2.gz")) )
226 11922 : return FALSE;
227 :
228 1 : return TRUE;
229 : }
230 :
231 : /************************************************************************/
232 : /* Open() */
233 : /************************************************************************/
234 :
235 1830 : GDALDataset *ACE2Dataset::Open( GDALOpenInfo * poOpenInfo )
236 :
237 : {
238 1830 : if (!Identify(poOpenInfo))
239 1829 : return NULL;
240 :
241 1 : const char* pszBasename = CPLGetBasename(poOpenInfo->pszFilename);
242 1 : int nXSize = 0, nYSize = 0;
243 :
244 1 : if (strlen(pszBasename) < 7)
245 0 : return NULL;
246 :
247 : /* Determine southwest coordinates from filename */
248 :
249 : /* e.g. 30S120W_5M.ACE2 */
250 : char pszLatLonValueString[4];
251 1 : memset(pszLatLonValueString, 0, 4);
252 1 : strncpy(pszLatLonValueString, &pszBasename[0], 2);
253 1 : int southWestLat = atoi(pszLatLonValueString);
254 1 : memset(pszLatLonValueString, 0, 4);
255 1 : strncpy(pszLatLonValueString, &pszBasename[3], 3);
256 1 : int southWestLon = atoi(pszLatLonValueString);
257 :
258 1 : if(pszBasename[2] == 'N' || pszBasename[2] == 'n')
259 : /*southWestLat = southWestLat*/;
260 0 : else if(pszBasename[2] == 'S' || pszBasename[2] == 's')
261 0 : southWestLat = southWestLat * -1;
262 : else
263 0 : return NULL;
264 :
265 1 : if(pszBasename[6] == 'E' || pszBasename[6] == 'e')
266 : /*southWestLon = southWestLon*/;
267 0 : else if(pszBasename[6] == 'W' || pszBasename[6] == 'w')
268 0 : southWestLon = southWestLon * -1;
269 : else
270 0 : return NULL;
271 :
272 :
273 : GDALDataType eDT;
274 1 : if (strstr(pszBasename, "_CONF_") ||
275 : strstr(pszBasename, "_QUALITY_") ||
276 : strstr(pszBasename, "_SOURCE_"))
277 0 : eDT = GDT_Int16;
278 : else
279 1 : eDT = GDT_Float32;
280 1 : int nWordSize = GDALGetDataTypeSize(eDT) / 8;
281 :
282 : VSIStatBufL sStat;
283 1 : if (strstr(pszBasename, "_5M"))
284 1 : sStat.st_size = 180 * 180 * nWordSize;
285 0 : else if (strstr(pszBasename, "_30S"))
286 0 : sStat.st_size = 1800 * 1800 * nWordSize;
287 0 : else if (strstr(pszBasename, "_9S"))
288 0 : sStat.st_size = 6000 * 6000 * nWordSize;
289 0 : else if (strstr(pszBasename, "_3S"))
290 0 : sStat.st_size = 18000 * 18000 * nWordSize;
291 : /* Check file size otherwise */
292 0 : else if(VSIStatL(poOpenInfo->pszFilename, &sStat) != 0)
293 : {
294 0 : return NULL;
295 : }
296 :
297 1 : double dfPixelSize = 0;
298 1 : if (sStat.st_size == 180 * 180 * nWordSize)
299 : {
300 : /* 5 minute */
301 1 : nXSize = nYSize = 180;
302 1 : dfPixelSize = 5. / 60;
303 : }
304 0 : else if (sStat.st_size == 1800 * 1800 * nWordSize)
305 : {
306 : /* 30 s */
307 0 : nXSize = nYSize = 1800;
308 0 : dfPixelSize = 30. / 3600;
309 : }
310 0 : else if (sStat.st_size == 6000 * 6000 * nWordSize)
311 : {
312 : /* 9 s */
313 0 : nXSize = nYSize = 6000;
314 0 : dfPixelSize = 9. / 3600;
315 : }
316 0 : else if (sStat.st_size == 18000 * 18000 * nWordSize)
317 : {
318 : /* 3 s */
319 0 : nXSize = nYSize = 18000;
320 0 : dfPixelSize = 3. / 3600;
321 : }
322 : else
323 0 : return NULL;
324 :
325 : /* -------------------------------------------------------------------- */
326 : /* Open file. */
327 : /* -------------------------------------------------------------------- */
328 :
329 1 : CPLString osFilename = poOpenInfo->pszFilename;
330 1 : if ((strstr(poOpenInfo->pszFilename, ".ACE2.gz") ||
331 : strstr(poOpenInfo->pszFilename, ".ace2.gz")) &&
332 : strncmp(poOpenInfo->pszFilename, "/vsigzip/", 9) != 0)
333 0 : osFilename = "/vsigzip/" + osFilename;
334 :
335 1 : VSILFILE* fpImage = VSIFOpenL( osFilename, "rb+" );
336 1 : if (fpImage == NULL)
337 0 : return NULL;
338 :
339 : /* -------------------------------------------------------------------- */
340 : /* Create the dataset. */
341 : /* -------------------------------------------------------------------- */
342 : ACE2Dataset *poDS;
343 :
344 1 : poDS = new ACE2Dataset();
345 1 : poDS->nRasterXSize = nXSize;
346 1 : poDS->nRasterYSize = nYSize;
347 :
348 1 : poDS->adfGeoTransform[0] = southWestLon;
349 1 : poDS->adfGeoTransform[1] = dfPixelSize;
350 1 : poDS->adfGeoTransform[2] = 0.0;
351 1 : poDS->adfGeoTransform[3] = southWestLat + nYSize * dfPixelSize;
352 1 : poDS->adfGeoTransform[4] = 0.0;
353 1 : poDS->adfGeoTransform[5] = -dfPixelSize;
354 :
355 : /* -------------------------------------------------------------------- */
356 : /* Create band information objects */
357 : /* -------------------------------------------------------------------- */
358 2 : poDS->SetBand( 1, new ACE2RasterBand(fpImage, eDT, nXSize, nYSize ) );
359 :
360 : /* -------------------------------------------------------------------- */
361 : /* Initialize any PAM information. */
362 : /* -------------------------------------------------------------------- */
363 1 : poDS->SetDescription( poOpenInfo->pszFilename );
364 1 : poDS->TryLoadXML();
365 :
366 : /* -------------------------------------------------------------------- */
367 : /* Check for overviews. */
368 : /* -------------------------------------------------------------------- */
369 1 : poDS->oOvManager.Initialize( poDS, poOpenInfo->pszFilename );
370 :
371 1 : return( poDS );
372 : }
373 :
374 : /************************************************************************/
375 : /* GDALRegister_ACE2() */
376 : /************************************************************************/
377 :
378 582 : void GDALRegister_ACE2()
379 :
380 : {
381 : GDALDriver *poDriver;
382 :
383 582 : if( GDALGetDriverByName( "ACE2" ) == NULL )
384 : {
385 561 : poDriver = new GDALDriver();
386 :
387 561 : poDriver->SetDescription( "ACE2" );
388 : poDriver->SetMetadataItem( GDAL_DMD_LONGNAME,
389 561 : "ACE2" );
390 : poDriver->SetMetadataItem( GDAL_DMD_HELPTOPIC,
391 561 : "frmt_various.html#ACE2" );
392 561 : poDriver->SetMetadataItem( GDAL_DMD_EXTENSION, "ACE2" );
393 561 : poDriver->SetMetadataItem( GDAL_DCAP_VIRTUALIO, "YES" );
394 :
395 561 : poDriver->pfnOpen = ACE2Dataset::Open;
396 561 : poDriver->pfnIdentify = ACE2Dataset::Identify;
397 :
398 561 : GetGDALDriverManager()->RegisterDriver( poDriver );
399 : }
400 582 : }
|