Line data Source code
1 : /******************************************************************************
2 : *
3 : * Project: CTG driver
4 : * Purpose: GDALDataset driver for CTG dataset.
5 : * Author: Even Rouault, <even dot rouault at spatialys.com>
6 : *
7 : ******************************************************************************
8 : * Copyright (c) 2011, Even Rouault <even dot rouault at spatialys.com>
9 : *
10 : * SPDX-License-Identifier: MIT
11 : ****************************************************************************/
12 :
13 : #include "gdal_frmts.h"
14 : #include "gdal_pam.h"
15 : #include "gdal_driver.h"
16 : #include "gdal_drivermanager.h"
17 : #include "gdal_openinfo.h"
18 : #include "gdal_cpp_functions.h"
19 : #include "ogr_spatialref.h"
20 :
21 : constexpr int HEADER_LINE_COUNT = 5;
22 :
23 : typedef struct
24 : {
25 : int nCode;
26 : const char *pszDesc;
27 : } LULCDescStruct;
28 :
29 : static const LULCDescStruct asLULCDesc[] = {
30 : {1, "Urban or Built-Up Land"},
31 : {2, "Agricultural Land"},
32 : {3, "Rangeland"},
33 : {4, "Forest Land"},
34 : {5, "Water"},
35 : {6, "Wetland"},
36 : {7, "Barren Land"},
37 : {8, "Tundra"},
38 : {9, "Perennial Snow and Ice"},
39 : {11, "Residential"},
40 : {12, "Commercial Services"},
41 : {13, "Industrial"},
42 : {14, "Transportation, Communications"},
43 : {15, "Industrial and Commercial"},
44 : {16, "Mixed Urban or Built-Up Land"},
45 : {17, "Other Urban or Built-Up Land"},
46 : {21, "Cropland and Pasture"},
47 : {22, "Orchards, Groves, Vineyards, Nurseries"},
48 : {23, "Confined Feeding Operations"},
49 : {24, "Other Agricultural Land"},
50 : {31, "Herbaceous Rangeland"},
51 : {32, "Shrub and Brush Rangeland"},
52 : {33, "Mixed Rangeland"},
53 : {41, "Deciduous Forest Land"},
54 : {42, "Evergreen Forest Land"},
55 : {43, "Mixed Forest Land"},
56 : {51, "Streams and Canals"},
57 : {52, "Lakes"},
58 : {53, "Reservoirs"},
59 : {54, "Bays and Estuaries"},
60 : {61, "Forested Wetlands"},
61 : {62, "Nonforested Wetlands"},
62 : {71, "Dry Salt Flats"},
63 : {72, "Beaches"},
64 : {73, "Sandy Areas Other than Beaches"},
65 : {74, "Bare Exposed Rock"},
66 : {75, "Strip Mines, Quarries, and Gravel Pits"},
67 : {76, "Transitional Areas"},
68 : {77, "Mixed Barren Land"},
69 : {81, "Shrub and Brush Tundra"},
70 : {82, "Herbaceous Tundra"},
71 : {83, "Bare Ground"},
72 : {84, "Wet Tundra"},
73 : {85, "Mixed Tundra"},
74 : {91, "Perennial Snowfields"},
75 : {92, "Glaciers"}};
76 :
77 : static const char *const apszBandDescription[] = {
78 : "Land Use and Land Cover",
79 : "Political units",
80 : "Census county subdivisions and SMSA tracts",
81 : "Hydrologic units",
82 : "Federal land ownership",
83 : "State land ownership"};
84 :
85 : /************************************************************************/
86 : /* ==================================================================== */
87 : /* CTGDataset */
88 : /* ==================================================================== */
89 : /************************************************************************/
90 :
91 : class CTGRasterBand;
92 :
93 : class CTGDataset final : public GDALPamDataset
94 : {
95 : friend class CTGRasterBand;
96 :
97 : VSILFILE *fp;
98 :
99 : int nNWEasting, nNWNorthing, nCellSize, nUTMZone;
100 : OGRSpatialReference m_oSRS{};
101 :
102 : int bHasReadImagery;
103 : GByte *pabyImage;
104 :
105 : int ReadImagery();
106 :
107 : static const char *ExtractField(char *szOutput, const char *pszBuffer,
108 : int nOffset, int nLength);
109 :
110 : public:
111 : CTGDataset();
112 : ~CTGDataset() override;
113 :
114 : CPLErr GetGeoTransform(GDALGeoTransform >) const override;
115 :
116 1 : const OGRSpatialReference *GetSpatialRef() const override
117 : {
118 1 : return &m_oSRS;
119 : }
120 :
121 : static GDALDataset *Open(GDALOpenInfo *);
122 : static int Identify(GDALOpenInfo *);
123 : };
124 :
125 : /************************************************************************/
126 : /* ==================================================================== */
127 : /* CTGRasterBand */
128 : /* ==================================================================== */
129 : /************************************************************************/
130 :
131 : class CTGRasterBand final : public GDALPamRasterBand
132 : {
133 : friend class CTGDataset;
134 :
135 : char **papszCategories;
136 :
137 : public:
138 : CTGRasterBand(CTGDataset *, int);
139 : ~CTGRasterBand() override;
140 :
141 : CPLErr IReadBlock(int, int, void *) override;
142 : double GetNoDataValue(int *pbSuccess = nullptr) override;
143 : char **GetCategoryNames() override;
144 : };
145 :
146 : /************************************************************************/
147 : /* CTGRasterBand() */
148 : /************************************************************************/
149 :
150 18 : CTGRasterBand::CTGRasterBand(CTGDataset *poDSIn, int nBandIn)
151 18 : : papszCategories(nullptr)
152 : {
153 18 : poDS = poDSIn;
154 18 : nBand = nBandIn;
155 :
156 18 : eDataType = GDT_Int32;
157 :
158 18 : nBlockXSize = poDS->GetRasterXSize();
159 18 : nBlockYSize = poDS->GetRasterYSize();
160 18 : }
161 :
162 : /************************************************************************/
163 : /* ~CTGRasterBand() */
164 : /************************************************************************/
165 :
166 36 : CTGRasterBand::~CTGRasterBand()
167 :
168 : {
169 18 : CSLDestroy(papszCategories);
170 36 : }
171 :
172 : /************************************************************************/
173 : /* IReadBlock() */
174 : /************************************************************************/
175 :
176 1 : CPLErr CTGRasterBand::IReadBlock(int /* nBlockXOff */, int /* nBlockYOff */,
177 : void *pImage)
178 : {
179 1 : CTGDataset *poGDS = cpl::down_cast<CTGDataset *>(poDS);
180 :
181 1 : poGDS->ReadImagery();
182 1 : memcpy(pImage,
183 1 : poGDS->pabyImage +
184 1 : sizeof(int) * (nBand - 1) * nBlockXSize * nBlockYSize,
185 1 : sizeof(int) * nBlockXSize * nBlockYSize);
186 :
187 1 : return CE_None;
188 : }
189 :
190 : /************************************************************************/
191 : /* GetNoDataValue() */
192 : /************************************************************************/
193 :
194 1 : double CTGRasterBand::GetNoDataValue(int *pbSuccess)
195 : {
196 1 : if (pbSuccess)
197 1 : *pbSuccess = TRUE;
198 :
199 1 : return 0.0;
200 : }
201 :
202 : /************************************************************************/
203 : /* GetCategoryNames() */
204 : /************************************************************************/
205 :
206 2 : char **CTGRasterBand::GetCategoryNames()
207 : {
208 2 : if (nBand != 1)
209 1 : return nullptr;
210 :
211 1 : if (papszCategories != nullptr)
212 0 : return papszCategories;
213 :
214 1 : int nasLULCDescSize = (int)(sizeof(asLULCDesc) / sizeof(asLULCDesc[0]));
215 1 : int nCategoriesSize = asLULCDesc[nasLULCDescSize - 1].nCode;
216 1 : papszCategories = (char **)CPLCalloc(nCategoriesSize + 2, sizeof(char *));
217 47 : for (int i = 0; i < nasLULCDescSize; i++)
218 : {
219 46 : papszCategories[asLULCDesc[i].nCode] = CPLStrdup(asLULCDesc[i].pszDesc);
220 : }
221 93 : for (int i = 0; i < nCategoriesSize; i++)
222 : {
223 92 : if (papszCategories[i] == nullptr)
224 47 : papszCategories[i] = CPLStrdup("");
225 : }
226 1 : papszCategories[nCategoriesSize + 1] = nullptr;
227 :
228 1 : return papszCategories;
229 : }
230 :
231 : /************************************************************************/
232 : /* ~CTGDataset() */
233 : /************************************************************************/
234 :
235 3 : CTGDataset::CTGDataset()
236 : : fp(nullptr), nNWEasting(0), nNWNorthing(0), nCellSize(0), nUTMZone(0),
237 3 : bHasReadImagery(FALSE), pabyImage(nullptr)
238 : {
239 3 : m_oSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
240 3 : }
241 :
242 : /************************************************************************/
243 : /* ~CTGDataset() */
244 : /************************************************************************/
245 :
246 6 : CTGDataset::~CTGDataset()
247 :
248 : {
249 3 : CPLFree(pabyImage);
250 3 : if (fp != nullptr)
251 3 : VSIFCloseL(fp);
252 6 : }
253 :
254 : /************************************************************************/
255 : /* ExtractField() */
256 : /************************************************************************/
257 :
258 69 : const char *CTGDataset::ExtractField(char *szField, const char *pszBuffer,
259 : int nOffset, int nLength)
260 : {
261 69 : CPLAssert(nLength <= 10);
262 69 : memcpy(szField, pszBuffer + nOffset, nLength);
263 69 : szField[nLength] = 0;
264 69 : return szField;
265 : }
266 :
267 : /************************************************************************/
268 : /* ReadImagery() */
269 : /************************************************************************/
270 :
271 1 : int CTGDataset::ReadImagery()
272 : {
273 1 : if (bHasReadImagery)
274 0 : return TRUE;
275 :
276 1 : bHasReadImagery = TRUE;
277 :
278 : char szLine[81];
279 : char szField[11];
280 1 : szLine[80] = 0;
281 1 : int nLine = HEADER_LINE_COUNT;
282 1 : VSIFSeekL(fp, nLine * 80, SEEK_SET);
283 1 : const int nCells = nRasterXSize * nRasterYSize;
284 2 : while (VSIFReadL(szLine, 1, 80, fp) == 80)
285 : {
286 1 : const int nZone = atoi(ExtractField(szField, szLine, 0, 3));
287 1 : if (nZone != nUTMZone)
288 : {
289 0 : CPLError(CE_Failure, CPLE_AppDefined,
290 : "Read error at line %d, %s. Did not expected UTM zone %d",
291 : nLine, szLine, nZone);
292 0 : return FALSE;
293 : }
294 : const int nX =
295 1 : atoi(ExtractField(szField, szLine, 3, 8)) - nCellSize / 2;
296 : const int nY =
297 1 : atoi(ExtractField(szField, szLine, 11, 8)) + nCellSize / 2;
298 1 : const GIntBig nDiffX = static_cast<GIntBig>(nX) - nNWEasting;
299 1 : const GIntBig nDiffY = static_cast<GIntBig>(nNWNorthing) - nY;
300 1 : if (nDiffX < 0 || (nDiffX % nCellSize) != 0 || nDiffY < 0 ||
301 1 : (nDiffY % nCellSize) != 0)
302 : {
303 0 : CPLError(CE_Failure, CPLE_AppDefined,
304 : "Read error at line %d, %s. Unexpected cell coordinates",
305 : nLine, szLine);
306 0 : return FALSE;
307 : }
308 1 : const GIntBig nCellX = nDiffX / nCellSize;
309 1 : const GIntBig nCellY = nDiffY / nCellSize;
310 1 : if (nCellX >= nRasterXSize || nCellY >= nRasterYSize)
311 : {
312 0 : CPLError(CE_Failure, CPLE_AppDefined,
313 : "Read error at line %d, %s. Unexpected cell coordinates",
314 : nLine, szLine);
315 0 : return FALSE;
316 : }
317 7 : for (int i = 0; i < 6; i++)
318 : {
319 6 : int nVal = atoi(ExtractField(szField, szLine, 20 + 10 * i, 10));
320 6 : if (nVal >= 2000000000)
321 0 : nVal = 0;
322 : reinterpret_cast<int *>(
323 6 : pabyImage)[i * nCells +
324 6 : static_cast<int>(nCellY) * nRasterXSize + nCellX] =
325 : nVal;
326 : }
327 :
328 1 : nLine++;
329 : }
330 :
331 1 : return TRUE;
332 : }
333 :
334 : /************************************************************************/
335 : /* Identify() */
336 : /************************************************************************/
337 :
338 58301 : int CTGDataset::Identify(GDALOpenInfo *poOpenInfo)
339 : {
340 116603 : CPLString osFilename; // let in that scope
341 :
342 58301 : GDALOpenInfo *poOpenInfoToDelete = nullptr;
343 : /* GZipped grid_cell.gz files are common, so automagically open them */
344 : /* if the /vsigzip/ has not been explicitly passed */
345 58301 : const char *pszFilename = CPLGetFilename(poOpenInfo->pszFilename);
346 58302 : if ((EQUAL(pszFilename, "grid_cell.gz") ||
347 58302 : EQUAL(pszFilename, "grid_cell1.gz") ||
348 58302 : EQUAL(pszFilename, "grid_cell2.gz")) &&
349 0 : !STARTS_WITH_CI(poOpenInfo->pszFilename, "/vsigzip/"))
350 : {
351 0 : osFilename = "/vsigzip/";
352 0 : osFilename += poOpenInfo->pszFilename;
353 0 : poOpenInfo = poOpenInfoToDelete = new GDALOpenInfo(
354 0 : osFilename.c_str(), GA_ReadOnly, poOpenInfo->GetSiblingFiles());
355 : }
356 :
357 58302 : if (poOpenInfo->nHeaderBytes < HEADER_LINE_COUNT * 80)
358 : {
359 55830 : delete poOpenInfoToDelete;
360 55830 : return FALSE;
361 : }
362 :
363 : /* -------------------------------------------------------------------- */
364 : /* Check that it looks roughly as a CTG dataset */
365 : /* -------------------------------------------------------------------- */
366 2472 : const char *pszData = (const char *)poOpenInfo->pabyHeader;
367 5100 : for (int i = 0; i < 4 * 80; i++)
368 : {
369 5093 : if (!((pszData[i] >= '0' && pszData[i] <= '9') || pszData[i] == ' ' ||
370 2467 : pszData[i] == '-'))
371 : {
372 2465 : delete poOpenInfoToDelete;
373 2465 : return FALSE;
374 : }
375 : }
376 :
377 : char szField[11];
378 7 : int nRows = atoi(ExtractField(szField, pszData, 0, 10));
379 7 : int nCols = atoi(ExtractField(szField, pszData, 20, 10));
380 7 : int nMinColIndex = atoi(ExtractField(szField, pszData + 80, 0, 5));
381 7 : int nMinRowIndex = atoi(ExtractField(szField, pszData + 80, 5, 5));
382 7 : int nMaxColIndex = atoi(ExtractField(szField, pszData + 80, 10, 5));
383 7 : int nMaxRowIndex = atoi(ExtractField(szField, pszData + 80, 15, 5));
384 :
385 7 : if (nRows <= 0 || nCols <= 0 || nMinColIndex != 1 || nMinRowIndex != 1 ||
386 6 : nMaxRowIndex != nRows || nMaxColIndex != nCols)
387 : {
388 1 : delete poOpenInfoToDelete;
389 1 : return FALSE;
390 : }
391 :
392 6 : delete poOpenInfoToDelete;
393 6 : return TRUE;
394 : }
395 :
396 : /************************************************************************/
397 : /* Open() */
398 : /************************************************************************/
399 :
400 3 : GDALDataset *CTGDataset::Open(GDALOpenInfo *poOpenInfo)
401 :
402 : {
403 3 : if (!Identify(poOpenInfo))
404 0 : return nullptr;
405 :
406 6 : CPLString osFilename(poOpenInfo->pszFilename);
407 :
408 : /* GZipped grid_cell.gz files are common, so automagically open them */
409 : /* if the /vsigzip/ has not been explicitly passed */
410 3 : const char *pszFilename = CPLGetFilename(poOpenInfo->pszFilename);
411 3 : if ((EQUAL(pszFilename, "grid_cell.gz") ||
412 3 : EQUAL(pszFilename, "grid_cell1.gz") ||
413 3 : EQUAL(pszFilename, "grid_cell2.gz")) &&
414 0 : !STARTS_WITH_CI(poOpenInfo->pszFilename, "/vsigzip/"))
415 : {
416 0 : osFilename = "/vsigzip/";
417 0 : osFilename += poOpenInfo->pszFilename;
418 : }
419 :
420 3 : if (poOpenInfo->eAccess == GA_Update)
421 : {
422 0 : ReportUpdateNotSupportedByDriver("CTG");
423 0 : return nullptr;
424 : }
425 :
426 : /* -------------------------------------------------------------------- */
427 : /* Find dataset characteristics */
428 : /* -------------------------------------------------------------------- */
429 3 : VSILFILE *fp = VSIFOpenL(osFilename.c_str(), "rb");
430 3 : if (fp == nullptr)
431 0 : return nullptr;
432 :
433 : char szHeader[HEADER_LINE_COUNT * 80 + 1];
434 3 : szHeader[HEADER_LINE_COUNT * 80] = 0;
435 3 : if (VSIFReadL(szHeader, 1, HEADER_LINE_COUNT * 80, fp) !=
436 : HEADER_LINE_COUNT * 80)
437 : {
438 0 : VSIFCloseL(fp);
439 0 : return nullptr;
440 : }
441 :
442 216 : for (int i = HEADER_LINE_COUNT * 80 - 1; i >= 0; i--)
443 : {
444 216 : if (szHeader[i] == ' ')
445 213 : szHeader[i] = 0;
446 : else
447 3 : break;
448 : }
449 :
450 : char szField[11];
451 3 : int nRows = atoi(ExtractField(szField, szHeader, 0, 10));
452 3 : int nCols = atoi(ExtractField(szField, szHeader, 20, 10));
453 :
454 : /* -------------------------------------------------------------------- */
455 : /* Create a corresponding GDALDataset. */
456 : /* -------------------------------------------------------------------- */
457 3 : CTGDataset *poDS = new CTGDataset();
458 3 : poDS->fp = fp;
459 3 : fp = nullptr;
460 3 : poDS->nRasterXSize = nCols;
461 3 : poDS->nRasterYSize = nRows;
462 :
463 3 : poDS->SetMetadataItem("TITLE", szHeader + 4 * 80);
464 :
465 3 : poDS->nCellSize = atoi(ExtractField(szField, szHeader, 35, 5));
466 3 : if (poDS->nCellSize <= 0 || poDS->nCellSize >= 10000)
467 : {
468 0 : delete poDS;
469 0 : return nullptr;
470 : }
471 3 : poDS->nNWEasting = atoi(ExtractField(szField, szHeader + 3 * 80, 40, 10));
472 3 : poDS->nNWNorthing = atoi(ExtractField(szField, szHeader + 3 * 80, 50, 10));
473 3 : poDS->nUTMZone = atoi(ExtractField(szField, szHeader, 50, 5));
474 3 : if (poDS->nUTMZone <= 0 || poDS->nUTMZone > 60)
475 : {
476 0 : delete poDS;
477 0 : return nullptr;
478 : }
479 :
480 3 : poDS->m_oSRS.importFromEPSG(32600 + poDS->nUTMZone);
481 :
482 3 : if (!GDALCheckDatasetDimensions(poDS->nRasterXSize, poDS->nRasterYSize))
483 : {
484 0 : delete poDS;
485 0 : return nullptr;
486 : }
487 :
488 : /* -------------------------------------------------------------------- */
489 : /* Read the imagery */
490 : /* -------------------------------------------------------------------- */
491 : GByte *pabyImage =
492 3 : (GByte *)VSICalloc(static_cast<size_t>(nCols) * nRows, 6 * sizeof(int));
493 3 : if (pabyImage == nullptr)
494 : {
495 0 : delete poDS;
496 0 : return nullptr;
497 : }
498 3 : poDS->pabyImage = pabyImage;
499 :
500 : /* -------------------------------------------------------------------- */
501 : /* Create band information objects. */
502 : /* -------------------------------------------------------------------- */
503 3 : poDS->nBands = 6;
504 21 : for (int i = 0; i < poDS->nBands; i++)
505 : {
506 18 : poDS->SetBand(i + 1, new CTGRasterBand(poDS, i + 1));
507 18 : poDS->GetRasterBand(i + 1)->SetDescription(apszBandDescription[i]);
508 : }
509 :
510 : /* -------------------------------------------------------------------- */
511 : /* Initialize any PAM information. */
512 : /* -------------------------------------------------------------------- */
513 3 : poDS->SetDescription(poOpenInfo->pszFilename);
514 3 : poDS->TryLoadXML();
515 :
516 : /* -------------------------------------------------------------------- */
517 : /* Support overviews. */
518 : /* -------------------------------------------------------------------- */
519 3 : poDS->oOvManager.Initialize(poDS, poOpenInfo->pszFilename);
520 :
521 3 : return poDS;
522 : }
523 :
524 : /************************************************************************/
525 : /* GetGeoTransform() */
526 : /************************************************************************/
527 :
528 1 : CPLErr CTGDataset::GetGeoTransform(GDALGeoTransform >) const
529 :
530 : {
531 1 : gt[0] = static_cast<double>(nNWEasting) - nCellSize / 2;
532 1 : gt[1] = nCellSize;
533 1 : gt[2] = 0;
534 1 : gt[3] = static_cast<double>(nNWNorthing) + nCellSize / 2;
535 1 : gt[4] = 0.;
536 1 : gt[5] = -nCellSize;
537 :
538 1 : return CE_None;
539 : }
540 :
541 : /************************************************************************/
542 : /* GDALRegister_CTG() */
543 : /************************************************************************/
544 :
545 2038 : void GDALRegister_CTG()
546 :
547 : {
548 2038 : if (GDALGetDriverByName("CTG") != nullptr)
549 283 : return;
550 :
551 1755 : GDALDriver *poDriver = new GDALDriver();
552 :
553 1755 : poDriver->SetDescription("CTG");
554 1755 : poDriver->SetMetadataItem(GDAL_DCAP_RASTER, "YES");
555 1755 : poDriver->SetMetadataItem(GDAL_DMD_LONGNAME,
556 1755 : "USGS LULC Composite Theme Grid");
557 1755 : poDriver->SetMetadataItem(GDAL_DMD_HELPTOPIC, "drivers/raster/ctg.html");
558 :
559 1755 : poDriver->SetMetadataItem(GDAL_DCAP_VIRTUALIO, "YES");
560 :
561 1755 : poDriver->pfnOpen = CTGDataset::Open;
562 1755 : poDriver->pfnIdentify = CTGDataset::Identify;
563 :
564 1755 : GetGDALDriverManager()->RegisterDriver(poDriver);
565 : }
|
