/****************************************************************************** * $Id: grass57dataset.cpp 11554 2007-05-18 17:09:03Z mloskot $ * * Project: GRASS Driver * Purpose: Implement GRASS raster read/write support * This version is for GRASS 5.7+ and uses GRASS libraries * directly instead of using libgrass. * Author: Frank Warmerdam * Radim Blazek * ****************************************************************************** * Copyright (c) 2000 Frank Warmerdam * * Permission is hereby granted, free of charge, to any person obtaining a * copy of this software and associated documentation files (the "Software"), * to deal in the Software without restriction, including without limitation * the rights to use, copy, modify, merge, publish, distribute, sublicense, * and/or sell copies of the Software, and to permit persons to whom the * Software is furnished to do so, subject to the following conditions: * * The above copyright notice and this permission notice shall be included * in all copies or substantial portions of the Software. * * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER * DEALINGS IN THE SOFTWARE. ****************************************************************************/ #include extern "C" { #ifdef __cplusplus #define class _class #endif #include #ifdef __cplusplus #undef class #endif #include #include #include } #include "gdal_priv.h" #include "cpl_string.h" #include "ogr_spatialref.h" #define GRASS_MAX_COLORS 100000 // what is the right value CPL_CVSID("$Id: grass57dataset.cpp 11554 2007-05-18 17:09:03Z mloskot $"); CPL_C_START void GDALRegister_GRASS(void); CPL_C_END /************************************************************************/ /* Grass2CPLErrorHook() */ /************************************************************************/ int Grass2CPLErrorHook( char * pszMessage, int bFatal ) { if( !bFatal ) //CPLDebug( "GRASS", "%s", pszMessage ); CPLError( CE_Warning, CPLE_AppDefined, "GRASS warning: %s", pszMessage ); else CPLError( CE_Warning, CPLE_AppDefined, "GRASS fatal error: %s", pszMessage ); return 0; } /************************************************************************/ /* ==================================================================== */ /* GRASSDataset */ /* ==================================================================== */ /************************************************************************/ class GRASSRasterBand; class GRASSDataset : public GDALDataset { friend class GRASSRasterBand; char *pszGisdbase; char *pszLocation; /* LOCATION_NAME */ char *pszElement; /* cellhd or group */ struct Cell_head sCellInfo; /* raster region */ char *pszProjection; double adfGeoTransform[6]; public: GRASSDataset(); ~GRASSDataset(); virtual const char *GetProjectionRef(void); virtual CPLErr GetGeoTransform( double * ); static GDALDataset *Open( GDALOpenInfo * ); private: static bool SplitPath ( char *, char **, char **, char **, char **, char ** ); }; /************************************************************************/ /* ==================================================================== */ /* GRASSRasterBand */ /* ==================================================================== */ /************************************************************************/ class GRASSRasterBand : public GDALRasterBand { friend class GRASSDataset; char *pszCellName; char *pszMapset; int hCell; int nGRSType; // GRASS raster type: CELL_TYPE, FCELL_TYPE, DCELL_TYPE bool nativeNulls; // use GRASS native NULL values struct Colors sGrassColors; GDALColorTable *poCT; struct Cell_head sOpenWindow; /* the region when the raster was opened */ int bHaveMinMax; double dfCellMin; double dfCellMax; double dfNoData; bool valid; public: GRASSRasterBand( GRASSDataset *, int, const char *, const char * ); virtual ~GRASSRasterBand(); virtual CPLErr IReadBlock( int, int, void * ); virtual CPLErr IRasterIO ( GDALRWFlag, int, int, int, int, void *, int, int, GDALDataType, int, int ); virtual GDALColorInterp GetColorInterpretation(); virtual GDALColorTable *GetColorTable(); virtual double GetMinimum( int *pbSuccess = NULL ); virtual double GetMaximum( int *pbSuccess = NULL ); virtual double GetNoDataValue( int *pbSuccess = NULL ); private: CPLErr ResetReading( struct Cell_head * ); }; /************************************************************************/ /* GRASSRasterBand() */ /************************************************************************/ GRASSRasterBand::GRASSRasterBand( GRASSDataset *poDS, int nBand, const char * pszMapset, const char * pszCellName ) { struct Cell_head sCellInfo; // Note: GISDBASE, LOCATION_NAME ans MAPSET was set in GRASSDataset::Open this->poDS = poDS; this->nBand = nBand; this->valid = false; this->pszCellName = G_store ( (char *) pszCellName ); this->pszMapset = G_store ( (char *) pszMapset ); G_get_cellhd( (char *) pszCellName, (char *) pszMapset, &sCellInfo ); nGRSType = G_raster_map_type( (char *) pszCellName, (char *) pszMapset ); /* -------------------------------------------------------------------- */ /* Get min/max values. */ /* -------------------------------------------------------------------- */ struct FPRange sRange; if( G_read_fp_range( (char *) pszCellName, (char *) pszMapset, &sRange ) == -1 ) { bHaveMinMax = FALSE; } else { bHaveMinMax = TRUE; G_get_fp_range_min_max( &sRange, &dfCellMin, &dfCellMax ); } /* -------------------------------------------------------------------- */ /* Setup band type, and preferred nodata value. */ /* -------------------------------------------------------------------- */ // Negative values are also (?) stored as 4 bytes (format = 3) // => raster with format < 3 has only positive values // GRASS modules usually do not waste space and only the format necessary to keep // full raster values range is used -> no checks if shorter type could be used if( nGRSType == CELL_TYPE ) { if ( sCellInfo.format == 0 ) { // 1 byte / cell -> possible range 0,255 if ( bHaveMinMax && dfCellMin > 0 ) { this->eDataType = GDT_Byte; dfNoData = 0.0; } else if ( bHaveMinMax && dfCellMax < 255 ) { this->eDataType = GDT_Byte; dfNoData = 255.0; } else { // maximum is not known or full range is used this->eDataType = GDT_UInt16; dfNoData = 256.0; } nativeNulls = false; } else if ( sCellInfo.format == 1 ) { // 2 bytes / cell -> possible range 0,65535 if ( bHaveMinMax && dfCellMin > 0 ) { this->eDataType = GDT_UInt16; dfNoData = 0.0; } else if ( bHaveMinMax && dfCellMax < 65535 ) { this->eDataType = GDT_UInt16; dfNoData = 65535; } else { // maximum is not known or full range is used CELL cval; this->eDataType = GDT_Int32; G_set_c_null_value ( &cval, 1); dfNoData = (double) cval; nativeNulls = true; } nativeNulls = false; } else { // 3-4 bytes CELL cval; this->eDataType = GDT_Int32; G_set_c_null_value ( &cval, 1); dfNoData = (double) cval; nativeNulls = true; } } else if( nGRSType == FCELL_TYPE ) { FCELL fval; this->eDataType = GDT_Float32; G_set_f_null_value ( &fval, 1); dfNoData = (double) fval; nativeNulls = true; } else if( nGRSType == DCELL_TYPE ) { DCELL dval; this->eDataType = GDT_Float64; G_set_d_null_value ( &dval, 1); dfNoData = (double) dval; nativeNulls = true; } nBlockXSize = poDS->nRasterXSize;; nBlockYSize = 1; G_set_window( &(((GRASSDataset *)poDS)->sCellInfo) ); if ( (hCell = G_open_cell_old((char *) pszCellName, (char *) pszMapset)) < 0 ) { CPLError( CE_Warning, CPLE_AppDefined, "GRASS: Cannot open raster '%s'", pszCellName ); return; } G_copy((void *) &sOpenWindow, (void *) &(((GRASSDataset *)poDS)->sCellInfo), sizeof(struct Cell_head)); /* -------------------------------------------------------------------- */ /* Do we have a color table? */ /* -------------------------------------------------------------------- */ poCT = NULL; if( G_read_colors( (char *) pszCellName, (char *) pszMapset, &sGrassColors ) == 1 ) { int maxcolor; CELL min, max; G_get_color_range ( &min, &max, &sGrassColors); if ( bHaveMinMax ) { if ( max < dfCellMax ) { maxcolor = max; } else { maxcolor = (int) ceil ( dfCellMax ); } if ( maxcolor > GRASS_MAX_COLORS ) { maxcolor = GRASS_MAX_COLORS; CPLDebug( "GRASS", "Too many values, color table cut to %d entries.", maxcolor ); } } else { if ( max < GRASS_MAX_COLORS ) { maxcolor = max; } else { maxcolor = GRASS_MAX_COLORS; CPLDebug( "GRASS", "Too many values, color table set to %d entries.", maxcolor ); } } poCT = new GDALColorTable(); for( int iColor = 0; iColor <= maxcolor; iColor++ ) { int nRed, nGreen, nBlue; GDALColorEntry sColor; if( G_get_color( iColor, &nRed, &nGreen, &nBlue, &sGrassColors ) ) { sColor.c1 = nRed; sColor.c2 = nGreen; sColor.c3 = nBlue; sColor.c4 = 255; poCT->SetColorEntry( iColor, &sColor ); } else { sColor.c1 = 0; sColor.c2 = 0; sColor.c3 = 0; sColor.c4 = 0; poCT->SetColorEntry( iColor, &sColor ); } } /* Create metadata enries for color table rules */ char key[200], value[200]; int rcount = G_colors_count ( &sGrassColors ); sprintf ( value, "%d", rcount ); this->SetMetadataItem( "COLOR_TABLE_RULES_COUNT", value ); /* Add the rules in reverse order */ for ( int i = rcount-1; i >= 0; i-- ) { DCELL val1, val2; unsigned char r1, g1, b1, r2, g2, b2; G_get_f_color_rule ( &val1, &r1, &g1, &b1, &val2, &r2, &g2, &b2, &sGrassColors, i ); sprintf ( key, "COLOR_TABLE_RULE_RGB_%d", rcount-i-1 ); sprintf ( value, "%e %e %d %d %d %d %d %d", val1, val2, r1, g1, b1, r2, g2, b2 ); this->SetMetadataItem( key, value ); } } else { this->SetMetadataItem( "COLOR_TABLE_RULES_COUNT", "0" ); } this->valid = true; } /************************************************************************/ /* ~GRASSRasterBand() */ /************************************************************************/ GRASSRasterBand::~GRASSRasterBand() { if( poCT != NULL ) { G_free_colors( &sGrassColors ); delete poCT; } if( hCell >= 0 ) G_close_cell( hCell ); if ( pszCellName ) free ( pszCellName ); if ( pszMapset ) free ( pszMapset ); } /************************************************************************/ /* ResetReading */ /* */ /* Reset current window and reopen cell if the window has changed, */ /* reset GRASS variables */ /* */ /* Returns CE_Failure if fails, otherwise CE_None */ /************************************************************************/ CPLErr GRASSRasterBand::ResetReading ( struct Cell_head *sNewWindow ) { /* Check if the window has changed */ if ( sNewWindow->north != sOpenWindow.north || sNewWindow->south != sOpenWindow.south || sNewWindow->east != sOpenWindow.east || sNewWindow->west != sOpenWindow.west || sNewWindow->ew_res != sOpenWindow.ew_res || sNewWindow->ns_res != sOpenWindow.ns_res || sNewWindow->rows != sOpenWindow.rows || sNewWindow->cols != sOpenWindow.cols ) { if( hCell >= 0 ) { G_close_cell( hCell ); hCell = -1; } /* Set window */ G_set_window( sNewWindow ); /* Open raster */ G__setenv( "GISDBASE", ((GRASSDataset *)poDS)->pszGisdbase ); G__setenv( "LOCATION_NAME", ((GRASSDataset *)poDS)->pszLocation ); G__setenv( "MAPSET", pszMapset); G_reset_mapsets(); G_add_mapset_to_search_path ( pszMapset ); if ( (hCell = G_open_cell_old( pszCellName, pszMapset)) < 0 ) { CPLError( CE_Warning, CPLE_AppDefined, "GRASS: Cannot open raster '%s'", pszCellName ); this->valid = false; return CE_Failure; } G_copy((void *) &sOpenWindow, (void *) sNewWindow, sizeof(struct Cell_head)); } else { /* The windows are identical, check current window */ struct Cell_head sCurrentWindow; G_get_window ( &sCurrentWindow ); if ( sNewWindow->north != sCurrentWindow.north || sNewWindow->south != sCurrentWindow.south || sNewWindow->east != sCurrentWindow.east || sNewWindow->west != sCurrentWindow.west || sNewWindow->ew_res != sCurrentWindow.ew_res || sNewWindow->ns_res != sCurrentWindow.ns_res || sNewWindow->rows != sCurrentWindow.rows || sNewWindow->cols != sCurrentWindow.cols ) { /* Reset window */ G_set_window( sNewWindow ); } } return CE_None; } /************************************************************************/ /* IReadBlock() */ /* */ /************************************************************************/ CPLErr GRASSRasterBand::IReadBlock( int nBlockXOff, int nBlockYOff, void * pImage ) { if ( ! this->valid ) return CE_Failure; // Reset window because IRasterIO could be previosly called if ( ResetReading ( &(((GRASSDataset *)poDS)->sCellInfo) ) != CE_None ) { return CE_Failure; } if ( eDataType == GDT_Byte || eDataType == GDT_UInt16 ) { CELL *cbuf; cbuf = G_allocate_c_raster_buf(); G_get_c_raster_row ( hCell, cbuf, nBlockYOff ); /* Reset NULLs */ for ( int col = 0; col < nBlockXSize; col++ ) { if ( G_is_c_null_value(&(cbuf[col])) ) cbuf[col] = (CELL) dfNoData; } GDALCopyWords ( (void *) cbuf, GDT_Int32, sizeof(CELL), pImage, eDataType, GDALGetDataTypeSize(eDataType)/8, nBlockXSize ); free ( cbuf ); } else if ( eDataType == GDT_Int32 ) { G_get_c_raster_row ( hCell, (CELL *) pImage, nBlockYOff ); } else if ( eDataType == GDT_Float32 ) { G_get_f_raster_row ( hCell, (FCELL *) pImage, nBlockYOff ); } else if ( eDataType == GDT_Float64 ) { G_get_d_raster_row ( hCell, (DCELL *) pImage, nBlockYOff ); } return CE_None; } /************************************************************************/ /* IRasterIO() */ /* */ /************************************************************************/ CPLErr GRASSRasterBand::IRasterIO ( GDALRWFlag eRWFlag, int nXOff, int nYOff, int nXSize, int nYSize, void * pData, int nBufXSize, int nBufYSize, GDALDataType eBufType, int nPixelSpace, int nLineSpace ) { /* GRASS library does that, we have only calculate and reset the region in map units * and if the region has changed, reopen the raster */ /* Calculate the region */ struct Cell_head sWindow; struct Cell_head *psDsWindow; if ( ! this->valid ) return CE_Failure; psDsWindow = &(((GRASSDataset *)poDS)->sCellInfo); sWindow.north = psDsWindow->north - nYOff * psDsWindow->ns_res; sWindow.south = sWindow.north - nYSize * psDsWindow->ns_res; sWindow.west = psDsWindow->west + nXOff * psDsWindow->ew_res; sWindow.east = sWindow.west + nXSize * psDsWindow->ew_res; sWindow.proj = psDsWindow->proj; sWindow.zone = psDsWindow->zone; sWindow.cols = nBufXSize; sWindow.rows = nBufYSize; /* Reset resolution */ G_adjust_Cell_head ( &sWindow, 1, 1); if ( ResetReading ( &sWindow ) != CE_None ) { return CE_Failure; } /* Read Data */ CELL *cbuf = NULL; FCELL *fbuf = NULL; DCELL *dbuf = NULL; bool direct = false; /* Reset space if default (0) */ if ( nPixelSpace == 0 ) nPixelSpace = GDALGetDataTypeSize ( eBufType ) / 8; if ( nLineSpace == 0 ) nLineSpace = nBufXSize * nPixelSpace; if ( nGRSType == CELL_TYPE && ( !nativeNulls || eBufType != GDT_Int32 || sizeof(CELL) != 4 || nPixelSpace != sizeof(CELL) ) ) { cbuf = G_allocate_c_raster_buf(); } else if( nGRSType == FCELL_TYPE && ( eBufType != GDT_Float32 || nPixelSpace != sizeof(FCELL) ) ) { fbuf = G_allocate_f_raster_buf(); } else if( nGRSType == DCELL_TYPE && ( eBufType != GDT_Float64 || nPixelSpace != sizeof(DCELL) ) ) { dbuf = G_allocate_d_raster_buf(); } else { direct = true; } for ( int row = 0; row < nBufYSize; row++ ) { char *pnt = (char *)pData + row * nLineSpace; if ( nGRSType == CELL_TYPE ) { if ( direct ) { G_get_c_raster_row ( hCell, (CELL *) pnt, row ); } else { G_get_c_raster_row ( hCell, cbuf, row ); /* Reset NULLs */ for ( int col = 0; col < nBufXSize; col++ ) { if ( G_is_c_null_value(&(cbuf[col])) ) cbuf[col] = (CELL) dfNoData; } GDALCopyWords ( (void *) cbuf, GDT_Int32, sizeof(CELL), (void *) pnt, eBufType, nPixelSpace, nBufXSize ); } } else if( nGRSType == FCELL_TYPE ) { if ( direct ) { G_get_f_raster_row ( hCell, (FCELL *) pnt, row ); } else { G_get_f_raster_row ( hCell, fbuf, row ); GDALCopyWords ( (void *) fbuf, GDT_Float32, sizeof(FCELL), (void *) pnt, eBufType, nPixelSpace, nBufXSize ); } } else if( nGRSType == DCELL_TYPE ) { if ( direct ) { G_get_d_raster_row ( hCell, (DCELL *) pnt, row ); } else { G_get_d_raster_row ( hCell, dbuf, row ); GDALCopyWords ( (void *) dbuf, GDT_Float64, sizeof(DCELL), (void *) pnt, eBufType, nPixelSpace, nBufXSize ); } } } if ( cbuf ) free ( cbuf ); if ( fbuf ) free ( fbuf ); if ( dbuf ) free ( dbuf ); return CE_None; } /************************************************************************/ /* GetColorInterpretation() */ /************************************************************************/ GDALColorInterp GRASSRasterBand::GetColorInterpretation() { if( poCT != NULL ) return GCI_PaletteIndex; else return GCI_GrayIndex; } /************************************************************************/ /* GetColorTable() */ /************************************************************************/ GDALColorTable *GRASSRasterBand::GetColorTable() { return poCT; } /************************************************************************/ /* GetMinimum() */ /************************************************************************/ double GRASSRasterBand::GetMinimum( int *pbSuccess ) { if( pbSuccess ) *pbSuccess = bHaveMinMax; if( bHaveMinMax ) return dfCellMin; else if( eDataType == GDT_Float32 || eDataType == GDT_Float64 ) return -4294967295.0; else return 0; } /************************************************************************/ /* GetMaximum() */ /************************************************************************/ double GRASSRasterBand::GetMaximum( int *pbSuccess ) { if( pbSuccess ) *pbSuccess = bHaveMinMax; if( bHaveMinMax ) return dfCellMax; else if( eDataType == GDT_Float32 || eDataType == GDT_Float64 ) return 4294967295.0; else if( eDataType == GDT_UInt32 ) return 4294967295.0; else if( eDataType == GDT_UInt16 ) return 65535; else return 255; } /************************************************************************/ /* GetNoDataValue() */ /************************************************************************/ double GRASSRasterBand::GetNoDataValue( int *pbSuccess ) { if( pbSuccess ) *pbSuccess = TRUE; return dfNoData; } /************************************************************************/ /* ==================================================================== */ /* GRASSDataset */ /* ==================================================================== */ /************************************************************************/ /************************************************************************/ /* GRASSDataset() */ /************************************************************************/ GRASSDataset::GRASSDataset() { pszProjection = NULL; adfGeoTransform[0] = 0.0; adfGeoTransform[1] = 1.0; adfGeoTransform[2] = 0.0; adfGeoTransform[3] = 0.0; adfGeoTransform[4] = 0.0; adfGeoTransform[5] = 1.0; } /************************************************************************/ /* ~GRASSDataset() */ /************************************************************************/ GRASSDataset::~GRASSDataset() { if ( pszGisdbase ) free ( pszGisdbase ); if ( pszLocation ) free ( pszLocation ); if ( pszElement ) free ( pszElement ); CPLFree( pszProjection ); } /************************************************************************/ /* GetProjectionRef() */ /************************************************************************/ const char *GRASSDataset::GetProjectionRef() { if( pszProjection == NULL ) return ""; else return pszProjection; } /************************************************************************/ /* GetGeoTransform() */ /************************************************************************/ CPLErr GRASSDataset::GetGeoTransform( double * padfGeoTransform ) { memcpy( padfGeoTransform, adfGeoTransform, sizeof(double) * 6 ); return CE_None; } /************************************************************************/ /* SplitPath() */ /* Split full path to cell or group to: */ /* gisdbase, location, mapset, element, name */ /* New string are allocated and should be freed when no longer needed. */ /* */ /* Returns: true - OK */ /* false - failed */ /************************************************************************/ bool GRASSDataset::SplitPath( char *path, char **gisdbase, char **location, char **mapset, char **element, char **name ) { char *p, *ptr[5], *tmp; int i = 0; *gisdbase = *location = *mapset = *element = *name = NULL; if ( !path || strlen(path) == 0 ) return false; tmp = G_store ( path ); while ( (p = strrchr(tmp,'/')) != NULL && i < 4 ) { *p = '\0'; if ( strlen(p+1) == 0 ) /* repeated '/' */ continue; ptr[i++] = p+1; } /* Note: empty GISDBASE == 0 is not accepted (relative path) */ if ( i != 4 ) { free ( tmp ); return false; } *gisdbase = G_store ( tmp ); *location = G_store ( ptr[3] ); *mapset = G_store ( ptr[2] ); *element = G_store ( ptr[1] ); *name = G_store ( ptr[0] ); free ( tmp ); return true; } /************************************************************************/ /* Open() */ /************************************************************************/ #if GRASS_VERSION_MAJOR >= 6 && GRASS_VERSION_MINOR >= 3 typedef int (*GrassErrorHandler)(const char *, int); #else typedef int (*GrassErrorHandler)(char *, int); #endif GDALDataset *GRASSDataset::Open( GDALOpenInfo * poOpenInfo ) { char *pszGisdb = NULL, *pszLoc = NULL; char *pszMapset = NULL, *pszElem = NULL, *pszName = NULL; char **papszCells = NULL; char **papszMapsets = NULL; static char fake_gisbase[50]; static bool hasGisbase; /* -------------------------------------------------------------------- */ /* Does this even look like a grass file path? */ /* -------------------------------------------------------------------- */ if( strstr(poOpenInfo->pszFilename,"/cellhd/") == NULL && strstr(poOpenInfo->pszFilename,"/group/") == NULL ) return NULL; /* Always init, if no rasters are opened G_no_gisinit resets the projection and * rasters in different projection may be then opened */ // Don't use GISRC file and read/write GRASS variables (from location G_VAR_GISRC) to memory only. G_set_gisrc_mode ( G_GISRC_MODE_MEMORY ); // Init GRASS libraries (required) G_no_gisinit(); // Doesn't check write permissions for mapset compare to G_gisinit // Set error function G_set_error_routine ( (GrassErrorHandler) Grass2CPLErrorHook ); if ( !getenv( "GISBASE" ) ) { // we are outside a GRASS session // TODO: use function instead of hardcoded path sprintf(fake_gisbase, "GISBASE=/usr/local/share/gdal/grass/" ); putenv( fake_gisbase ); hasGisbase = false; } else { hasGisbase = true; } if ( !SplitPath( poOpenInfo->pszFilename, &pszGisdb, &pszLoc, &pszMapset, &pszElem, &pszName) ) { return NULL; } /* -------------------------------------------------------------------- */ /* Check element name */ /* -------------------------------------------------------------------- */ if ( strcmp(pszElem,"cellhd") != 0 && strcmp(pszElem,"group") != 0 ) { free(pszGisdb); free(pszLoc); free(pszMapset); free(pszElem); free(pszName); return NULL; } /* -------------------------------------------------------------------- */ /* Set GRASS variables */ /* -------------------------------------------------------------------- */ G__setenv( "GISDBASE", pszGisdb ); G__setenv( "LOCATION_NAME", pszLoc ); G__setenv( "MAPSET", pszMapset); // group is searched only in current mapset G_reset_mapsets(); G_add_mapset_to_search_path ( pszMapset ); /* -------------------------------------------------------------------- */ /* Check if this is a valid grass cell. */ /* -------------------------------------------------------------------- */ if ( strcmp(pszElem,"cellhd") == 0 ) { if ( G_find_file2("cell", pszName, pszMapset) == NULL ) { free(pszGisdb); free(pszLoc); free(pszMapset); free(pszElem); free(pszName); return NULL; } papszMapsets = CSLAddString( papszMapsets, pszMapset ); papszCells = CSLAddString( papszCells, pszName ); } /* -------------------------------------------------------------------- */ /* Check if this is a valid GRASS imagery group. */ /* -------------------------------------------------------------------- */ else { struct Ref ref; I_init_group_ref( &ref ); if ( I_get_group_ref( pszName, &ref ) == 0 ) { free(pszGisdb); free(pszLoc); free(pszMapset); free(pszElem); free(pszName); return NULL; } for( int iRef = 0; iRef < ref.nfiles; iRef++ ) { papszCells = CSLAddString( papszCells, ref.file[iRef].name ); papszMapsets = CSLAddString( papszMapsets, ref.file[iRef].mapset ); G_add_mapset_to_search_path ( ref.file[iRef].mapset ); } I_free_group_ref( &ref ); } free( pszMapset ); free( pszName ); /* -------------------------------------------------------------------- */ /* Create a corresponding GDALDataset. */ /* -------------------------------------------------------------------- */ GRASSDataset *poDS; poDS = new GRASSDataset(); /* notdef: should only allow read access to an existing cell, right? */ poDS->eAccess = poOpenInfo->eAccess; poDS->pszGisdbase = pszGisdb; poDS->pszLocation = pszLoc; poDS->pszElement = pszElem; /* -------------------------------------------------------------------- */ /* Capture some information from the file that is of interest. */ /* -------------------------------------------------------------------- */ if( G_get_cellhd( papszCells[0], papszMapsets[0], &(poDS->sCellInfo) ) != 0 ) { CPLError( CE_Warning, CPLE_AppDefined, "GRASS: Cannot open raster header"); /* TODO: delete poDS ? */ return NULL; } poDS->nRasterXSize = poDS->sCellInfo.cols; poDS->nRasterYSize = poDS->sCellInfo.rows; poDS->adfGeoTransform[0] = poDS->sCellInfo.west; poDS->adfGeoTransform[1] = poDS->sCellInfo.ew_res; poDS->adfGeoTransform[2] = 0.0; poDS->adfGeoTransform[3] = poDS->sCellInfo.north; poDS->adfGeoTransform[4] = 0.0; poDS->adfGeoTransform[5] = -1 * poDS->sCellInfo.ns_res; /* -------------------------------------------------------------------- */ /* Try to get a projection definition. */ /* -------------------------------------------------------------------- */ struct Key_Value *projinfo, *projunits; if ( hasGisbase ) { projinfo = G_get_projinfo(); projunits = G_get_projunits(); poDS->pszProjection = GPJ_grass_to_wkt ( projinfo, projunits, 0, 0); } /* -------------------------------------------------------------------- */ /* Create band information objects. */ /* -------------------------------------------------------------------- */ for( int iBand = 0; papszCells[iBand] != NULL; iBand++ ) { GRASSRasterBand *rb = new GRASSRasterBand( poDS, iBand+1, papszMapsets[iBand], papszCells[iBand] ); if ( !rb->valid ) { CPLError( CE_Warning, CPLE_AppDefined, "GRASS: Cannot open raster band %d", iBand); // TODO: delete poDS ? return NULL; } poDS->SetBand( iBand+1, rb ); } return poDS; } /************************************************************************/ /* GDALRegister_GRASS() */ /************************************************************************/ void GDALRegister_GRASS() { GDALDriver *poDriver; if( GDALGetDriverByName( "GRASS" ) == NULL ) { poDriver = new GDALDriver(); poDriver->SetDescription( "GRASS" ); poDriver->SetMetadataItem( GDAL_DMD_LONGNAME, "GRASS Database Rasters (5.7+)" ); poDriver->SetMetadataItem( GDAL_DMD_HELPTOPIC, "frmt_grass.html" ); poDriver->pfnOpen = GRASSDataset::Open; GetGDALDriverManager()->RegisterDriver( poDriver ); } }