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