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