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