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