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