Line data Source code
1 : /******************************************************************************
2 : *
3 : * Project: GXF Reader
4 : * Purpose: GDAL binding for GXF reader.
5 : * Author: Frank Warmerdam, warmerdam@pobox.com
6 : *
7 : ******************************************************************************
8 : * Copyright (c) 1998, Frank Warmerdam
9 : * Copyright (c) 2008-2012, Even Rouault <even dot rouault at spatialys.com>
10 : *
11 : * SPDX-License-Identifier: MIT
12 : ****************************************************************************/
13 :
14 : #include "gdal_frmts.h"
15 : #include "gdal_pam.h"
16 : #include "gxfopen.h"
17 :
18 : /************************************************************************/
19 : /* ==================================================================== */
20 : /* GXFDataset */
21 : /* ==================================================================== */
22 : /************************************************************************/
23 :
24 : class GXFRasterBand;
25 :
26 : class GXFDataset final : public GDALPamDataset
27 : {
28 : friend class GXFRasterBand;
29 :
30 : GXFHandle hGXF;
31 :
32 : OGRSpatialReference m_oSRS{};
33 : double dfNoDataValue;
34 : GDALDataType eDataType;
35 :
36 : public:
37 : GXFDataset();
38 : ~GXFDataset();
39 :
40 : static GDALDataset *Open(GDALOpenInfo *);
41 :
42 : CPLErr GetGeoTransform(double *padfTransform) override;
43 : const OGRSpatialReference *GetSpatialRef() const override;
44 : };
45 :
46 : /************************************************************************/
47 : /* ==================================================================== */
48 : /* GXFRasterBand */
49 : /* ==================================================================== */
50 : /************************************************************************/
51 :
52 : class GXFRasterBand final : public GDALPamRasterBand
53 : {
54 : friend class GXFDataset;
55 :
56 : public:
57 : GXFRasterBand(GXFDataset *, int);
58 : double GetNoDataValue(int *bGotNoDataValue) override;
59 :
60 : virtual CPLErr IReadBlock(int, int, void *) override;
61 : };
62 :
63 : /************************************************************************/
64 : /* GXFRasterBand() */
65 : /************************************************************************/
66 :
67 4 : GXFRasterBand::GXFRasterBand(GXFDataset *poDSIn, int nBandIn)
68 :
69 : {
70 4 : poDS = poDSIn;
71 4 : nBand = nBandIn;
72 :
73 4 : eDataType = poDSIn->eDataType;
74 :
75 4 : nBlockXSize = poDS->GetRasterXSize();
76 4 : nBlockYSize = 1;
77 4 : }
78 :
79 : /************************************************************************/
80 : /* GetNoDataValue() */
81 : /************************************************************************/
82 :
83 0 : double GXFRasterBand::GetNoDataValue(int *bGotNoDataValue)
84 :
85 : {
86 0 : GXFDataset *poGXF_DS = (GXFDataset *)poDS;
87 0 : if (bGotNoDataValue)
88 0 : *bGotNoDataValue = (fabs(poGXF_DS->dfNoDataValue - -1e12) > .1);
89 0 : if (eDataType == GDT_Float32)
90 0 : return (double)(float)poGXF_DS->dfNoDataValue;
91 :
92 0 : return poGXF_DS->dfNoDataValue;
93 : }
94 :
95 : /************************************************************************/
96 : /* IReadBlock() */
97 : /************************************************************************/
98 :
99 11 : CPLErr GXFRasterBand::IReadBlock(CPL_UNUSED int nBlockXOff, int nBlockYOff,
100 : void *pImage)
101 : {
102 11 : GXFDataset *const poGXF_DS = (GXFDataset *)poDS;
103 :
104 11 : if (eDataType == GDT_Float32)
105 : {
106 11 : double *padfBuffer = (double *)VSIMalloc2(sizeof(double), nBlockXSize);
107 11 : if (padfBuffer == nullptr)
108 0 : return CE_Failure;
109 : const CPLErr eErr =
110 11 : GXFGetScanline(poGXF_DS->hGXF, nBlockYOff, padfBuffer);
111 :
112 11 : float *pafBuffer = (float *)pImage;
113 103 : for (int i = 0; i < nBlockXSize; i++)
114 92 : pafBuffer[i] = (float)padfBuffer[i];
115 :
116 11 : CPLFree(padfBuffer);
117 :
118 11 : return eErr;
119 : }
120 :
121 : const CPLErr eErr =
122 0 : eDataType == GDT_Float64
123 0 : ? GXFGetScanline(poGXF_DS->hGXF, nBlockYOff, (double *)pImage)
124 0 : : CE_Failure;
125 :
126 0 : return eErr;
127 : }
128 :
129 : /************************************************************************/
130 : /* ==================================================================== */
131 : /* GXFDataset */
132 : /* ==================================================================== */
133 : /************************************************************************/
134 :
135 : /************************************************************************/
136 : /* GXFDataset() */
137 : /************************************************************************/
138 :
139 4 : GXFDataset::GXFDataset()
140 4 : : hGXF(nullptr), dfNoDataValue(0), eDataType(GDT_Float32)
141 : {
142 4 : m_oSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
143 4 : }
144 :
145 : /************************************************************************/
146 : /* ~GXFDataset() */
147 : /************************************************************************/
148 :
149 8 : GXFDataset::~GXFDataset()
150 :
151 : {
152 4 : FlushCache(true);
153 4 : if (hGXF != nullptr)
154 4 : GXFClose(hGXF);
155 8 : }
156 :
157 : /************************************************************************/
158 : /* GetGeoTransform() */
159 : /************************************************************************/
160 :
161 0 : CPLErr GXFDataset::GetGeoTransform(double *padfTransform)
162 :
163 : {
164 0 : double dfXOrigin = 0.0;
165 0 : double dfYOrigin = 0.0;
166 0 : double dfXSize = 0.0;
167 0 : double dfYSize = 0.0;
168 0 : double dfRotation = 0.0;
169 :
170 0 : const CPLErr eErr = GXFGetPosition(hGXF, &dfXOrigin, &dfYOrigin, &dfXSize,
171 : &dfYSize, &dfRotation);
172 :
173 0 : if (eErr != CE_None)
174 0 : return eErr;
175 :
176 : // Transform to radians.
177 0 : dfRotation = (dfRotation / 360.0) * 2.0 * M_PI;
178 :
179 0 : padfTransform[1] = dfXSize * cos(dfRotation);
180 0 : padfTransform[2] = dfYSize * sin(dfRotation);
181 0 : padfTransform[4] = dfXSize * sin(dfRotation);
182 0 : padfTransform[5] = -1 * dfYSize * cos(dfRotation);
183 :
184 : // take into account that GXF is point or center of pixel oriented.
185 0 : padfTransform[0] =
186 0 : dfXOrigin - 0.5 * padfTransform[1] - 0.5 * padfTransform[2];
187 0 : padfTransform[3] =
188 0 : dfYOrigin - 0.5 * padfTransform[4] - 0.5 * padfTransform[5];
189 :
190 0 : return CE_None;
191 : }
192 :
193 : /************************************************************************/
194 : /* GetSpatialRef() */
195 : /************************************************************************/
196 :
197 1 : const OGRSpatialReference *GXFDataset::GetSpatialRef() const
198 : {
199 1 : return m_oSRS.IsEmpty() ? nullptr : &m_oSRS;
200 : }
201 :
202 : /************************************************************************/
203 : /* Open() */
204 : /************************************************************************/
205 :
206 30245 : GDALDataset *GXFDataset::Open(GDALOpenInfo *poOpenInfo)
207 :
208 : {
209 : /* -------------------------------------------------------------------- */
210 : /* Before trying GXFOpen() we first verify that there is at */
211 : /* least one "\n#keyword" type signature in the first chunk of */
212 : /* the file. */
213 : /* -------------------------------------------------------------------- */
214 30245 : if (poOpenInfo->nHeaderBytes < 50 || poOpenInfo->fpL == nullptr)
215 27415 : return nullptr;
216 :
217 2830 : bool bFoundKeyword = false;
218 2830 : bool bFoundIllegal = false;
219 710605 : for (int i = 0; i < poOpenInfo->nHeaderBytes - 1; i++)
220 : {
221 709805 : if ((poOpenInfo->pabyHeader[i] == 10 ||
222 692498 : poOpenInfo->pabyHeader[i] == 13) &&
223 18431 : poOpenInfo->pabyHeader[i + 1] == '#')
224 : {
225 159 : if (STARTS_WITH((const char *)poOpenInfo->pabyHeader + i + 2,
226 : "include"))
227 0 : return nullptr;
228 159 : if (STARTS_WITH((const char *)poOpenInfo->pabyHeader + i + 2,
229 : "define"))
230 0 : return nullptr;
231 159 : if (STARTS_WITH((const char *)poOpenInfo->pabyHeader + i + 2,
232 : "ifdef"))
233 0 : return nullptr;
234 159 : bFoundKeyword = true;
235 : }
236 709805 : if (poOpenInfo->pabyHeader[i] == 0)
237 : {
238 2030 : bFoundIllegal = true;
239 2030 : break;
240 : }
241 : }
242 :
243 2830 : if (!bFoundKeyword || bFoundIllegal)
244 2790 : return nullptr;
245 :
246 : /* -------------------------------------------------------------------- */
247 : /* At this point it is plausible that this is a GXF file, but */
248 : /* we also now verify that there is a #GRID keyword before */
249 : /* passing it off to GXFOpen(). We check in the first 50K. */
250 : /* -------------------------------------------------------------------- */
251 40 : CPL_IGNORE_RET_VAL(poOpenInfo->TryToIngest(50000));
252 40 : bool bGotGrid = false;
253 :
254 40 : const char *pszBigBuf = (const char *)poOpenInfo->pabyHeader;
255 9209 : for (int i = 0; i < poOpenInfo->nHeaderBytes - 5 && !bGotGrid; i++)
256 : {
257 9169 : if (pszBigBuf[i] == '#' && STARTS_WITH_CI(pszBigBuf + i + 1, "GRID"))
258 4 : bGotGrid = true;
259 : }
260 :
261 40 : if (!bGotGrid)
262 36 : return nullptr;
263 :
264 4 : VSIFCloseL(poOpenInfo->fpL);
265 4 : poOpenInfo->fpL = nullptr;
266 :
267 : /* -------------------------------------------------------------------- */
268 : /* Try opening the dataset. */
269 : /* -------------------------------------------------------------------- */
270 :
271 4 : GXFHandle l_hGXF = GXFOpen(poOpenInfo->pszFilename);
272 :
273 4 : if (l_hGXF == nullptr)
274 0 : return nullptr;
275 :
276 : /* -------------------------------------------------------------------- */
277 : /* Confirm the requested access is supported. */
278 : /* -------------------------------------------------------------------- */
279 4 : if (poOpenInfo->eAccess == GA_Update)
280 : {
281 0 : GXFClose(l_hGXF);
282 0 : CPLError(CE_Failure, CPLE_NotSupported,
283 : "The GXF driver does not support update access to existing"
284 : " datasets.");
285 0 : return nullptr;
286 : }
287 :
288 : /* -------------------------------------------------------------------- */
289 : /* Create a corresponding GDALDataset. */
290 : /* -------------------------------------------------------------------- */
291 4 : GXFDataset *poDS = new GXFDataset();
292 :
293 4 : const char *pszGXFDataType = CPLGetConfigOption("GXF_DATATYPE", "Float32");
294 4 : GDALDataType eDT = GDALGetDataTypeByName(pszGXFDataType);
295 4 : if (!(eDT == GDT_Float32 || eDT == GDT_Float64))
296 : {
297 0 : CPLError(CE_Warning, CPLE_NotSupported,
298 : "Unsupported value for GXF_DATATYPE : %s", pszGXFDataType);
299 0 : eDT = GDT_Float32;
300 : }
301 :
302 4 : poDS->hGXF = l_hGXF;
303 4 : poDS->eDataType = eDT;
304 :
305 : /* -------------------------------------------------------------------- */
306 : /* Establish the projection. */
307 : /* -------------------------------------------------------------------- */
308 4 : char *pszProjection = GXFGetMapProjectionAsOGCWKT(l_hGXF);
309 4 : if (pszProjection && pszProjection[0] != '\0')
310 2 : poDS->m_oSRS.importFromWkt(pszProjection);
311 4 : CPLFree(pszProjection);
312 :
313 : /* -------------------------------------------------------------------- */
314 : /* Capture some information from the file that is of interest. */
315 : /* -------------------------------------------------------------------- */
316 4 : GXFGetRawInfo(l_hGXF, &(poDS->nRasterXSize), &(poDS->nRasterYSize), nullptr,
317 : nullptr, nullptr, &(poDS->dfNoDataValue));
318 :
319 4 : if (poDS->nRasterXSize <= 0 || poDS->nRasterYSize <= 0)
320 : {
321 0 : CPLError(CE_Failure, CPLE_AppDefined, "Invalid dimensions : %d x %d",
322 : poDS->nRasterXSize, poDS->nRasterYSize);
323 0 : delete poDS;
324 0 : return nullptr;
325 : }
326 :
327 : /* -------------------------------------------------------------------- */
328 : /* Create band information objects. */
329 : /* -------------------------------------------------------------------- */
330 4 : poDS->nBands = 1;
331 4 : poDS->SetBand(1, new GXFRasterBand(poDS, 1));
332 :
333 : /* -------------------------------------------------------------------- */
334 : /* Initialize any PAM information. */
335 : /* -------------------------------------------------------------------- */
336 4 : poDS->SetDescription(poOpenInfo->pszFilename);
337 4 : poDS->TryLoadXML();
338 :
339 : /* -------------------------------------------------------------------- */
340 : /* Check for external overviews. */
341 : /* -------------------------------------------------------------------- */
342 8 : poDS->oOvManager.Initialize(poDS, poOpenInfo->pszFilename,
343 4 : poOpenInfo->GetSiblingFiles());
344 :
345 4 : return poDS;
346 : }
347 :
348 : /************************************************************************/
349 : /* GDALRegister_GXF() */
350 : /************************************************************************/
351 :
352 1595 : void GDALRegister_GXF()
353 :
354 : {
355 1595 : if (GDALGetDriverByName("GXF") != nullptr)
356 302 : return;
357 :
358 1293 : GDALDriver *poDriver = new GDALDriver();
359 :
360 1293 : poDriver->SetDescription("GXF");
361 1293 : poDriver->SetMetadataItem(GDAL_DCAP_RASTER, "YES");
362 1293 : poDriver->SetMetadataItem(GDAL_DMD_LONGNAME,
363 1293 : "GeoSoft Grid Exchange Format");
364 1293 : poDriver->SetMetadataItem(GDAL_DMD_HELPTOPIC, "drivers/raster/gxf.html");
365 1293 : poDriver->SetMetadataItem(GDAL_DMD_EXTENSION, "gxf");
366 1293 : poDriver->SetMetadataItem(GDAL_DCAP_VIRTUALIO, "YES");
367 :
368 1293 : poDriver->pfnOpen = GXFDataset::Open;
369 :
370 1293 : GetGDALDriverManager()->RegisterDriver(poDriver);
371 : }
|
