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