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