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 : * Permission is hereby granted, free of charge, to any person obtaining a
12 : * copy of this software and associated documentation files (the "Software"),
13 : * to deal in the Software without restriction, including without limitation
14 : * the rights to use, copy, modify, merge, publish, distribute, sublicense,
15 : * and/or sell copies of the Software, and to permit persons to whom the
16 : * Software is furnished to do so, subject to the following conditions:
17 : *
18 : * The above copyright notice and this permission notice shall be included
19 : * in all copies or substantial portions of the Software.
20 : *
21 : * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
22 : * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
23 : * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
24 : * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
25 : * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
26 : * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
27 : * DEALINGS IN THE SOFTWARE.
28 : ****************************************************************************/
29 :
30 : #include "cpl_string.h"
31 : #include "gdal_frmts.h"
32 : #include "ogr_srs_api.h"
33 : #include "rawdataset.h"
34 :
35 : #include <algorithm>
36 :
37 : /**
38 :
39 : NOAA .LOS/.LAS Datum Grid Shift Format
40 :
41 : Also used for .geo file from https://geodesy.noaa.gov/GEOID/MEXICO97/
42 :
43 : All values are little endian
44 :
45 : Header
46 : ------
47 :
48 : char[56] "NADCON EXTRACTED REGION" or
49 : "GEOID EXTRACTED REGION "
50 : char[8] "NADGRD " or
51 : "GEOGRD "
52 : int32 grid width
53 : int32 grid height
54 : int32 z count (1)
55 : float32 origin longitude
56 : float32 grid cell width longitude
57 : float32 origin latitude
58 : float32 grid cell height latitude
59 : float32 angle (0.0)
60 :
61 : Data
62 : ----
63 :
64 : int32 ? always 0
65 : float32*gridwidth offset in arcseconds (or in metres for geoids)
66 :
67 : Note that the record length is always gridwidth*4 + 4, and
68 : even the header record is this length though it means some waste.
69 :
70 : **/
71 :
72 : /************************************************************************/
73 : /* ==================================================================== */
74 : /* LOSLASDataset */
75 : /* ==================================================================== */
76 : /************************************************************************/
77 :
78 : class LOSLASDataset final : public RawDataset
79 : {
80 : VSILFILE *fpImage; // image data file.
81 :
82 : int nRecordLength;
83 :
84 : OGRSpatialReference m_oSRS{};
85 : double adfGeoTransform[6];
86 :
87 : CPL_DISALLOW_COPY_ASSIGN(LOSLASDataset)
88 :
89 : CPLErr Close() override;
90 :
91 : public:
92 : LOSLASDataset();
93 : ~LOSLASDataset() override;
94 :
95 : CPLErr GetGeoTransform(double *padfTransform) override;
96 :
97 1 : const OGRSpatialReference *GetSpatialRef() const override
98 : {
99 1 : return &m_oSRS;
100 : }
101 :
102 : static GDALDataset *Open(GDALOpenInfo *);
103 : static int Identify(GDALOpenInfo *);
104 : };
105 :
106 : /************************************************************************/
107 : /* ==================================================================== */
108 : /* LOSLASDataset */
109 : /* ==================================================================== */
110 : /************************************************************************/
111 :
112 : /************************************************************************/
113 : /* LOSLASDataset() */
114 : /************************************************************************/
115 :
116 2 : LOSLASDataset::LOSLASDataset() : fpImage(nullptr), nRecordLength(0)
117 : {
118 2 : m_oSRS.SetFromUserInput(SRS_WKT_WGS84_LAT_LONG);
119 2 : m_oSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
120 :
121 2 : memset(adfGeoTransform, 0, sizeof(adfGeoTransform));
122 2 : }
123 :
124 : /************************************************************************/
125 : /* ~LOSLASDataset() */
126 : /************************************************************************/
127 :
128 4 : LOSLASDataset::~LOSLASDataset()
129 :
130 : {
131 2 : LOSLASDataset::Close();
132 4 : }
133 :
134 : /************************************************************************/
135 : /* Close() */
136 : /************************************************************************/
137 :
138 4 : CPLErr LOSLASDataset::Close()
139 : {
140 4 : CPLErr eErr = CE_None;
141 4 : if (nOpenFlags != OPEN_FLAGS_CLOSED)
142 : {
143 2 : if (LOSLASDataset::FlushCache(true) != CE_None)
144 0 : eErr = CE_Failure;
145 :
146 2 : if (fpImage)
147 : {
148 2 : if (VSIFCloseL(fpImage) != 0)
149 : {
150 0 : CPLError(CE_Failure, CPLE_FileIO, "I/O error");
151 0 : eErr = CE_Failure;
152 : }
153 : }
154 :
155 2 : if (GDALPamDataset::Close() != CE_None)
156 0 : eErr = CE_Failure;
157 : }
158 4 : return eErr;
159 : }
160 :
161 : /************************************************************************/
162 : /* Identify() */
163 : /************************************************************************/
164 :
165 48756 : int LOSLASDataset::Identify(GDALOpenInfo *poOpenInfo)
166 :
167 : {
168 48756 : if (poOpenInfo->nHeaderBytes < 64)
169 45524 : return FALSE;
170 :
171 : #ifndef FUZZING_BUILD_MODE_UNSAFE_FOR_PRODUCTION
172 3232 : const char *pszExt = CPLGetExtension(poOpenInfo->pszFilename);
173 3232 : if (!EQUAL(pszExt, "las") && !EQUAL(pszExt, "los") && !EQUAL(pszExt, "geo"))
174 3228 : return FALSE;
175 : #endif
176 :
177 4 : if (!STARTS_WITH_CI((const char *)poOpenInfo->pabyHeader + 56, "NADGRD") &&
178 0 : !STARTS_WITH_CI((const char *)poOpenInfo->pabyHeader + 56, "GEOGRD"))
179 0 : return FALSE;
180 :
181 4 : return TRUE;
182 : }
183 :
184 : /************************************************************************/
185 : /* Open() */
186 : /************************************************************************/
187 :
188 2 : GDALDataset *LOSLASDataset::Open(GDALOpenInfo *poOpenInfo)
189 :
190 : {
191 2 : if (!Identify(poOpenInfo) || poOpenInfo->fpL == nullptr)
192 0 : return nullptr;
193 :
194 : /* -------------------------------------------------------------------- */
195 : /* Confirm the requested access is supported. */
196 : /* -------------------------------------------------------------------- */
197 2 : if (poOpenInfo->eAccess == GA_Update)
198 : {
199 0 : CPLError(CE_Failure, CPLE_NotSupported,
200 : "The LOSLAS driver does not support update access to existing"
201 : " datasets.");
202 0 : return nullptr;
203 : }
204 :
205 : /* -------------------------------------------------------------------- */
206 : /* Create a corresponding GDALDataset. */
207 : /* -------------------------------------------------------------------- */
208 4 : auto poDS = std::make_unique<LOSLASDataset>();
209 2 : std::swap(poDS->fpImage, poOpenInfo->fpL);
210 :
211 : /* -------------------------------------------------------------------- */
212 : /* Read the header. */
213 : /* -------------------------------------------------------------------- */
214 2 : CPL_IGNORE_RET_VAL(VSIFSeekL(poDS->fpImage, 64, SEEK_SET));
215 :
216 2 : CPL_IGNORE_RET_VAL(VSIFReadL(&(poDS->nRasterXSize), 4, 1, poDS->fpImage));
217 2 : CPL_IGNORE_RET_VAL(VSIFReadL(&(poDS->nRasterYSize), 4, 1, poDS->fpImage));
218 :
219 2 : CPL_LSBPTR32(&(poDS->nRasterXSize));
220 2 : CPL_LSBPTR32(&(poDS->nRasterYSize));
221 :
222 4 : if (!GDALCheckDatasetDimensions(poDS->nRasterXSize, poDS->nRasterYSize) ||
223 2 : poDS->nRasterXSize > (INT_MAX - 4) / 4)
224 : {
225 0 : return nullptr;
226 : }
227 :
228 2 : CPL_IGNORE_RET_VAL(VSIFSeekL(poDS->fpImage, 76, SEEK_SET));
229 :
230 : float min_lon, min_lat, delta_lon, delta_lat;
231 :
232 2 : CPL_IGNORE_RET_VAL(VSIFReadL(&min_lon, 4, 1, poDS->fpImage));
233 2 : CPL_IGNORE_RET_VAL(VSIFReadL(&delta_lon, 4, 1, poDS->fpImage));
234 2 : CPL_IGNORE_RET_VAL(VSIFReadL(&min_lat, 4, 1, poDS->fpImage));
235 2 : CPL_IGNORE_RET_VAL(VSIFReadL(&delta_lat, 4, 1, poDS->fpImage));
236 :
237 2 : CPL_LSBPTR32(&min_lon);
238 2 : CPL_LSBPTR32(&delta_lon);
239 2 : CPL_LSBPTR32(&min_lat);
240 2 : CPL_LSBPTR32(&delta_lat);
241 :
242 2 : poDS->nRecordLength = poDS->nRasterXSize * 4 + 4;
243 :
244 : /* -------------------------------------------------------------------- */
245 : /* Create band information object. */
246 : /* */
247 : /* Note we are setting up to read from the last image record to */
248 : /* the first since the data comes with the southern most record */
249 : /* first, not the northernmost like we would want. */
250 : /* -------------------------------------------------------------------- */
251 : auto poBand = RawRasterBand::Create(
252 4 : poDS.get(), 1, poDS->fpImage,
253 2 : static_cast<vsi_l_offset>(poDS->nRasterYSize) * poDS->nRecordLength + 4,
254 2 : 4, -1 * poDS->nRecordLength, GDT_Float32,
255 : RawRasterBand::ByteOrder::ORDER_LITTLE_ENDIAN,
256 6 : RawRasterBand::OwnFP::NO);
257 2 : if (!poBand)
258 0 : return nullptr;
259 2 : poDS->SetBand(1, std::move(poBand));
260 :
261 2 : if (EQUAL(CPLGetExtension(poOpenInfo->pszFilename), "las"))
262 : {
263 0 : poDS->GetRasterBand(1)->SetDescription("Latitude Offset (arc seconds)");
264 : }
265 2 : else if (EQUAL(CPLGetExtension(poOpenInfo->pszFilename), "los"))
266 : {
267 2 : poDS->GetRasterBand(1)->SetDescription(
268 2 : "Longitude Offset (arc seconds)");
269 2 : poDS->GetRasterBand(1)->SetMetadataItem("positive_value", "west");
270 : }
271 0 : else if (EQUAL(CPLGetExtension(poOpenInfo->pszFilename), "geo"))
272 : {
273 0 : poDS->GetRasterBand(1)->SetDescription("Geoid undulation (meters)");
274 : }
275 :
276 : /* -------------------------------------------------------------------- */
277 : /* Setup georeferencing. */
278 : /* -------------------------------------------------------------------- */
279 2 : poDS->adfGeoTransform[0] = min_lon - delta_lon * 0.5;
280 2 : poDS->adfGeoTransform[1] = delta_lon;
281 2 : poDS->adfGeoTransform[2] = 0.0;
282 2 : poDS->adfGeoTransform[3] = min_lat + (poDS->nRasterYSize - 0.5) * delta_lat;
283 2 : poDS->adfGeoTransform[4] = 0.0;
284 2 : poDS->adfGeoTransform[5] = -1.0 * delta_lat;
285 :
286 : /* -------------------------------------------------------------------- */
287 : /* Initialize any PAM information. */
288 : /* -------------------------------------------------------------------- */
289 2 : poDS->SetDescription(poOpenInfo->pszFilename);
290 2 : poDS->TryLoadXML();
291 :
292 : /* -------------------------------------------------------------------- */
293 : /* Check for overviews. */
294 : /* -------------------------------------------------------------------- */
295 2 : poDS->oOvManager.Initialize(poDS.get(), poOpenInfo->pszFilename);
296 :
297 2 : return poDS.release();
298 : }
299 :
300 : /************************************************************************/
301 : /* GetGeoTransform() */
302 : /************************************************************************/
303 :
304 1 : CPLErr LOSLASDataset::GetGeoTransform(double *padfTransform)
305 :
306 : {
307 1 : memcpy(padfTransform, adfGeoTransform, sizeof(double) * 6);
308 1 : return CE_None;
309 : }
310 :
311 : /************************************************************************/
312 : /* GDALRegister_LOSLAS() */
313 : /************************************************************************/
314 :
315 1511 : void GDALRegister_LOSLAS()
316 :
317 : {
318 1511 : if (GDALGetDriverByName("LOSLAS") != nullptr)
319 295 : return;
320 :
321 1216 : GDALDriver *poDriver = new GDALDriver();
322 :
323 1216 : poDriver->SetDescription("LOSLAS");
324 1216 : poDriver->SetMetadataItem(GDAL_DCAP_RASTER, "YES");
325 1216 : poDriver->SetMetadataItem(GDAL_DMD_LONGNAME,
326 1216 : "NADCON .los/.las Datum Grid Shift");
327 1216 : poDriver->SetMetadataItem(GDAL_DCAP_VIRTUALIO, "YES");
328 :
329 1216 : poDriver->pfnOpen = LOSLASDataset::Open;
330 1216 : poDriver->pfnIdentify = LOSLASDataset::Identify;
331 :
332 1216 : GetGDALDriverManager()->RegisterDriver(poDriver);
333 : }
|
