Line data Source code
1 : /******************************************************************************
2 : *
3 : * Project: JDEM Reader
4 : * Purpose: All code for Japanese DEM Reader
5 : * Author: Frank Warmerdam, warmerdam@pobox.com
6 : *
7 : ******************************************************************************
8 : * Copyright (c) 2000, Frank Warmerdam <warmerdam@pobox.com>
9 : * Copyright (c) 2009-2012, Even Rouault <even dot rouault at spatialys.com>
10 : *
11 : * SPDX-License-Identifier: MIT
12 : ****************************************************************************/
13 :
14 : #include "cpl_port.h"
15 : #include "gdal_frmts.h"
16 : #include "gdal_pam.h"
17 : #include "gdal_driver.h"
18 : #include "gdal_drivermanager.h"
19 : #include "gdal_openinfo.h"
20 : #include "gdal_cpp_functions.h"
21 :
22 : #include <algorithm>
23 :
24 : constexpr int HEADER_SIZE = 1011;
25 :
26 : /************************************************************************/
27 : /* JDEMGetField() */
28 : /************************************************************************/
29 :
30 26 : static int JDEMGetField(const char *pszField, int nWidth)
31 :
32 : {
33 26 : char szWork[32] = {};
34 26 : CPLAssert(nWidth < static_cast<int>(sizeof(szWork)));
35 :
36 26 : strncpy(szWork, pszField, nWidth);
37 26 : szWork[nWidth] = '\0';
38 :
39 52 : return atoi(szWork);
40 : }
41 :
42 : /************************************************************************/
43 : /* JDEMGetAngle() */
44 : /************************************************************************/
45 :
46 16 : static double JDEMGetAngle(const char *pszField)
47 :
48 : {
49 16 : const int nAngle = JDEMGetField(pszField, 7);
50 :
51 : // Note, this isn't very general purpose, but it would appear
52 : // from the field widths that angles are never negative. Nice
53 : // to be a country in the "first quadrant".
54 :
55 16 : const int nDegree = nAngle / 10000;
56 16 : const int nMin = (nAngle / 100) % 100;
57 16 : const int nSec = nAngle % 100;
58 :
59 16 : return nDegree + nMin / 60.0 + nSec / 3600.0;
60 : }
61 :
62 : /************************************************************************/
63 : /* ==================================================================== */
64 : /* JDEMDataset */
65 : /* ==================================================================== */
66 : /************************************************************************/
67 :
68 : class JDEMRasterBand;
69 :
70 : class JDEMDataset final : public GDALPamDataset
71 : {
72 : friend class JDEMRasterBand;
73 :
74 : VSILFILE *m_fp = nullptr;
75 : GByte m_abyHeader[HEADER_SIZE];
76 : OGRSpatialReference m_oSRS{};
77 :
78 : public:
79 : JDEMDataset();
80 : ~JDEMDataset() override;
81 :
82 : static GDALDataset *Open(GDALOpenInfo *);
83 : static int Identify(GDALOpenInfo *);
84 :
85 : CPLErr GetGeoTransform(GDALGeoTransform >) const override;
86 : const OGRSpatialReference *GetSpatialRef() const override;
87 : };
88 :
89 : /************************************************************************/
90 : /* ==================================================================== */
91 : /* JDEMRasterBand */
92 : /* ==================================================================== */
93 : /************************************************************************/
94 :
95 : class JDEMRasterBand final : public GDALPamRasterBand
96 : {
97 : friend class JDEMDataset;
98 :
99 : int m_nRecordSize = 0;
100 : char *m_pszRecord = nullptr;
101 : bool m_bBufferAllocFailed = false;
102 :
103 : public:
104 : JDEMRasterBand(JDEMDataset *, int);
105 : ~JDEMRasterBand() override;
106 :
107 : CPLErr IReadBlock(int, int, void *) override;
108 : };
109 :
110 : /************************************************************************/
111 : /* JDEMRasterBand() */
112 : /************************************************************************/
113 :
114 2 : JDEMRasterBand::JDEMRasterBand(JDEMDataset *poDSIn, int nBandIn)
115 : : // Cannot overflow as nBlockXSize <= 999.
116 2 : m_nRecordSize(poDSIn->GetRasterXSize() * 5 + 9 + 2)
117 : {
118 2 : poDS = poDSIn;
119 2 : nBand = nBandIn;
120 :
121 2 : eDataType = GDT_Float32;
122 :
123 2 : nBlockXSize = poDS->GetRasterXSize();
124 2 : nBlockYSize = 1;
125 2 : }
126 :
127 : /************************************************************************/
128 : /* ~JDEMRasterBand() */
129 : /************************************************************************/
130 :
131 4 : JDEMRasterBand::~JDEMRasterBand()
132 : {
133 2 : VSIFree(m_pszRecord);
134 4 : }
135 :
136 : /************************************************************************/
137 : /* IReadBlock() */
138 : /************************************************************************/
139 :
140 2 : CPLErr JDEMRasterBand::IReadBlock(int /* nBlockXOff */, int nBlockYOff,
141 : void *pImage)
142 :
143 : {
144 2 : JDEMDataset *poGDS = cpl::down_cast<JDEMDataset *>(poDS);
145 :
146 2 : if (m_pszRecord == nullptr)
147 : {
148 1 : if (m_bBufferAllocFailed)
149 0 : return CE_Failure;
150 :
151 1 : m_pszRecord = static_cast<char *>(VSI_MALLOC_VERBOSE(m_nRecordSize));
152 1 : if (m_pszRecord == nullptr)
153 : {
154 0 : m_bBufferAllocFailed = true;
155 0 : return CE_Failure;
156 : }
157 : }
158 :
159 2 : CPL_IGNORE_RET_VAL(
160 2 : VSIFSeekL(poGDS->m_fp, 1011 + m_nRecordSize * nBlockYOff, SEEK_SET));
161 :
162 2 : if (VSIFReadL(m_pszRecord, m_nRecordSize, 1, poGDS->m_fp) != 1)
163 : {
164 0 : CPLError(CE_Failure, CPLE_AppDefined, "Cannot read scanline %d",
165 : nBlockYOff);
166 0 : return CE_Failure;
167 : }
168 :
169 2 : if (!EQUALN(reinterpret_cast<char *>(poGDS->m_abyHeader), m_pszRecord, 6))
170 : {
171 0 : CPLError(CE_Failure, CPLE_AppDefined,
172 : "JDEM Scanline corrupt. Perhaps file was not transferred "
173 : "in binary mode?");
174 0 : return CE_Failure;
175 : }
176 :
177 2 : if (JDEMGetField(m_pszRecord + 6, 3) != nBlockYOff + 1)
178 : {
179 0 : CPLError(CE_Failure, CPLE_AppDefined,
180 : "JDEM scanline out of order, JDEM driver does not "
181 : "currently support partial datasets.");
182 0 : return CE_Failure;
183 : }
184 :
185 6 : for (int i = 0; i < nBlockXSize; i++)
186 4 : static_cast<float *>(pImage)[i] =
187 4 : JDEMGetField(m_pszRecord + 9 + 5 * i, 5) * 0.1f;
188 :
189 2 : return CE_None;
190 : }
191 :
192 : /************************************************************************/
193 : /* ==================================================================== */
194 : /* JDEMDataset */
195 : /* ==================================================================== */
196 : /************************************************************************/
197 :
198 : /************************************************************************/
199 : /* JDEMDataset() */
200 : /************************************************************************/
201 :
202 2 : JDEMDataset::JDEMDataset()
203 : {
204 2 : std::fill_n(m_abyHeader, CPL_ARRAYSIZE(m_abyHeader), static_cast<GByte>(0));
205 2 : m_oSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
206 2 : m_oSRS.importFromEPSG(4301); // Tokyo geographic CRS
207 2 : }
208 :
209 : /************************************************************************/
210 : /* ~JDEMDataset() */
211 : /************************************************************************/
212 :
213 4 : JDEMDataset::~JDEMDataset()
214 :
215 : {
216 2 : FlushCache(true);
217 2 : if (m_fp != nullptr)
218 2 : CPL_IGNORE_RET_VAL(VSIFCloseL(m_fp));
219 4 : }
220 :
221 : /************************************************************************/
222 : /* GetGeoTransform() */
223 : /************************************************************************/
224 :
225 0 : CPLErr JDEMDataset::GetGeoTransform(GDALGeoTransform >) const
226 :
227 : {
228 0 : const char *psHeader = reinterpret_cast<const char *>(m_abyHeader);
229 :
230 0 : const double dfLLLat = JDEMGetAngle(psHeader + 29);
231 0 : const double dfLLLong = JDEMGetAngle(psHeader + 36);
232 0 : const double dfURLat = JDEMGetAngle(psHeader + 43);
233 0 : const double dfURLong = JDEMGetAngle(psHeader + 50);
234 :
235 0 : gt[0] = dfLLLong;
236 0 : gt[3] = dfURLat;
237 0 : gt[1] = (dfURLong - dfLLLong) / GetRasterXSize();
238 0 : gt[2] = 0.0;
239 :
240 0 : gt[4] = 0.0;
241 0 : gt[5] = -1 * (dfURLat - dfLLLat) / GetRasterYSize();
242 :
243 0 : return CE_None;
244 : }
245 :
246 : /************************************************************************/
247 : /* GetSpatialRef() */
248 : /************************************************************************/
249 :
250 0 : const OGRSpatialReference *JDEMDataset::GetSpatialRef() const
251 :
252 : {
253 0 : return &m_oSRS;
254 : }
255 :
256 : /************************************************************************/
257 : /* Identify() */
258 : /************************************************************************/
259 :
260 62895 : int JDEMDataset::Identify(GDALOpenInfo *poOpenInfo)
261 :
262 : {
263 62895 : if (poOpenInfo->nHeaderBytes < HEADER_SIZE)
264 58555 : return FALSE;
265 :
266 : // Confirm that the header has what appears to be dates in the
267 : // expected locations.
268 : // Check if century values seem reasonable.
269 4340 : const char *psHeader = reinterpret_cast<char *>(poOpenInfo->pabyHeader);
270 4340 : if ((!STARTS_WITH_CI(psHeader + 11, "19") &&
271 4336 : !STARTS_WITH_CI(psHeader + 11, "20")) ||
272 9 : (!STARTS_WITH_CI(psHeader + 15, "19") &&
273 5 : !STARTS_WITH_CI(psHeader + 15, "20")) ||
274 4 : (!STARTS_WITH_CI(psHeader + 19, "19") &&
275 0 : !STARTS_WITH_CI(psHeader + 19, "20")))
276 : {
277 4336 : return FALSE;
278 : }
279 :
280 : // Check the extent too. In particular, that we are in the first quadrant,
281 : // as this is only for Japan.
282 4 : const double dfLLLat = JDEMGetAngle(psHeader + 29);
283 4 : const double dfLLLong = JDEMGetAngle(psHeader + 36);
284 4 : const double dfURLat = JDEMGetAngle(psHeader + 43);
285 4 : const double dfURLong = JDEMGetAngle(psHeader + 50);
286 4 : if (dfLLLat > 90 || dfLLLat < 0 || dfLLLong > 180 || dfLLLong < 0 ||
287 4 : dfURLat > 90 || dfURLat < 0 || dfURLong > 180 || dfURLong < 0 ||
288 4 : dfLLLat > dfURLat || dfLLLong > dfURLong)
289 : {
290 0 : return FALSE;
291 : }
292 :
293 4 : return TRUE;
294 : }
295 :
296 : /************************************************************************/
297 : /* Open() */
298 : /************************************************************************/
299 :
300 2 : GDALDataset *JDEMDataset::Open(GDALOpenInfo *poOpenInfo)
301 :
302 : {
303 : // Confirm that the header is compatible with a JDEM dataset.
304 2 : if (!Identify(poOpenInfo))
305 0 : return nullptr;
306 :
307 : // Confirm the requested access is supported.
308 2 : if (poOpenInfo->eAccess == GA_Update)
309 : {
310 0 : ReportUpdateNotSupportedByDriver("JDEM");
311 0 : return nullptr;
312 : }
313 :
314 : // Check that the file pointer from GDALOpenInfo* is available.
315 2 : if (poOpenInfo->fpL == nullptr)
316 : {
317 0 : return nullptr;
318 : }
319 :
320 : // Create a corresponding GDALDataset.
321 4 : auto poDS = std::make_unique<JDEMDataset>();
322 :
323 : // Borrow the file pointer from GDALOpenInfo*.
324 2 : std::swap(poDS->m_fp, poOpenInfo->fpL);
325 :
326 : // Store the header (we have already checked it is at least HEADER_SIZE
327 : // byte large).
328 2 : memcpy(poDS->m_abyHeader, poOpenInfo->pabyHeader, HEADER_SIZE);
329 :
330 2 : const char *psHeader = reinterpret_cast<const char *>(poDS->m_abyHeader);
331 2 : poDS->nRasterXSize = JDEMGetField(psHeader + 23, 3);
332 2 : poDS->nRasterYSize = JDEMGetField(psHeader + 26, 3);
333 2 : if (!GDALCheckDatasetDimensions(poDS->nRasterXSize, poDS->nRasterYSize))
334 : {
335 0 : return nullptr;
336 : }
337 :
338 : // Create band information objects.
339 2 : poDS->SetBand(1, new JDEMRasterBand(poDS.get(), 1));
340 :
341 : // Initialize any PAM information.
342 2 : poDS->SetDescription(poOpenInfo->pszFilename);
343 2 : poDS->TryLoadXML();
344 :
345 : // Check for overviews.
346 2 : poDS->oOvManager.Initialize(poDS.get(), poOpenInfo->pszFilename);
347 :
348 2 : return poDS.release();
349 : }
350 :
351 : /************************************************************************/
352 : /* GDALRegister_JDEM() */
353 : /************************************************************************/
354 :
355 2033 : void GDALRegister_JDEM()
356 :
357 : {
358 2033 : if (GDALGetDriverByName("JDEM") != nullptr)
359 283 : return;
360 :
361 1750 : GDALDriver *poDriver = new GDALDriver();
362 :
363 1750 : poDriver->SetDescription("JDEM");
364 1750 : poDriver->SetMetadataItem(GDAL_DCAP_RASTER, "YES");
365 1750 : poDriver->SetMetadataItem(GDAL_DMD_LONGNAME, "Japanese DEM (.mem)");
366 1750 : poDriver->SetMetadataItem(GDAL_DMD_HELPTOPIC, "drivers/raster/jdem.html");
367 1750 : poDriver->SetMetadataItem(GDAL_DMD_EXTENSION, "mem");
368 1750 : poDriver->SetMetadataItem(GDAL_DCAP_VIRTUALIO, "YES");
369 :
370 1750 : poDriver->pfnOpen = JDEMDataset::Open;
371 1750 : poDriver->pfnIdentify = JDEMDataset::Identify;
372 :
373 1750 : GetGDALDriverManager()->RegisterDriver(poDriver);
374 : }
|
