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