Line data Source code
1 : /******************************************************************************
2 : *
3 : * Project: SRTM HGT Driver
4 : * Purpose: SRTM HGT File Read Support.
5 : * http://dds.cr.usgs.gov/srtm/version2_1/Documentation/SRTM_Topo.pdf
6 : * http://www2.jpl.nasa.gov/srtm/faq.html
7 : * http://dds.cr.usgs.gov/srtm/version2_1
8 : * Authors: Michael Mazzella, Even Rouault
9 : *
10 : ******************************************************************************
11 : * Copyright (c) 2005, Frank Warmerdam <warmerdam@pobox.com>
12 : * Copyright (c) 2007-2011, Even Rouault <even dot rouault at spatialys.com>
13 : *
14 : * Permission is hereby granted, free of charge, to any person obtaining a
15 : * copy of this software and associated documentation files (the "Software"),
16 : * to deal in the Software without restriction, including without limitation
17 : * the rights to use, copy, modify, merge, publish, distribute, sublicense,
18 : * and/or sell copies of the Software, and to permit persons to whom the
19 : * Software is furnished to do so, subject to the following conditions:
20 : *
21 : * The above copyright notice and this permission notice shall be included
22 : * in all copies or substantial portions of the Software.
23 : *
24 : * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
25 : * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
26 : * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
27 : * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
28 : * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
29 : * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
30 : * DEALINGS IN THE SOFTWARE.
31 : ****************************************************************************/
32 :
33 : #include "cpl_port.h"
34 : #include "cpl_string.h"
35 : #include "gdal_frmts.h"
36 : #include "gdal_pam.h"
37 : #include "ogr_spatialref.h"
38 :
39 : #include <cmath>
40 :
41 : constexpr GInt16 SRTMHG_NODATA_VALUE = -32768;
42 :
43 : /************************************************************************/
44 : /* ==================================================================== */
45 : /* SRTMHGTDataset */
46 : /* ==================================================================== */
47 : /************************************************************************/
48 :
49 : class SRTMHGTRasterBand;
50 :
51 : class SRTMHGTDataset final : public GDALPamDataset
52 : {
53 : friend class SRTMHGTRasterBand;
54 :
55 : VSILFILE *fpImage = nullptr;
56 : double adfGeoTransform[6];
57 : GByte *pabyBuffer = nullptr;
58 : OGRSpatialReference m_oSRS{};
59 :
60 : public:
61 : SRTMHGTDataset();
62 : virtual ~SRTMHGTDataset();
63 :
64 : const OGRSpatialReference *GetSpatialRef() const override;
65 : virtual CPLErr GetGeoTransform(double *) override;
66 :
67 : static int Identify(GDALOpenInfo *poOpenInfo);
68 : static GDALDataset *Open(GDALOpenInfo *);
69 : static GDALDataset *CreateCopy(const char *pszFilename,
70 : GDALDataset *poSrcDS, int bStrict,
71 : char **papszOptions,
72 : GDALProgressFunc pfnProgress,
73 : void *pProgressData);
74 : };
75 :
76 : /************************************************************************/
77 : /* ==================================================================== */
78 : /* SRTMHGTRasterBand */
79 : /* ==================================================================== */
80 : /************************************************************************/
81 :
82 : class SRTMHGTRasterBand final : public GDALPamRasterBand
83 : {
84 : friend class SRTMHGTDataset;
85 :
86 : int bNoDataSet;
87 : double dfNoDataValue;
88 :
89 : public:
90 : SRTMHGTRasterBand(SRTMHGTDataset *, int, GDALDataType);
91 :
92 : virtual CPLErr IReadBlock(int, int, void *) override;
93 : virtual CPLErr IWriteBlock(int nBlockXOff, int nBlockYOff,
94 : void *pImage) override;
95 :
96 : virtual GDALColorInterp GetColorInterpretation() override;
97 :
98 : virtual double GetNoDataValue(int *pbSuccess = nullptr) override;
99 :
100 : virtual const char *GetUnitType() override;
101 : };
102 :
103 : /************************************************************************/
104 : /* SRTMHGTRasterBand() */
105 : /************************************************************************/
106 :
107 12 : SRTMHGTRasterBand::SRTMHGTRasterBand(SRTMHGTDataset *poDSIn, int nBandIn,
108 12 : GDALDataType eDT)
109 12 : : bNoDataSet(TRUE), dfNoDataValue(SRTMHG_NODATA_VALUE)
110 : {
111 12 : poDS = poDSIn;
112 12 : nBand = nBandIn;
113 12 : eDataType = eDT;
114 12 : nBlockXSize = poDSIn->nRasterXSize;
115 12 : nBlockYSize = 1;
116 12 : }
117 :
118 : /************************************************************************/
119 : /* IReadBlock() */
120 : /************************************************************************/
121 :
122 14408 : CPLErr SRTMHGTRasterBand::IReadBlock(int /*nBlockXOff*/, int nBlockYOff,
123 : void *pImage)
124 : {
125 14408 : SRTMHGTDataset *poGDS = reinterpret_cast<SRTMHGTDataset *>(poDS);
126 :
127 : /* -------------------------------------------------------------------- */
128 : /* Load the desired data into the working buffer. */
129 : /* -------------------------------------------------------------------- */
130 14408 : const int nDTSize = GDALGetDataTypeSizeBytes(eDataType);
131 14408 : VSIFSeekL(poGDS->fpImage,
132 14408 : static_cast<size_t>(nBlockYOff) * nBlockXSize * nDTSize,
133 : SEEK_SET);
134 14408 : VSIFReadL((unsigned char *)pImage, nBlockXSize, nDTSize, poGDS->fpImage);
135 : #ifdef CPL_LSB
136 14408 : GDALSwapWords(pImage, nDTSize, nBlockXSize, nDTSize);
137 : #endif
138 :
139 14408 : return CE_None;
140 : }
141 :
142 : /************************************************************************/
143 : /* IWriteBlock() */
144 : /************************************************************************/
145 :
146 1201 : CPLErr SRTMHGTRasterBand::IWriteBlock(int /*nBlockXOff*/, int nBlockYOff,
147 : void *pImage)
148 : {
149 1201 : SRTMHGTDataset *poGDS = reinterpret_cast<SRTMHGTDataset *>(poDS);
150 :
151 1201 : if (poGDS->eAccess != GA_Update)
152 0 : return CE_Failure;
153 :
154 1201 : const int nDTSize = GDALGetDataTypeSizeBytes(eDataType);
155 1201 : VSIFSeekL(poGDS->fpImage,
156 1201 : static_cast<size_t>(nBlockYOff) * nBlockXSize * nDTSize,
157 : SEEK_SET);
158 :
159 : #ifdef CPL_LSB
160 1201 : if (nDTSize > 1)
161 : {
162 1201 : memcpy(poGDS->pabyBuffer, pImage,
163 1201 : static_cast<size_t>(nBlockXSize) * nDTSize);
164 1201 : GDALSwapWords(poGDS->pabyBuffer, nDTSize, nBlockXSize, nDTSize);
165 1201 : VSIFWriteL(reinterpret_cast<unsigned char *>(poGDS->pabyBuffer),
166 1201 : nBlockXSize, nDTSize, poGDS->fpImage);
167 : }
168 : else
169 : #endif
170 : {
171 0 : VSIFWriteL(reinterpret_cast<unsigned char *>(pImage), nBlockXSize,
172 : nDTSize, poGDS->fpImage);
173 : }
174 :
175 1201 : return CE_None;
176 : }
177 :
178 : /************************************************************************/
179 : /* GetNoDataValue() */
180 : /************************************************************************/
181 :
182 13 : double SRTMHGTRasterBand::GetNoDataValue(int *pbSuccess)
183 :
184 : {
185 13 : if (eDataType == GDT_Byte)
186 0 : return GDALPamRasterBand::GetNoDataValue(pbSuccess);
187 :
188 13 : if (pbSuccess)
189 10 : *pbSuccess = bNoDataSet;
190 :
191 13 : return dfNoDataValue;
192 : }
193 :
194 : /************************************************************************/
195 : /* GetUnitType() */
196 : /************************************************************************/
197 :
198 6 : const char *SRTMHGTRasterBand::GetUnitType()
199 : {
200 6 : const char *pszExt = CPLGetExtension(poDS->GetDescription());
201 6 : if (EQUAL(pszExt, "err") || EQUAL(pszExt, "img") || EQUAL(pszExt, "num") ||
202 6 : EQUAL(pszExt, "swb"))
203 : {
204 0 : return "";
205 : }
206 6 : return "m";
207 : }
208 :
209 : /************************************************************************/
210 : /* GetColorInterpretation() */
211 : /************************************************************************/
212 :
213 5 : GDALColorInterp SRTMHGTRasterBand::GetColorInterpretation()
214 : {
215 5 : return GCI_Undefined;
216 : }
217 :
218 : /************************************************************************/
219 : /* ==================================================================== */
220 : /* SRTMHGTDataset */
221 : /* ==================================================================== */
222 : /************************************************************************/
223 :
224 : /************************************************************************/
225 : /* SRTMHGTDataset() */
226 : /************************************************************************/
227 :
228 12 : SRTMHGTDataset::SRTMHGTDataset()
229 : {
230 12 : m_oSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
231 12 : if (CPLTestBool(CPLGetConfigOption("REPORT_COMPD_CS", "NO")))
232 : {
233 0 : m_oSRS.importFromWkt(
234 : "COMPD_CS[\"WGS 84 + EGM96 geoid height\", GEOGCS[\"WGS 84\", "
235 : "DATUM[\"WGS_1984\", SPHEROID[\"WGS 84\",6378137,298.257223563, "
236 : "AUTHORITY[\"EPSG\",\"7030\"]], AUTHORITY[\"EPSG\",\"6326\"]], "
237 : "PRIMEM[\"Greenwich\",0, AUTHORITY[\"EPSG\",\"8901\"]], "
238 : "UNIT[\"degree\",0.0174532925199433, "
239 : "AUTHORITY[\"EPSG\",\"9122\"]], AUTHORITY[\"EPSG\",\"4326\"]], "
240 : "VERT_CS[\"EGM96 geoid height\", VERT_DATUM[\"EGM96 geoid\",2005, "
241 : "AUTHORITY[\"EPSG\",\"5171\"]], UNIT[\"metre\",1, "
242 : "AUTHORITY[\"EPSG\",\"9001\"]], AXIS[\"Up\",UP], "
243 : "AUTHORITY[\"EPSG\",\"5773\"]]]");
244 : }
245 : else
246 : {
247 12 : m_oSRS.importFromWkt(SRS_WKT_WGS84_LAT_LONG);
248 : }
249 12 : adfGeoTransform[0] = 0.0;
250 12 : adfGeoTransform[1] = 1.0;
251 12 : adfGeoTransform[2] = 0.0;
252 12 : adfGeoTransform[3] = 0.0;
253 12 : adfGeoTransform[4] = 0.0;
254 12 : adfGeoTransform[5] = 1.0;
255 12 : }
256 :
257 : /************************************************************************/
258 : /* ~SRTMHGTDataset() */
259 : /************************************************************************/
260 :
261 24 : SRTMHGTDataset::~SRTMHGTDataset()
262 : {
263 12 : FlushCache(true);
264 12 : if (fpImage != nullptr)
265 12 : VSIFCloseL(fpImage);
266 12 : CPLFree(pabyBuffer);
267 24 : }
268 :
269 : /************************************************************************/
270 : /* GetGeoTransform() */
271 : /************************************************************************/
272 :
273 9 : CPLErr SRTMHGTDataset::GetGeoTransform(double *padfTransform)
274 : {
275 9 : memcpy(padfTransform, adfGeoTransform, sizeof(double) * 6);
276 9 : return CE_None;
277 : }
278 :
279 : /************************************************************************/
280 : /* GetSpatialRef() */
281 : /************************************************************************/
282 :
283 9 : const OGRSpatialReference *SRTMHGTDataset::GetSpatialRef() const
284 :
285 : {
286 9 : return &m_oSRS;
287 : }
288 :
289 : /************************************************************************/
290 : /* Identify() */
291 : /************************************************************************/
292 :
293 52694 : int SRTMHGTDataset::Identify(GDALOpenInfo *poOpenInfo)
294 :
295 : {
296 52694 : const char *fileName = CPLGetFilename(poOpenInfo->pszFilename);
297 52694 : if (strlen(fileName) < 11 || fileName[7] != '.')
298 49439 : return FALSE;
299 6510 : CPLString osLCFilename(CPLString(fileName).tolower());
300 3703 : if ((osLCFilename[0] != 'n' && osLCFilename[0] != 's') ||
301 448 : (osLCFilename[3] != 'e' && osLCFilename[3] != 'w'))
302 3156 : return FALSE;
303 99 : if (!STARTS_WITH(fileName, "/vsizip/") && osLCFilename.endsWith(".hgt.zip"))
304 : {
305 4 : CPLString osNewName("/vsizip/");
306 2 : osNewName += poOpenInfo->pszFilename;
307 2 : osNewName += "/";
308 2 : osNewName += CPLString(fileName).substr(0, 7);
309 2 : osNewName += ".hgt";
310 4 : GDALOpenInfo oOpenInfo(osNewName, GA_ReadOnly);
311 2 : return Identify(&oOpenInfo);
312 : }
313 :
314 194 : if (!STARTS_WITH(fileName, "/vsizip/") &&
315 194 : osLCFilename.endsWith(".srtmswbd.raw.zip"))
316 : {
317 4 : CPLString osNewName("/vsizip/");
318 2 : osNewName += poOpenInfo->pszFilename;
319 2 : osNewName += "/";
320 2 : osNewName += CPLString(fileName).substr(0, 7);
321 2 : osNewName += ".raw";
322 4 : GDALOpenInfo oOpenInfo(osNewName, GA_ReadOnly);
323 2 : return Identify(&oOpenInfo);
324 : }
325 :
326 : // .hgts and .err files from
327 : // https://e4ftl01.cr.usgs.gov/MEASURES/NASADEM_SHHP.001/2000.02.11/ .img
328 : // and .img.num files from
329 : // https://e4ftl01.cr.usgs.gov/MEASURES/NASADEM_SIM.001/2000.02.11/
330 195 : if (!osLCFilename.endsWith(".hgt") && !osLCFilename.endsWith(".hgts") &&
331 98 : !osLCFilename.endsWith(".err") && !osLCFilename.endsWith(".img") &&
332 98 : !osLCFilename.endsWith(".num") && // .img.num or .num
333 98 : !osLCFilename.endsWith(".raw") &&
334 95 : !osLCFilename.endsWith(
335 195 : ".swb") && // https://e4ftl01.cr.usgs.gov/MEASURES/NASADEM_HGT.001/2000.02.11/
336 95 : !osLCFilename.endsWith(".hgt.gz"))
337 0 : return FALSE;
338 :
339 : /* -------------------------------------------------------------------- */
340 : /* We check the file size to see if it is */
341 : /* SRTM1 (below or above lat 50) or SRTM 3 */
342 : /* -------------------------------------------------------------------- */
343 : VSIStatBufL fileStat;
344 :
345 95 : if (VSIStatL(poOpenInfo->pszFilename, &fileStat) != 0)
346 69 : return FALSE;
347 26 : if (fileStat.st_size != 3601 * 3601 &&
348 23 : fileStat.st_size != 3601 * 3601 * 2 &&
349 23 : fileStat.st_size != 3601 * 3601 * 4 && // .hgts
350 21 : fileStat.st_size != 1801 * 3601 * 2 &&
351 21 : fileStat.st_size != 1201 * 1201 * 2)
352 0 : return FALSE;
353 :
354 26 : return TRUE;
355 : }
356 :
357 : /************************************************************************/
358 : /* Open() */
359 : /************************************************************************/
360 :
361 14 : GDALDataset *SRTMHGTDataset::Open(GDALOpenInfo *poOpenInfo)
362 : {
363 14 : if (!Identify(poOpenInfo))
364 0 : return nullptr;
365 :
366 14 : const char *fileName = CPLGetFilename(poOpenInfo->pszFilename);
367 28 : CPLString osLCFilename(CPLString(fileName).tolower());
368 14 : if (!STARTS_WITH(fileName, "/vsizip/") && osLCFilename.endsWith(".hgt.zip"))
369 : {
370 2 : CPLString osFilename("/vsizip/");
371 1 : osFilename += poOpenInfo->pszFilename;
372 1 : osFilename += "/";
373 1 : osFilename += CPLString(fileName).substr(0, 7);
374 1 : osFilename += ".hgt";
375 1 : GDALOpenInfo oOpenInfo(osFilename, poOpenInfo->eAccess);
376 1 : GDALDataset *poDS = Open(&oOpenInfo);
377 1 : if (poDS != nullptr)
378 : {
379 : // override description with the main one
380 1 : poDS->SetDescription(poOpenInfo->pszFilename);
381 : }
382 1 : return poDS;
383 : }
384 :
385 26 : if (!STARTS_WITH(fileName, "/vsizip/") &&
386 26 : osLCFilename.endsWith(".srtmswbd.raw.zip"))
387 : {
388 2 : CPLString osFilename("/vsizip/");
389 1 : osFilename += poOpenInfo->pszFilename;
390 1 : osFilename += "/";
391 1 : osFilename += CPLString(fileName).substr(0, 7);
392 1 : osFilename += ".raw";
393 1 : GDALOpenInfo oOpenInfo(osFilename, poOpenInfo->eAccess);
394 1 : GDALDataset *poDS = Open(&oOpenInfo);
395 1 : if (poDS != nullptr)
396 : {
397 : // override description with the main one
398 1 : poDS->SetDescription(poOpenInfo->pszFilename);
399 : }
400 1 : return poDS;
401 : }
402 :
403 : char latLonValueString[4];
404 12 : memset(latLonValueString, 0, 4);
405 12 : strncpy(latLonValueString, &fileName[1], 2);
406 12 : int southWestLat = atoi(latLonValueString);
407 12 : memset(latLonValueString, 0, 4);
408 : // cppcheck-suppress redundantCopy
409 12 : strncpy(latLonValueString, &fileName[4], 3);
410 12 : int southWestLon = atoi(latLonValueString);
411 :
412 12 : if (fileName[0] == 'N' || fileName[0] == 'n')
413 : /*southWestLat = southWestLat */;
414 0 : else if (fileName[0] == 'S' || fileName[0] == 's')
415 0 : southWestLat = southWestLat * -1;
416 : else
417 0 : return nullptr;
418 :
419 12 : if (fileName[3] == 'E' || fileName[3] == 'e')
420 : /*southWestLon = southWestLon */;
421 11 : else if (fileName[3] == 'W' || fileName[3] == 'w')
422 11 : southWestLon = southWestLon * -1;
423 : else
424 0 : return nullptr;
425 :
426 : /* -------------------------------------------------------------------- */
427 : /* Create a corresponding GDALDataset. */
428 : /* -------------------------------------------------------------------- */
429 12 : SRTMHGTDataset *poDS = new SRTMHGTDataset();
430 :
431 12 : poDS->fpImage = poOpenInfo->fpL;
432 12 : poOpenInfo->fpL = nullptr;
433 :
434 : VSIStatBufL fileStat;
435 12 : if (VSIStatL(poOpenInfo->pszFilename, &fileStat) != 0)
436 : {
437 0 : delete poDS;
438 0 : return nullptr;
439 : }
440 :
441 : int numPixels_x, numPixels_y;
442 :
443 12 : GDALDataType eDT = GDT_Int16;
444 12 : switch (fileStat.st_size)
445 : {
446 1 : case 3601 * 3601:
447 1 : numPixels_x = numPixels_y = 3601;
448 1 : eDT = GDT_Byte;
449 1 : break;
450 0 : case 3601 * 3601 * 2:
451 0 : numPixels_x = numPixels_y = 3601;
452 0 : break;
453 1 : case 3601 * 3601 * 4: // .hgts
454 1 : numPixels_x = numPixels_y = 3601;
455 1 : eDT = GDT_Float32;
456 1 : break;
457 0 : case 1801 * 3601 * 2:
458 0 : numPixels_x = 1801;
459 0 : numPixels_y = 3601;
460 0 : break;
461 10 : case 1201 * 1201 * 2:
462 10 : numPixels_x = numPixels_y = 1201;
463 10 : break;
464 0 : default:
465 0 : numPixels_x = numPixels_y = 0;
466 0 : break;
467 : }
468 :
469 12 : poDS->eAccess = poOpenInfo->eAccess;
470 : #ifdef CPL_LSB
471 12 : if (poDS->eAccess == GA_Update && eDT != GDT_Byte)
472 : {
473 1 : poDS->pabyBuffer =
474 1 : static_cast<GByte *>(CPLMalloc(numPixels_x * sizeof(eDT)));
475 : }
476 : #endif
477 :
478 : /* -------------------------------------------------------------------- */
479 : /* Capture some information from the file that is of interest. */
480 : /* -------------------------------------------------------------------- */
481 12 : poDS->nRasterXSize = numPixels_x;
482 12 : poDS->nRasterYSize = numPixels_y;
483 12 : poDS->nBands = 1;
484 :
485 12 : poDS->adfGeoTransform[0] = southWestLon - 0.5 / (numPixels_x - 1);
486 12 : poDS->adfGeoTransform[1] = 1.0 / (numPixels_x - 1);
487 12 : poDS->adfGeoTransform[2] = 0.0;
488 12 : poDS->adfGeoTransform[3] = southWestLat + 1 + 0.5 / (numPixels_y - 1);
489 12 : poDS->adfGeoTransform[4] = 0.0;
490 12 : poDS->adfGeoTransform[5] = -1.0 / (numPixels_y - 1);
491 :
492 12 : poDS->SetMetadataItem(GDALMD_AREA_OR_POINT, GDALMD_AOP_POINT);
493 :
494 : /* -------------------------------------------------------------------- */
495 : /* Create band information object. */
496 : /* -------------------------------------------------------------------- */
497 12 : SRTMHGTRasterBand *tmpBand = new SRTMHGTRasterBand(poDS, 1, eDT);
498 12 : poDS->SetBand(1, tmpBand);
499 :
500 : /* -------------------------------------------------------------------- */
501 : /* Initialize any PAM information. */
502 : /* -------------------------------------------------------------------- */
503 12 : poDS->SetDescription(poOpenInfo->pszFilename);
504 12 : poDS->TryLoadXML();
505 :
506 : /* -------------------------------------------------------------------- */
507 : /* Support overviews. */
508 : /* -------------------------------------------------------------------- */
509 12 : poDS->oOvManager.Initialize(poDS, poOpenInfo->pszFilename);
510 :
511 12 : return poDS;
512 : }
513 :
514 : /************************************************************************/
515 : /* CreateCopy() */
516 : /************************************************************************/
517 :
518 23 : GDALDataset *SRTMHGTDataset::CreateCopy(const char *pszFilename,
519 : GDALDataset *poSrcDS, int bStrict,
520 : char ** /* papszOptions*/,
521 : GDALProgressFunc pfnProgress,
522 : void *pProgressData)
523 : {
524 : /* -------------------------------------------------------------------- */
525 : /* Some some rudimentary checks */
526 : /* -------------------------------------------------------------------- */
527 23 : const int nBands = poSrcDS->GetRasterCount();
528 23 : if (nBands == 0)
529 : {
530 1 : CPLError(
531 : CE_Failure, CPLE_NotSupported,
532 : "SRTMHGT driver does not support source dataset with zero band.\n");
533 1 : return nullptr;
534 : }
535 22 : else if (nBands != 1)
536 : {
537 4 : CPLError((bStrict) ? CE_Failure : CE_Warning, CPLE_NotSupported,
538 : "SRTMHGT driver only uses the first band of the dataset.\n");
539 4 : if (bStrict)
540 4 : return nullptr;
541 : }
542 :
543 : /* -------------------------------------------------------------------- */
544 : /* Checks the input SRS */
545 : /* -------------------------------------------------------------------- */
546 36 : OGRSpatialReference ogrsr_input;
547 18 : ogrsr_input.importFromWkt(poSrcDS->GetProjectionRef());
548 :
549 36 : OGRSpatialReference ogrsr_wgs84;
550 18 : ogrsr_wgs84.SetWellKnownGeogCS("WGS84");
551 :
552 18 : if (ogrsr_input.IsSameGeogCS(&ogrsr_wgs84) == FALSE)
553 : {
554 0 : CPLError(CE_Warning, CPLE_AppDefined,
555 : "The source projection coordinate system is %s. Only WGS 84 "
556 : "is supported.\nThe SRTMHGT driver will generate a file as "
557 : "if the source was WGS 84 projection coordinate system.",
558 : poSrcDS->GetProjectionRef());
559 : }
560 :
561 : /* -------------------------------------------------------------------- */
562 : /* Work out the LL origin. */
563 : /* -------------------------------------------------------------------- */
564 : double adfGeoTransform[6];
565 18 : if (poSrcDS->GetGeoTransform(adfGeoTransform) != CE_None)
566 : {
567 0 : CPLError(CE_Failure, CPLE_AppDefined,
568 : "Source image must have a geo transform matrix.");
569 0 : return nullptr;
570 : }
571 :
572 : const int nLLOriginLat = static_cast<int>(
573 36 : std::floor(adfGeoTransform[3] +
574 18 : poSrcDS->GetRasterYSize() * adfGeoTransform[5] + 0.5));
575 :
576 18 : int nLLOriginLong = static_cast<int>(std::floor(adfGeoTransform[0] + 0.5));
577 :
578 54 : if (std::abs(nLLOriginLat -
579 18 : (adfGeoTransform[3] + (poSrcDS->GetRasterYSize() - 0.5) *
580 23 : adfGeoTransform[5])) > 1e-10 ||
581 5 : std::abs(nLLOriginLong -
582 5 : (adfGeoTransform[0] + 0.5 * adfGeoTransform[1])) > 1e-10)
583 : {
584 13 : CPLError(CE_Warning, CPLE_AppDefined,
585 : "The corner coordinates of the source are not properly "
586 : "aligned on plain latitude/longitude boundaries.");
587 : }
588 :
589 : /* -------------------------------------------------------------------- */
590 : /* Check image dimensions. */
591 : /* -------------------------------------------------------------------- */
592 18 : const int nXSize = poSrcDS->GetRasterXSize();
593 18 : const int nYSize = poSrcDS->GetRasterYSize();
594 :
595 18 : if (!((nXSize == 1201 && nYSize == 1201) ||
596 0 : (nXSize == 3601 && nYSize == 3601) ||
597 0 : (nXSize == 1801 && nYSize == 3601)))
598 : {
599 11 : CPLError(
600 : CE_Failure, CPLE_AppDefined,
601 : "Image dimensions should be 1201x1201, 3601x3601 or 1801x3601.");
602 11 : return nullptr;
603 : }
604 :
605 : /* -------------------------------------------------------------------- */
606 : /* Check filename. */
607 : /* -------------------------------------------------------------------- */
608 : char expectedFileName[12];
609 :
610 14 : CPLsnprintf(expectedFileName, sizeof(expectedFileName), "%c%02d%c%03d.HGT",
611 : (nLLOriginLat >= 0) ? 'N' : 'S',
612 7 : (nLLOriginLat >= 0) ? nLLOriginLat : -nLLOriginLat,
613 : (nLLOriginLong >= 0) ? 'E' : 'W',
614 : (nLLOriginLong >= 0) ? nLLOriginLong : -nLLOriginLong);
615 :
616 7 : if (!EQUAL(expectedFileName, CPLGetFilename(pszFilename)))
617 : {
618 0 : CPLError(CE_Warning, CPLE_AppDefined, "Expected output filename is %s.",
619 : expectedFileName);
620 : }
621 :
622 : /* -------------------------------------------------------------------- */
623 : /* Write output file. */
624 : /* -------------------------------------------------------------------- */
625 7 : VSILFILE *fp = VSIFOpenL(pszFilename, "wb");
626 7 : if (fp == nullptr)
627 : {
628 2 : CPLError(CE_Failure, CPLE_FileIO, "Cannot create file %s", pszFilename);
629 2 : return nullptr;
630 : }
631 :
632 : GInt16 *panData =
633 5 : reinterpret_cast<GInt16 *>(CPLMalloc(sizeof(GInt16) * nXSize));
634 5 : GDALRasterBand *poSrcBand = poSrcDS->GetRasterBand(1);
635 :
636 : int bSrcBandHasNoData;
637 5 : double srcBandNoData = poSrcBand->GetNoDataValue(&bSrcBandHasNoData);
638 :
639 6010 : for (int iY = 0; iY < nYSize; iY++)
640 : {
641 6005 : if (poSrcBand->RasterIO(GF_Read, 0, iY, nXSize, 1,
642 : reinterpret_cast<void *>(panData), nXSize, 1,
643 6005 : GDT_Int16, 0, 0, nullptr) != CE_None)
644 : {
645 0 : VSIFCloseL(fp);
646 0 : CPLFree(panData);
647 0 : return nullptr;
648 : }
649 :
650 : /* Translate nodata values */
651 6005 : if (bSrcBandHasNoData && srcBandNoData != SRTMHG_NODATA_VALUE)
652 : {
653 0 : for (int iX = 0; iX < nXSize; iX++)
654 : {
655 0 : if (panData[iX] == srcBandNoData)
656 0 : panData[iX] = SRTMHG_NODATA_VALUE;
657 : }
658 : }
659 :
660 : #ifdef CPL_LSB
661 6005 : GDALSwapWords(panData, 2, nXSize, 2);
662 : #endif
663 :
664 6005 : if (VSIFWriteL(panData, sizeof(GInt16) * nXSize, 1, fp) != 1)
665 : {
666 0 : CPLError(CE_Failure, CPLE_FileIO,
667 : "Failed to write line %d in SRTMHGT dataset.\n", iY);
668 0 : VSIFCloseL(fp);
669 0 : CPLFree(panData);
670 0 : return nullptr;
671 : }
672 :
673 6005 : if (pfnProgress && !pfnProgress((iY + 1) / static_cast<double>(nYSize),
674 : nullptr, pProgressData))
675 : {
676 0 : CPLError(CE_Failure, CPLE_UserInterrupt,
677 : "User terminated CreateCopy()");
678 0 : VSIFCloseL(fp);
679 0 : CPLFree(panData);
680 0 : return nullptr;
681 : }
682 : }
683 :
684 5 : CPLFree(panData);
685 5 : VSIFCloseL(fp);
686 :
687 : /* -------------------------------------------------------------------- */
688 : /* Reopen and copy missing information into a PAM file. */
689 : /* -------------------------------------------------------------------- */
690 : GDALPamDataset *poDS =
691 5 : reinterpret_cast<GDALPamDataset *>(GDALOpen(pszFilename, GA_ReadOnly));
692 :
693 5 : if (poDS)
694 5 : poDS->CloneInfo(poSrcDS, GCIF_PAM_DEFAULT);
695 :
696 5 : return poDS;
697 : }
698 :
699 : /************************************************************************/
700 : /* GDALRegister_SRTMHGT() */
701 : /************************************************************************/
702 1520 : void GDALRegister_SRTMHGT()
703 : {
704 1520 : if (GDALGetDriverByName("SRTMHGT") != nullptr)
705 301 : return;
706 :
707 1219 : GDALDriver *poDriver = new GDALDriver();
708 :
709 1219 : poDriver->SetDescription("SRTMHGT");
710 1219 : poDriver->SetMetadataItem(GDAL_DCAP_RASTER, "YES");
711 1219 : poDriver->SetMetadataItem(GDAL_DMD_LONGNAME, "SRTMHGT File Format");
712 1219 : poDriver->SetMetadataItem(GDAL_DMD_EXTENSION, "hgt");
713 1219 : poDriver->SetMetadataItem(GDAL_DMD_HELPTOPIC,
714 1219 : "drivers/raster/srtmhgt.html");
715 1219 : poDriver->SetMetadataItem(GDAL_DMD_CREATIONDATATYPES, "Byte Int16 UInt16");
716 1219 : poDriver->SetMetadataItem(GDAL_DCAP_VIRTUALIO, "YES");
717 :
718 1219 : poDriver->pfnIdentify = SRTMHGTDataset::Identify;
719 1219 : poDriver->pfnOpen = SRTMHGTDataset::Open;
720 1219 : poDriver->pfnCreateCopy = SRTMHGTDataset::CreateCopy;
721 :
722 1219 : GetGDALDriverManager()->RegisterDriver(poDriver);
723 : }
|
