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