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