Line data Source code
1 : /******************************************************************************
2 : *
3 : * Project: Horizontal Datum Formats
4 : * Purpose: Implementation of NOAA/NADCON los/las datum shift format.
5 : * Author: Frank Warmerdam, warmerdam@pobox.com
6 : * Financial Support: i-cubed (http://www.i-cubed.com)
7 : *
8 : ******************************************************************************
9 : * Copyright (c) 2010, Frank Warmerdam
10 : *
11 : * SPDX-License-Identifier: MIT
12 : ****************************************************************************/
13 :
14 : #include "cpl_string.h"
15 : #include "gdal_frmts.h"
16 : #include "ogr_srs_api.h"
17 : #include "rawdataset.h"
18 :
19 : #include <algorithm>
20 :
21 : /**
22 :
23 : NOAA .LOS/.LAS Datum Grid Shift Format
24 :
25 : Also used for .geo file from https://geodesy.noaa.gov/GEOID/MEXICO97/
26 :
27 : All values are little endian
28 :
29 : Header
30 : ------
31 :
32 : char[56] "NADCON EXTRACTED REGION" or
33 : "GEOID EXTRACTED REGION "
34 : char[8] "NADGRD " or
35 : "GEOGRD "
36 : int32 grid width
37 : int32 grid height
38 : int32 z count (1)
39 : float32 origin longitude
40 : float32 grid cell width longitude
41 : float32 origin latitude
42 : float32 grid cell height latitude
43 : float32 angle (0.0)
44 :
45 : Data
46 : ----
47 :
48 : int32 ? always 0
49 : float32*gridwidth offset in arcseconds (or in metres for geoids)
50 :
51 : Note that the record length is always gridwidth*4 + 4, and
52 : even the header record is this length though it means some waste.
53 :
54 : **/
55 :
56 : /************************************************************************/
57 : /* ==================================================================== */
58 : /* LOSLASDataset */
59 : /* ==================================================================== */
60 : /************************************************************************/
61 :
62 : class LOSLASDataset final : public RawDataset
63 : {
64 : VSILFILE *fpImage; // image data file.
65 :
66 : int nRecordLength;
67 :
68 : OGRSpatialReference m_oSRS{};
69 : double adfGeoTransform[6];
70 :
71 : CPL_DISALLOW_COPY_ASSIGN(LOSLASDataset)
72 :
73 : CPLErr Close() override;
74 :
75 : public:
76 : LOSLASDataset();
77 : ~LOSLASDataset() override;
78 :
79 : CPLErr GetGeoTransform(double *padfTransform) override;
80 :
81 1 : const OGRSpatialReference *GetSpatialRef() const override
82 : {
83 1 : return &m_oSRS;
84 : }
85 :
86 : static GDALDataset *Open(GDALOpenInfo *);
87 : static int Identify(GDALOpenInfo *);
88 : };
89 :
90 : /************************************************************************/
91 : /* ==================================================================== */
92 : /* LOSLASDataset */
93 : /* ==================================================================== */
94 : /************************************************************************/
95 :
96 : /************************************************************************/
97 : /* LOSLASDataset() */
98 : /************************************************************************/
99 :
100 2 : LOSLASDataset::LOSLASDataset() : fpImage(nullptr), nRecordLength(0)
101 : {
102 2 : m_oSRS.SetFromUserInput(SRS_WKT_WGS84_LAT_LONG);
103 2 : m_oSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
104 :
105 2 : memset(adfGeoTransform, 0, sizeof(adfGeoTransform));
106 2 : }
107 :
108 : /************************************************************************/
109 : /* ~LOSLASDataset() */
110 : /************************************************************************/
111 :
112 4 : LOSLASDataset::~LOSLASDataset()
113 :
114 : {
115 2 : LOSLASDataset::Close();
116 4 : }
117 :
118 : /************************************************************************/
119 : /* Close() */
120 : /************************************************************************/
121 :
122 4 : CPLErr LOSLASDataset::Close()
123 : {
124 4 : CPLErr eErr = CE_None;
125 4 : if (nOpenFlags != OPEN_FLAGS_CLOSED)
126 : {
127 2 : if (LOSLASDataset::FlushCache(true) != CE_None)
128 0 : eErr = CE_Failure;
129 :
130 2 : if (fpImage)
131 : {
132 2 : if (VSIFCloseL(fpImage) != 0)
133 : {
134 0 : CPLError(CE_Failure, CPLE_FileIO, "I/O error");
135 0 : eErr = CE_Failure;
136 : }
137 : }
138 :
139 2 : if (GDALPamDataset::Close() != CE_None)
140 0 : eErr = CE_Failure;
141 : }
142 4 : return eErr;
143 : }
144 :
145 : /************************************************************************/
146 : /* Identify() */
147 : /************************************************************************/
148 :
149 51446 : int LOSLASDataset::Identify(GDALOpenInfo *poOpenInfo)
150 :
151 : {
152 51446 : if (poOpenInfo->nHeaderBytes < 64)
153 48054 : return FALSE;
154 :
155 : #ifndef FUZZING_BUILD_MODE_UNSAFE_FOR_PRODUCTION
156 3392 : const char *pszExt = CPLGetExtension(poOpenInfo->pszFilename);
157 3392 : if (!EQUAL(pszExt, "las") && !EQUAL(pszExt, "los") && !EQUAL(pszExt, "geo"))
158 3388 : return FALSE;
159 : #endif
160 :
161 4 : if (!STARTS_WITH_CI((const char *)poOpenInfo->pabyHeader + 56, "NADGRD") &&
162 0 : !STARTS_WITH_CI((const char *)poOpenInfo->pabyHeader + 56, "GEOGRD"))
163 0 : return FALSE;
164 :
165 4 : return TRUE;
166 : }
167 :
168 : /************************************************************************/
169 : /* Open() */
170 : /************************************************************************/
171 :
172 2 : GDALDataset *LOSLASDataset::Open(GDALOpenInfo *poOpenInfo)
173 :
174 : {
175 2 : if (!Identify(poOpenInfo) || poOpenInfo->fpL == nullptr)
176 0 : return nullptr;
177 :
178 : /* -------------------------------------------------------------------- */
179 : /* Confirm the requested access is supported. */
180 : /* -------------------------------------------------------------------- */
181 2 : if (poOpenInfo->eAccess == GA_Update)
182 : {
183 0 : CPLError(CE_Failure, CPLE_NotSupported,
184 : "The LOSLAS driver does not support update access to existing"
185 : " datasets.");
186 0 : return nullptr;
187 : }
188 :
189 : /* -------------------------------------------------------------------- */
190 : /* Create a corresponding GDALDataset. */
191 : /* -------------------------------------------------------------------- */
192 4 : auto poDS = std::make_unique<LOSLASDataset>();
193 2 : std::swap(poDS->fpImage, poOpenInfo->fpL);
194 :
195 : /* -------------------------------------------------------------------- */
196 : /* Read the header. */
197 : /* -------------------------------------------------------------------- */
198 2 : CPL_IGNORE_RET_VAL(VSIFSeekL(poDS->fpImage, 64, SEEK_SET));
199 :
200 2 : CPL_IGNORE_RET_VAL(VSIFReadL(&(poDS->nRasterXSize), 4, 1, poDS->fpImage));
201 2 : CPL_IGNORE_RET_VAL(VSIFReadL(&(poDS->nRasterYSize), 4, 1, poDS->fpImage));
202 :
203 2 : CPL_LSBPTR32(&(poDS->nRasterXSize));
204 2 : CPL_LSBPTR32(&(poDS->nRasterYSize));
205 :
206 4 : if (!GDALCheckDatasetDimensions(poDS->nRasterXSize, poDS->nRasterYSize) ||
207 2 : poDS->nRasterXSize > (INT_MAX - 4) / 4)
208 : {
209 0 : return nullptr;
210 : }
211 :
212 2 : CPL_IGNORE_RET_VAL(VSIFSeekL(poDS->fpImage, 76, SEEK_SET));
213 :
214 : float min_lon, min_lat, delta_lon, delta_lat;
215 :
216 2 : CPL_IGNORE_RET_VAL(VSIFReadL(&min_lon, 4, 1, poDS->fpImage));
217 2 : CPL_IGNORE_RET_VAL(VSIFReadL(&delta_lon, 4, 1, poDS->fpImage));
218 2 : CPL_IGNORE_RET_VAL(VSIFReadL(&min_lat, 4, 1, poDS->fpImage));
219 2 : CPL_IGNORE_RET_VAL(VSIFReadL(&delta_lat, 4, 1, poDS->fpImage));
220 :
221 2 : CPL_LSBPTR32(&min_lon);
222 2 : CPL_LSBPTR32(&delta_lon);
223 2 : CPL_LSBPTR32(&min_lat);
224 2 : CPL_LSBPTR32(&delta_lat);
225 :
226 2 : poDS->nRecordLength = poDS->nRasterXSize * 4 + 4;
227 :
228 : /* -------------------------------------------------------------------- */
229 : /* Create band information object. */
230 : /* */
231 : /* Note we are setting up to read from the last image record to */
232 : /* the first since the data comes with the southern most record */
233 : /* first, not the northernmost like we would want. */
234 : /* -------------------------------------------------------------------- */
235 : auto poBand = RawRasterBand::Create(
236 4 : poDS.get(), 1, poDS->fpImage,
237 2 : static_cast<vsi_l_offset>(poDS->nRasterYSize) * poDS->nRecordLength + 4,
238 2 : 4, -1 * poDS->nRecordLength, GDT_Float32,
239 : RawRasterBand::ByteOrder::ORDER_LITTLE_ENDIAN,
240 6 : RawRasterBand::OwnFP::NO);
241 2 : if (!poBand)
242 0 : return nullptr;
243 2 : poDS->SetBand(1, std::move(poBand));
244 :
245 2 : if (EQUAL(CPLGetExtension(poOpenInfo->pszFilename), "las"))
246 : {
247 0 : poDS->GetRasterBand(1)->SetDescription("Latitude Offset (arc seconds)");
248 : }
249 2 : else if (EQUAL(CPLGetExtension(poOpenInfo->pszFilename), "los"))
250 : {
251 2 : poDS->GetRasterBand(1)->SetDescription(
252 2 : "Longitude Offset (arc seconds)");
253 2 : poDS->GetRasterBand(1)->SetMetadataItem("positive_value", "west");
254 : }
255 0 : else if (EQUAL(CPLGetExtension(poOpenInfo->pszFilename), "geo"))
256 : {
257 0 : poDS->GetRasterBand(1)->SetDescription("Geoid undulation (meters)");
258 : }
259 :
260 : /* -------------------------------------------------------------------- */
261 : /* Setup georeferencing. */
262 : /* -------------------------------------------------------------------- */
263 2 : poDS->adfGeoTransform[0] = min_lon - delta_lon * 0.5;
264 2 : poDS->adfGeoTransform[1] = delta_lon;
265 2 : poDS->adfGeoTransform[2] = 0.0;
266 2 : poDS->adfGeoTransform[3] = min_lat + (poDS->nRasterYSize - 0.5) * delta_lat;
267 2 : poDS->adfGeoTransform[4] = 0.0;
268 2 : poDS->adfGeoTransform[5] = -1.0 * delta_lat;
269 :
270 : /* -------------------------------------------------------------------- */
271 : /* Initialize any PAM information. */
272 : /* -------------------------------------------------------------------- */
273 2 : poDS->SetDescription(poOpenInfo->pszFilename);
274 2 : poDS->TryLoadXML();
275 :
276 : /* -------------------------------------------------------------------- */
277 : /* Check for overviews. */
278 : /* -------------------------------------------------------------------- */
279 2 : poDS->oOvManager.Initialize(poDS.get(), poOpenInfo->pszFilename);
280 :
281 2 : return poDS.release();
282 : }
283 :
284 : /************************************************************************/
285 : /* GetGeoTransform() */
286 : /************************************************************************/
287 :
288 1 : CPLErr LOSLASDataset::GetGeoTransform(double *padfTransform)
289 :
290 : {
291 1 : memcpy(padfTransform, adfGeoTransform, sizeof(double) * 6);
292 1 : return CE_None;
293 : }
294 :
295 : /************************************************************************/
296 : /* GDALRegister_LOSLAS() */
297 : /************************************************************************/
298 :
299 1595 : void GDALRegister_LOSLAS()
300 :
301 : {
302 1595 : if (GDALGetDriverByName("LOSLAS") != nullptr)
303 302 : return;
304 :
305 1293 : GDALDriver *poDriver = new GDALDriver();
306 :
307 1293 : poDriver->SetDescription("LOSLAS");
308 1293 : poDriver->SetMetadataItem(GDAL_DCAP_RASTER, "YES");
309 1293 : poDriver->SetMetadataItem(GDAL_DMD_LONGNAME,
310 1293 : "NADCON .los/.las Datum Grid Shift");
311 1293 : poDriver->SetMetadataItem(GDAL_DCAP_VIRTUALIO, "YES");
312 :
313 1293 : poDriver->pfnOpen = LOSLASDataset::Open;
314 1293 : poDriver->pfnIdentify = LOSLASDataset::Identify;
315 :
316 1293 : GetGDALDriverManager()->RegisterDriver(poDriver);
317 : }
|
