1 : /******************************************************************************
2 : * $Id: ctgdataset.cpp 21774 2011-02-21 19:33:56Z 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
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 21774 2011-02-21 19:33:56Z 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 : char szField[11];
122 :
123 : int nNWEasting, nNWNorthing, nCellSize, nUTMZone;
124 : char *pszProjection;
125 :
126 : int bHasReadImagery;
127 : GByte *pabyImage;
128 :
129 : int ReadImagery();
130 :
131 : static const char* ExtractField(char* szOutput, const char* pszBuffer,
132 : int nOffset, int nLength);
133 :
134 : public:
135 : CTGDataset();
136 : virtual ~CTGDataset();
137 :
138 : virtual CPLErr GetGeoTransform( double * );
139 : virtual const char* GetProjectionRef();
140 :
141 : static GDALDataset *Open( GDALOpenInfo * );
142 : static int Identify( GDALOpenInfo * );
143 : };
144 :
145 : /************************************************************************/
146 : /* ==================================================================== */
147 : /* CTGRasterBand */
148 : /* ==================================================================== */
149 : /************************************************************************/
150 :
151 : class CTGRasterBand : public GDALPamRasterBand
152 : {
153 : friend class CTGDataset;
154 :
155 : char** papszCategories;
156 :
157 : public:
158 :
159 : CTGRasterBand( CTGDataset *, int );
160 : ~CTGRasterBand();
161 :
162 : virtual CPLErr IReadBlock( int, int, void * );
163 : virtual double GetNoDataValue( int *pbSuccess = NULL );
164 : virtual char **GetCategoryNames();
165 : };
166 :
167 :
168 : /************************************************************************/
169 : /* CTGRasterBand() */
170 : /************************************************************************/
171 :
172 24 : CTGRasterBand::CTGRasterBand( CTGDataset *poDS, int nBand )
173 :
174 : {
175 24 : this->poDS = poDS;
176 24 : this->nBand = nBand;
177 :
178 24 : eDataType = GDT_Int32;
179 :
180 24 : nBlockXSize = poDS->GetRasterXSize();
181 24 : nBlockYSize = poDS->GetRasterYSize();
182 :
183 24 : papszCategories = NULL;
184 24 : }
185 :
186 : /************************************************************************/
187 : /* ~CTGRasterBand() */
188 : /************************************************************************/
189 :
190 24 : CTGRasterBand::~CTGRasterBand()
191 :
192 : {
193 24 : CSLDestroy(papszCategories);
194 24 : }
195 : /************************************************************************/
196 : /* IReadBlock() */
197 : /************************************************************************/
198 :
199 2 : CPLErr CTGRasterBand::IReadBlock( int nBlockXOff, int nBlockYOff,
200 : void * pImage )
201 :
202 : {
203 2 : CTGDataset* poGDS = (CTGDataset* ) poDS;
204 :
205 2 : poGDS->ReadImagery();
206 : memcpy(pImage,
207 : poGDS->pabyImage + (nBand - 1) * nBlockXSize * nBlockYSize * sizeof(int),
208 2 : nBlockXSize * nBlockYSize * sizeof(int));
209 :
210 2 : return CE_None;
211 : }
212 :
213 : /************************************************************************/
214 : /* GetNoDataValue() */
215 : /************************************************************************/
216 :
217 2 : double CTGRasterBand::GetNoDataValue( int *pbSuccess )
218 : {
219 2 : if (pbSuccess)
220 2 : *pbSuccess = TRUE;
221 :
222 2 : return 0.;
223 : }
224 :
225 : /************************************************************************/
226 : /* GetCategoryNames() */
227 : /************************************************************************/
228 :
229 4 : char **CTGRasterBand::GetCategoryNames()
230 : {
231 4 : if (nBand != 1)
232 2 : return NULL;
233 :
234 2 : if (papszCategories != NULL)
235 0 : return papszCategories;
236 :
237 : int i;
238 2 : int nasLULCDescSize = (int)(sizeof(asLULCDesc) / sizeof(asLULCDesc[0]));
239 2 : int nCategoriesSize = asLULCDesc[nasLULCDescSize - 1].nCode;
240 2 : papszCategories = (char**)CPLCalloc(nCategoriesSize + 2, sizeof(char*));
241 94 : for(i=0;i<nasLULCDescSize;i++)
242 : {
243 92 : papszCategories[asLULCDesc[i].nCode] = CPLStrdup(asLULCDesc[i].pszDesc);
244 : }
245 186 : for(i=0;i<nCategoriesSize;i++)
246 : {
247 184 : if (papszCategories[i] == NULL)
248 94 : papszCategories[i] = CPLStrdup("");
249 : }
250 2 : papszCategories[nCategoriesSize + 1] = NULL;
251 :
252 2 : return papszCategories;
253 : }
254 :
255 : /************************************************************************/
256 : /* ~CTGDataset() */
257 : /************************************************************************/
258 :
259 4 : CTGDataset::CTGDataset()
260 : {
261 4 : pszProjection = NULL;
262 4 : bHasReadImagery = FALSE;
263 4 : pabyImage = NULL;
264 4 : fp = NULL;
265 4 : }
266 :
267 : /************************************************************************/
268 : /* ~CTGDataset() */
269 : /************************************************************************/
270 :
271 4 : CTGDataset::~CTGDataset()
272 :
273 : {
274 4 : CPLFree(pszProjection);
275 4 : CPLFree(pabyImage);
276 4 : if (fp != NULL)
277 4 : VSIFCloseL(fp);
278 4 : }
279 :
280 : /************************************************************************/
281 : /* ExtractField() */
282 : /************************************************************************/
283 :
284 66 : const char* CTGDataset::ExtractField(char* szField, const char* pszBuffer,
285 : int nOffset, int nLength)
286 : {
287 66 : CPLAssert(nLength <= 10);
288 66 : memcpy(szField, pszBuffer + nOffset, nLength);
289 66 : szField[nLength] = 0;
290 66 : return szField;
291 : }
292 :
293 : /************************************************************************/
294 : /* ReadImagery() */
295 : /************************************************************************/
296 :
297 2 : int CTGDataset::ReadImagery()
298 : {
299 2 : if (bHasReadImagery)
300 0 : return TRUE;
301 :
302 2 : bHasReadImagery = TRUE;
303 :
304 : char szLine[81];
305 : char szField[11];
306 2 : szLine[80] = 0;
307 2 : int nLine = HEADER_LINE_COUNT;
308 2 : VSIFSeekL(fp, nLine * 80, SEEK_SET);
309 2 : int nCells = nRasterXSize * nRasterYSize;
310 6 : while(VSIFReadL(szLine, 1, 80, fp) == 80)
311 : {
312 2 : int nZone = atoi(ExtractField(szField, szLine, 0, 3));
313 2 : if (nZone != nUTMZone)
314 : {
315 : CPLError(CE_Failure, CPLE_AppDefined,
316 : "Read error at line %d, %s. Did not expected UTM zone %d",
317 0 : nLine, szLine, nZone);
318 0 : return FALSE;
319 : }
320 2 : int nX = atoi(ExtractField(szField, szLine, 3, 8)) - nCellSize / 2;
321 2 : int nY = atoi(ExtractField(szField, szLine, 11, 8)) + nCellSize / 2;
322 2 : int nDiffX = nX - nNWEasting;
323 2 : int nDiffY = nNWNorthing - nY;
324 2 : if (nDiffX < 0 || (nDiffX % nCellSize) != 0 ||
325 : nDiffY < 0 || (nDiffY % nCellSize) != 0)
326 : {
327 : CPLError(CE_Failure, CPLE_AppDefined,
328 : "Read error at line %d, %s. Unexpected cell coordinates",
329 0 : nLine, szLine);
330 0 : return FALSE;
331 : }
332 2 : int nCellX = nDiffX / nCellSize;
333 2 : int nCellY = nDiffY / nCellSize;
334 2 : if (nCellX >= nRasterXSize || nCellY >= nRasterYSize)
335 : {
336 : CPLError(CE_Failure, CPLE_AppDefined,
337 : "Read error at line %d, %s. Unexpected cell coordinates",
338 0 : nLine, szLine);
339 0 : return FALSE;
340 : }
341 : int i;
342 14 : for(i=0;i<6;i++)
343 : {
344 12 : int nVal = atoi(ExtractField(szField, szLine, 20 + 10*i, 10));
345 12 : if (nVal >= 2000000000)
346 0 : nVal = 0;
347 12 : ((int*)pabyImage)[i * nCells + nCellY * nRasterXSize + nCellX] = nVal;
348 : }
349 :
350 2 : nLine ++;
351 : }
352 :
353 2 : return TRUE;
354 : }
355 :
356 : /************************************************************************/
357 : /* Identify() */
358 : /************************************************************************/
359 :
360 21502 : int CTGDataset::Identify( GDALOpenInfo * poOpenInfo )
361 : {
362 21502 : CPLString osFilename(poOpenInfo->pszFilename);
363 :
364 21502 : GDALOpenInfo* poOpenInfoToDelete = NULL;
365 : /* GZipped grid_cell.gz files are common, so automagically open them */
366 : /* if the /vsigzip/ has not been explicitely passed */
367 21502 : const char* pszFilename = CPLGetFilename(poOpenInfo->pszFilename);
368 21502 : if ((EQUAL(pszFilename, "grid_cell.gz") ||
369 : EQUAL(pszFilename, "grid_cell1.gz") ||
370 : EQUAL(pszFilename, "grid_cell2.gz")) &&
371 : !EQUALN(poOpenInfo->pszFilename, "/vsigzip/", 9))
372 : {
373 0 : osFilename = "/vsigzip/";
374 0 : osFilename += poOpenInfo->pszFilename;
375 : poOpenInfo = poOpenInfoToDelete =
376 : new GDALOpenInfo(osFilename.c_str(), GA_ReadOnly,
377 0 : poOpenInfo->papszSiblingFiles);
378 : }
379 :
380 21502 : if (poOpenInfo->nHeaderBytes < HEADER_LINE_COUNT * 80)
381 : {
382 21208 : delete poOpenInfoToDelete;
383 21208 : return FALSE;
384 : }
385 :
386 : /* -------------------------------------------------------------------- */
387 : /* Chech that it looks roughly as a CTG dataset */
388 : /* -------------------------------------------------------------------- */
389 294 : const char* pszData = (const char*)poOpenInfo->pabyHeader;
390 : int i;
391 1574 : for(i=0;i<4 * 80;i++)
392 : {
393 2794 : if (!((pszData[i] >= '0' && pszData[i] <= '9') ||
394 1224 : pszData[i] == ' ' || pszData[i] == '-'))
395 : {
396 290 : delete poOpenInfoToDelete;
397 290 : return FALSE;
398 : }
399 : }
400 :
401 : char szField[11];
402 4 : int nRows = atoi(ExtractField(szField, pszData, 0, 10));
403 4 : int nCols = atoi(ExtractField(szField, pszData, 20, 10));
404 4 : int nMinColIndex = atoi(ExtractField(szField, pszData+80, 0, 5));
405 4 : int nMinRowIndex = atoi(ExtractField(szField, pszData+80, 5, 5));
406 4 : int nMaxColIndex = atoi(ExtractField(szField, pszData+80, 10, 5));
407 4 : int nMaxRowIndex = atoi(ExtractField(szField, pszData+80, 15, 5));
408 :
409 4 : if (nRows <= 0 || nCols <= 0 ||
410 : nMinColIndex != 1 || nMinRowIndex != 1 ||
411 : nMaxRowIndex != nRows || nMaxColIndex != nCols)
412 : {
413 0 : delete poOpenInfoToDelete;
414 0 : return FALSE;
415 : }
416 :
417 4 : delete poOpenInfoToDelete;
418 4 : return TRUE;
419 : }
420 :
421 :
422 : /************************************************************************/
423 : /* Open() */
424 : /************************************************************************/
425 :
426 2608 : GDALDataset *CTGDataset::Open( GDALOpenInfo * poOpenInfo )
427 :
428 : {
429 : int i;
430 :
431 2608 : if (!Identify(poOpenInfo))
432 2604 : return NULL;
433 :
434 4 : CPLString osFilename(poOpenInfo->pszFilename);
435 :
436 : /* GZipped grid_cell.gz files are common, so automagically open them */
437 : /* if the /vsigzip/ has not been explicitely passed */
438 4 : const char* pszFilename = CPLGetFilename(poOpenInfo->pszFilename);
439 4 : if ((EQUAL(pszFilename, "grid_cell.gz") ||
440 : EQUAL(pszFilename, "grid_cell1.gz") ||
441 : EQUAL(pszFilename, "grid_cell2.gz")) &&
442 : !EQUALN(poOpenInfo->pszFilename, "/vsigzip/", 9))
443 : {
444 0 : osFilename = "/vsigzip/";
445 0 : osFilename += poOpenInfo->pszFilename;
446 : }
447 :
448 4 : if (poOpenInfo->eAccess == GA_Update)
449 : {
450 : CPLError( CE_Failure, CPLE_NotSupported,
451 : "The CTG driver does not support update access to existing"
452 0 : " datasets.\n" );
453 0 : return NULL;
454 : }
455 :
456 : /* -------------------------------------------------------------------- */
457 : /* Find dataset characteristics */
458 : /* -------------------------------------------------------------------- */
459 4 : VSILFILE* fp = VSIFOpenL(osFilename.c_str(), "rb");
460 4 : if (fp == NULL)
461 0 : return NULL;
462 :
463 : char szHeader[HEADER_LINE_COUNT * 80+1];
464 4 : szHeader[HEADER_LINE_COUNT * 80] = 0;
465 4 : if (VSIFReadL(szHeader, 1, HEADER_LINE_COUNT * 80, fp) != HEADER_LINE_COUNT * 80)
466 : {
467 0 : VSIFCloseL(fp);
468 0 : return NULL;
469 : }
470 :
471 576 : for(i=HEADER_LINE_COUNT * 80 - 1;i>=0;i--)
472 : {
473 288 : if (szHeader[i] == ' ')
474 284 : szHeader[i] = 0;
475 : else
476 4 : break;
477 : }
478 :
479 : char szField[11];
480 4 : int nRows = atoi(ExtractField(szField, szHeader, 0, 10));
481 4 : int nCols = atoi(ExtractField(szField, szHeader, 20, 10));
482 :
483 : /* -------------------------------------------------------------------- */
484 : /* Create a corresponding GDALDataset. */
485 : /* -------------------------------------------------------------------- */
486 : CTGDataset *poDS;
487 :
488 4 : poDS = new CTGDataset();
489 4 : poDS->fp = fp;
490 4 : fp = NULL;
491 4 : poDS->nRasterXSize = nCols;
492 4 : poDS->nRasterYSize = nRows;
493 :
494 4 : poDS->SetMetadataItem("TITLE", szHeader + 4 * 80);
495 :
496 4 : poDS->nCellSize = atoi(ExtractField(szField, szHeader, 35, 5));
497 4 : if (poDS->nCellSize <= 0 || poDS->nCellSize >= 10000)
498 : {
499 0 : delete poDS;
500 0 : return NULL;
501 : }
502 4 : poDS->nNWEasting = atoi(ExtractField(szField, szHeader + 3*80, 40, 10));
503 4 : poDS->nNWNorthing = atoi(ExtractField(szField, szHeader + 3*80, 50, 10));
504 4 : poDS->nUTMZone = atoi(ExtractField(szField, szHeader, 50, 5));
505 4 : if (poDS->nUTMZone <= 0 || poDS->nUTMZone > 60)
506 : {
507 0 : delete poDS;
508 0 : return NULL;
509 : }
510 :
511 4 : OGRSpatialReference oSRS;
512 4 : oSRS.importFromEPSG(32600 + poDS->nUTMZone);
513 4 : oSRS.exportToWkt(&poDS->pszProjection);
514 :
515 4 : if (!GDALCheckDatasetDimensions(poDS->nRasterXSize, poDS->nRasterYSize))
516 : {
517 0 : delete poDS;
518 0 : return NULL;
519 : }
520 :
521 : /* -------------------------------------------------------------------- */
522 : /* Read the imagery */
523 : /* -------------------------------------------------------------------- */
524 4 : GByte* pabyImage = (GByte*)VSICalloc(nCols * nRows, 6 * sizeof(int));
525 4 : if (pabyImage == NULL)
526 : {
527 0 : delete poDS;
528 0 : return NULL;
529 : }
530 4 : poDS->pabyImage = pabyImage;
531 :
532 : /* -------------------------------------------------------------------- */
533 : /* Create band information objects. */
534 : /* -------------------------------------------------------------------- */
535 4 : poDS->nBands = 6;
536 28 : for( i = 0; i < poDS->nBands; i++ )
537 : {
538 24 : poDS->SetBand( i+1, new CTGRasterBand( poDS, i+1 ) );
539 24 : poDS->GetRasterBand(i+1)->SetDescription(apszBandDescription[i]);
540 : }
541 :
542 : /* -------------------------------------------------------------------- */
543 : /* Initialize any PAM information. */
544 : /* -------------------------------------------------------------------- */
545 4 : poDS->SetDescription( poOpenInfo->pszFilename );
546 4 : poDS->TryLoadXML();
547 :
548 : /* -------------------------------------------------------------------- */
549 : /* Support overviews. */
550 : /* -------------------------------------------------------------------- */
551 4 : poDS->oOvManager.Initialize( poDS, poOpenInfo->pszFilename );
552 :
553 4 : return( poDS );
554 : }
555 :
556 : /************************************************************************/
557 : /* GetGeoTransform() */
558 : /************************************************************************/
559 :
560 2 : CPLErr CTGDataset::GetGeoTransform( double * padfTransform )
561 :
562 : {
563 2 : padfTransform[0] = nNWEasting - nCellSize / 2;
564 2 : padfTransform[1] = nCellSize;
565 2 : padfTransform[2] = 0;
566 2 : padfTransform[3] = nNWNorthing + nCellSize / 2;
567 2 : padfTransform[4] = 0.;
568 2 : padfTransform[5] = -nCellSize;
569 :
570 2 : return( CE_None );
571 : }
572 :
573 : /************************************************************************/
574 : /* GetProjectionRef() */
575 : /************************************************************************/
576 :
577 2 : const char* CTGDataset::GetProjectionRef()
578 :
579 : {
580 2 : return pszProjection;
581 : }
582 :
583 : /************************************************************************/
584 : /* GDALRegister_CTG() */
585 : /************************************************************************/
586 :
587 1135 : void GDALRegister_CTG()
588 :
589 : {
590 : GDALDriver *poDriver;
591 :
592 1135 : if( GDALGetDriverByName( "CTG" ) == NULL )
593 : {
594 1093 : poDriver = new GDALDriver();
595 :
596 1093 : poDriver->SetDescription( "CTG" );
597 : poDriver->SetMetadataItem( GDAL_DMD_LONGNAME,
598 1093 : "USGS LULC Composite Theme Grid" );
599 : poDriver->SetMetadataItem( GDAL_DMD_HELPTOPIC,
600 1093 : "frmt_various.html#CTG" );
601 :
602 1093 : poDriver->SetMetadataItem( GDAL_DCAP_VIRTUALIO, "YES" );
603 :
604 1093 : poDriver->pfnOpen = CTGDataset::Open;
605 1093 : poDriver->pfnIdentify = CTGDataset::Identify;
606 :
607 1093 : GetGDALDriverManager()->RegisterDriver( poDriver );
608 : }
609 1135 : }
610 :
|
