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 12 : CTGRasterBand::CTGRasterBand( CTGDataset *poDS, int nBand )
173 :
174 : {
175 12 : this->poDS = poDS;
176 12 : this->nBand = nBand;
177 :
178 12 : eDataType = GDT_Int32;
179 :
180 12 : nBlockXSize = poDS->GetRasterXSize();
181 12 : nBlockYSize = poDS->GetRasterYSize();
182 :
183 12 : papszCategories = NULL;
184 12 : }
185 :
186 : /************************************************************************/
187 : /* ~CTGRasterBand() */
188 : /************************************************************************/
189 :
190 12 : CTGRasterBand::~CTGRasterBand()
191 :
192 : {
193 12 : CSLDestroy(papszCategories);
194 12 : }
195 : /************************************************************************/
196 : /* IReadBlock() */
197 : /************************************************************************/
198 :
199 1 : CPLErr CTGRasterBand::IReadBlock( int nBlockXOff, int nBlockYOff,
200 : void * pImage )
201 :
202 : {
203 1 : CTGDataset* poGDS = (CTGDataset* ) poDS;
204 :
205 1 : poGDS->ReadImagery();
206 : memcpy(pImage,
207 : poGDS->pabyImage + (nBand - 1) * nBlockXSize * nBlockYSize * sizeof(int),
208 1 : nBlockXSize * nBlockYSize * sizeof(int));
209 :
210 1 : return CE_None;
211 : }
212 :
213 : /************************************************************************/
214 : /* GetNoDataValue() */
215 : /************************************************************************/
216 :
217 1 : double CTGRasterBand::GetNoDataValue( int *pbSuccess )
218 : {
219 1 : if (pbSuccess)
220 1 : *pbSuccess = TRUE;
221 :
222 1 : return 0.;
223 : }
224 :
225 : /************************************************************************/
226 : /* GetCategoryNames() */
227 : /************************************************************************/
228 :
229 2 : char **CTGRasterBand::GetCategoryNames()
230 : {
231 2 : if (nBand != 1)
232 1 : return NULL;
233 :
234 1 : if (papszCategories != NULL)
235 0 : return papszCategories;
236 :
237 : int i;
238 1 : int nasLULCDescSize = (int)(sizeof(asLULCDesc) / sizeof(asLULCDesc[0]));
239 1 : int nCategoriesSize = asLULCDesc[nasLULCDescSize - 1].nCode;
240 1 : papszCategories = (char**)CPLCalloc(nCategoriesSize + 2, sizeof(char*));
241 47 : for(i=0;i<nasLULCDescSize;i++)
242 : {
243 46 : papszCategories[asLULCDesc[i].nCode] = CPLStrdup(asLULCDesc[i].pszDesc);
244 : }
245 93 : for(i=0;i<nCategoriesSize;i++)
246 : {
247 92 : if (papszCategories[i] == NULL)
248 47 : papszCategories[i] = CPLStrdup("");
249 : }
250 1 : papszCategories[nCategoriesSize + 1] = NULL;
251 :
252 1 : return papszCategories;
253 : }
254 :
255 : /************************************************************************/
256 : /* ~CTGDataset() */
257 : /************************************************************************/
258 :
259 2 : CTGDataset::CTGDataset()
260 : {
261 2 : pszProjection = NULL;
262 2 : bHasReadImagery = FALSE;
263 2 : pabyImage = NULL;
264 2 : fp = NULL;
265 2 : }
266 :
267 : /************************************************************************/
268 : /* ~CTGDataset() */
269 : /************************************************************************/
270 :
271 2 : CTGDataset::~CTGDataset()
272 :
273 : {
274 2 : CPLFree(pszProjection);
275 2 : CPLFree(pabyImage);
276 2 : if (fp != NULL)
277 2 : VSIFCloseL(fp);
278 2 : }
279 :
280 : /************************************************************************/
281 : /* ExtractField() */
282 : /************************************************************************/
283 :
284 33 : const char* CTGDataset::ExtractField(char* szField, const char* pszBuffer,
285 : int nOffset, int nLength)
286 : {
287 33 : CPLAssert(nLength <= 10);
288 33 : memcpy(szField, pszBuffer + nOffset, nLength);
289 33 : szField[nLength] = 0;
290 33 : return szField;
291 : }
292 :
293 : /************************************************************************/
294 : /* ReadImagery() */
295 : /************************************************************************/
296 :
297 1 : int CTGDataset::ReadImagery()
298 : {
299 1 : if (bHasReadImagery)
300 0 : return TRUE;
301 :
302 1 : bHasReadImagery = TRUE;
303 :
304 : char szLine[81];
305 : char szField[11];
306 1 : szLine[80] = 0;
307 1 : int nLine = HEADER_LINE_COUNT;
308 1 : VSIFSeekL(fp, nLine * 80, SEEK_SET);
309 1 : int nCells = nRasterXSize * nRasterYSize;
310 3 : while(VSIFReadL(szLine, 1, 80, fp) == 80)
311 : {
312 1 : int nZone = atoi(ExtractField(szField, szLine, 0, 3));
313 1 : 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 1 : int nX = atoi(ExtractField(szField, szLine, 3, 8)) - nCellSize / 2;
321 1 : int nY = atoi(ExtractField(szField, szLine, 11, 8)) + nCellSize / 2;
322 1 : int nDiffX = nX - nNWEasting;
323 1 : int nDiffY = nNWNorthing - nY;
324 1 : 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 1 : int nCellX = nDiffX / nCellSize;
333 1 : int nCellY = nDiffY / nCellSize;
334 1 : 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 7 : for(i=0;i<6;i++)
343 : {
344 6 : int nVal = atoi(ExtractField(szField, szLine, 20 + 10*i, 10));
345 6 : if (nVal >= 2000000000)
346 0 : nVal = 0;
347 6 : ((int*)pabyImage)[i * nCells + nCellY * nRasterXSize + nCellX] = nVal;
348 : }
349 :
350 1 : nLine ++;
351 : }
352 :
353 1 : return TRUE;
354 : }
355 :
356 : /************************************************************************/
357 : /* Identify() */
358 : /************************************************************************/
359 :
360 10500 : int CTGDataset::Identify( GDALOpenInfo * poOpenInfo )
361 : {
362 10500 : CPLString osFilename(poOpenInfo->pszFilename);
363 :
364 10500 : 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 10500 : const char* pszFilename = CPLGetFilename(poOpenInfo->pszFilename);
368 10500 : 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 10500 : if (poOpenInfo->nHeaderBytes < HEADER_LINE_COUNT * 80)
381 : {
382 10360 : delete poOpenInfoToDelete;
383 10360 : return FALSE;
384 : }
385 :
386 : /* -------------------------------------------------------------------- */
387 : /* Chech that it looks roughly as a CTG dataset */
388 : /* -------------------------------------------------------------------- */
389 140 : const char* pszData = (const char*)poOpenInfo->pabyHeader;
390 : int i;
391 780 : for(i=0;i<4 * 80;i++)
392 : {
393 1376 : if (!((pszData[i] >= '0' && pszData[i] <= '9') ||
394 598 : pszData[i] == ' ' || pszData[i] == '-'))
395 : {
396 138 : delete poOpenInfoToDelete;
397 138 : return FALSE;
398 : }
399 : }
400 :
401 : char szField[11];
402 2 : int nRows = atoi(ExtractField(szField, pszData, 0, 10));
403 2 : int nCols = atoi(ExtractField(szField, pszData, 20, 10));
404 2 : int nMinColIndex = atoi(ExtractField(szField, pszData+80, 0, 5));
405 2 : int nMinRowIndex = atoi(ExtractField(szField, pszData+80, 5, 5));
406 2 : int nMaxColIndex = atoi(ExtractField(szField, pszData+80, 10, 5));
407 2 : int nMaxRowIndex = atoi(ExtractField(szField, pszData+80, 15, 5));
408 :
409 2 : 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 2 : delete poOpenInfoToDelete;
418 2 : return TRUE;
419 : }
420 :
421 :
422 : /************************************************************************/
423 : /* Open() */
424 : /************************************************************************/
425 :
426 1262 : GDALDataset *CTGDataset::Open( GDALOpenInfo * poOpenInfo )
427 :
428 : {
429 : int i;
430 :
431 1262 : if (!Identify(poOpenInfo))
432 1260 : return NULL;
433 :
434 2 : 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 2 : const char* pszFilename = CPLGetFilename(poOpenInfo->pszFilename);
439 2 : 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 2 : 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 2 : VSILFILE* fp = VSIFOpenL(osFilename.c_str(), "rb");
460 2 : if (fp == NULL)
461 0 : return NULL;
462 :
463 : char szHeader[HEADER_LINE_COUNT * 80+1];
464 2 : szHeader[HEADER_LINE_COUNT * 80] = 0;
465 2 : 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 288 : for(i=HEADER_LINE_COUNT * 80 - 1;i>=0;i--)
472 : {
473 144 : if (szHeader[i] == ' ')
474 142 : szHeader[i] = 0;
475 : else
476 2 : break;
477 : }
478 :
479 : char szField[11];
480 2 : int nRows = atoi(ExtractField(szField, szHeader, 0, 10));
481 2 : int nCols = atoi(ExtractField(szField, szHeader, 20, 10));
482 :
483 : /* -------------------------------------------------------------------- */
484 : /* Create a corresponding GDALDataset. */
485 : /* -------------------------------------------------------------------- */
486 : CTGDataset *poDS;
487 :
488 2 : poDS = new CTGDataset();
489 2 : poDS->fp = fp;
490 2 : fp = NULL;
491 2 : poDS->nRasterXSize = nCols;
492 2 : poDS->nRasterYSize = nRows;
493 :
494 2 : poDS->SetMetadataItem("TITLE", szHeader + 4 * 80);
495 :
496 2 : poDS->nCellSize = atoi(ExtractField(szField, szHeader, 35, 5));
497 2 : if (poDS->nCellSize <= 0 || poDS->nCellSize >= 10000)
498 : {
499 0 : delete poDS;
500 0 : return NULL;
501 : }
502 2 : poDS->nNWEasting = atoi(ExtractField(szField, szHeader + 3*80, 40, 10));
503 2 : poDS->nNWNorthing = atoi(ExtractField(szField, szHeader + 3*80, 50, 10));
504 2 : poDS->nUTMZone = atoi(ExtractField(szField, szHeader, 50, 5));
505 2 : if (poDS->nUTMZone <= 0 || poDS->nUTMZone > 60)
506 : {
507 0 : delete poDS;
508 0 : return NULL;
509 : }
510 :
511 2 : OGRSpatialReference oSRS;
512 2 : oSRS.importFromEPSG(32600 + poDS->nUTMZone);
513 2 : oSRS.exportToWkt(&poDS->pszProjection);
514 :
515 2 : if (!GDALCheckDatasetDimensions(poDS->nRasterXSize, poDS->nRasterYSize))
516 : {
517 0 : delete poDS;
518 0 : return NULL;
519 : }
520 :
521 : /* -------------------------------------------------------------------- */
522 : /* Read the imagery */
523 : /* -------------------------------------------------------------------- */
524 2 : GByte* pabyImage = (GByte*)VSICalloc(nCols * nRows, 6 * sizeof(int));
525 2 : if (pabyImage == NULL)
526 : {
527 0 : delete poDS;
528 0 : return NULL;
529 : }
530 2 : poDS->pabyImage = pabyImage;
531 :
532 : /* -------------------------------------------------------------------- */
533 : /* Create band information objects. */
534 : /* -------------------------------------------------------------------- */
535 2 : poDS->nBands = 6;
536 14 : for( i = 0; i < poDS->nBands; i++ )
537 : {
538 12 : poDS->SetBand( i+1, new CTGRasterBand( poDS, i+1 ) );
539 12 : poDS->GetRasterBand(i+1)->SetDescription(apszBandDescription[i]);
540 : }
541 :
542 : /* -------------------------------------------------------------------- */
543 : /* Initialize any PAM information. */
544 : /* -------------------------------------------------------------------- */
545 2 : poDS->SetDescription( poOpenInfo->pszFilename );
546 2 : poDS->TryLoadXML();
547 :
548 : /* -------------------------------------------------------------------- */
549 : /* Support overviews. */
550 : /* -------------------------------------------------------------------- */
551 2 : poDS->oOvManager.Initialize( poDS, poOpenInfo->pszFilename );
552 :
553 2 : return( poDS );
554 : }
555 :
556 : /************************************************************************/
557 : /* GetGeoTransform() */
558 : /************************************************************************/
559 :
560 1 : CPLErr CTGDataset::GetGeoTransform( double * padfTransform )
561 :
562 : {
563 1 : padfTransform[0] = nNWEasting - nCellSize / 2;
564 1 : padfTransform[1] = nCellSize;
565 1 : padfTransform[2] = 0;
566 1 : padfTransform[3] = nNWNorthing + nCellSize / 2;
567 1 : padfTransform[4] = 0.;
568 1 : padfTransform[5] = -nCellSize;
569 :
570 1 : return( CE_None );
571 : }
572 :
573 : /************************************************************************/
574 : /* GetProjectionRef() */
575 : /************************************************************************/
576 :
577 1 : const char* CTGDataset::GetProjectionRef()
578 :
579 : {
580 1 : return pszProjection;
581 : }
582 :
583 : /************************************************************************/
584 : /* GDALRegister_CTG() */
585 : /************************************************************************/
586 :
587 558 : void GDALRegister_CTG()
588 :
589 : {
590 : GDALDriver *poDriver;
591 :
592 558 : if( GDALGetDriverByName( "CTG" ) == NULL )
593 : {
594 537 : poDriver = new GDALDriver();
595 :
596 537 : poDriver->SetDescription( "CTG" );
597 : poDriver->SetMetadataItem( GDAL_DMD_LONGNAME,
598 537 : "USGS LULC Composite Theme Grid" );
599 : poDriver->SetMetadataItem( GDAL_DMD_HELPTOPIC,
600 537 : "frmt_various.html#CTG" );
601 :
602 537 : poDriver->SetMetadataItem( GDAL_DCAP_VIRTUALIO, "YES" );
603 :
604 537 : poDriver->pfnOpen = CTGDataset::Open;
605 537 : poDriver->pfnIdentify = CTGDataset::Identify;
606 :
607 537 : GetGDALDriverManager()->RegisterDriver( poDriver );
608 : }
609 558 : }
610 :
|
