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