|
osgEarth 2.1.1
|
|
osgEarth 2.1.1
|
#include <osgEarth/TileSource>#include <osgEarth/FileUtils>#include <osgEarth/Registry>#include <osgEarth/ImageUtils>#include <osgDB/FileNameUtils>#include <osgDB/FileUtils>#include <osgDB/Registry>#include <osgDB/ReadFile>#include <osgDB/WriteFile>#include <osgDB/ImageOptions>#include <sstream>#include <stdlib.h>#include <memory.h>#include <gdal_priv.h>#include <gdalwarper.h>#include <ogr_spatialref.h>#include "GDALOptions"#include <cpl_string.h>#include <gdal_vrt.h>
Include dependency graph for ReaderWriterGDAL.cpp:Go to the source code of this file.
Classes | |
| struct | DatasetProperty |
| struct | BandProperty |
| class | GDALTileSource |
| class | ReaderWriterGDALTile |
Defines | |
| #define | LC "[GDAL driver] " |
| #define | GEOTRSFRM_TOPLEFT_X 0 |
| #define | GEOTRSFRM_WE_RES 1 |
| #define | GEOTRSFRM_ROTATION_PARAM1 2 |
| #define | GEOTRSFRM_TOPLEFT_Y 3 |
| #define | GEOTRSFRM_ROTATION_PARAM2 4 |
| #define | GEOTRSFRM_NS_RES 5 |
Enumerations | |
| enum | ResolutionStrategy { LOWEST_RESOLUTION, HIGHEST_RESOLUTION, AVERAGE_RESOLUTION } |
Functions | |
| float | Hue_2_RGB (float v1, float v2, float vH) |
| static void | getFiles (const std::string &file, const std::vector< std::string > &exts, std::vector< std::string > &files) |
| static GDALDatasetH | build_vrt (std::vector< std::string > &files, ResolutionStrategy resolutionStrategy) |
| static GDALDatasetH | GDALAutoCreateWarpedVRTforPolarStereographic (GDALDatasetH hSrcDS, const char *pszSrcWKT, const char *pszDstWKT, GDALResampleAlg eResampleAlg, double dfMaxError, const GDALWarpOptions *psOptionsIn) |
| #define GEOTRSFRM_NS_RES 5 |
Definition at line 86 of file ReaderWriterGDAL.cpp.
| #define GEOTRSFRM_ROTATION_PARAM1 2 |
Definition at line 83 of file ReaderWriterGDAL.cpp.
| #define GEOTRSFRM_ROTATION_PARAM2 4 |
Definition at line 85 of file ReaderWriterGDAL.cpp.
| #define GEOTRSFRM_TOPLEFT_X 0 |
Definition at line 81 of file ReaderWriterGDAL.cpp.
| #define GEOTRSFRM_TOPLEFT_Y 3 |
Definition at line 84 of file ReaderWriterGDAL.cpp.
| #define GEOTRSFRM_WE_RES 1 |
Definition at line 82 of file ReaderWriterGDAL.cpp.
| #define LC "[GDAL driver] " |
Definition at line 42 of file ReaderWriterGDAL.cpp.
| enum ResolutionStrategy |
Definition at line 90 of file ReaderWriterGDAL.cpp.
| static GDALDatasetH build_vrt | ( | std::vector< std::string > & | files, |
| ResolutionStrategy | resolutionStrategy | ||
| ) | [static] |
Definition at line 184 of file ReaderWriterGDAL.cpp.
{
GDAL_SCOPED_LOCK;
char* projectionRef = NULL;
int nBands = 0;
BandProperty* bandProperties = NULL;
double minX = 0, minY = 0, maxX = 0, maxY = 0;
int i,j;
double we_res = 0;
double ns_res = 0;
int rasterXSize;
int rasterYSize;
int nCount = 0;
int bFirst = TRUE;
VRTDatasetH hVRTDS = NULL;
int nInputFiles = files.size();
DatasetProperty* psDatasetProperties =
(DatasetProperty*) CPLMalloc(nInputFiles*sizeof(DatasetProperty));
for(i=0;i<nInputFiles;i++)
{
const char* dsFileName = files[i].c_str();
GDALTermProgress( 1.0 * (i+1) / nInputFiles, NULL, NULL);
GDALDatasetH hDS = GDALOpen(dsFileName, GA_ReadOnly );
psDatasetProperties[i].isFileOK = FALSE;
if (hDS)
{
const char* proj = GDALGetProjectionRef(hDS);
if (!proj || strlen(proj) == 0)
{
std::string prjLocation = osgDB::getNameLessExtension( std::string(dsFileName) ) + std::string(".prj");
std::string wkt;
if ( HTTPClient::readString( prjLocation, wkt ) == HTTPClient::RESULT_OK )
{
proj = CPLStrdup( wkt.c_str() );
}
}
GDALGetGeoTransform(hDS, psDatasetProperties[i].adfGeoTransform);
if (psDatasetProperties[i].adfGeoTransform[GEOTRSFRM_ROTATION_PARAM1] != 0 ||
psDatasetProperties[i].adfGeoTransform[GEOTRSFRM_ROTATION_PARAM2] != 0)
{
fprintf( stderr, "GDAL Driver does not support rotated geo transforms. Skipping %s\n",
dsFileName);
GDALClose(hDS);
continue;
}
if (psDatasetProperties[i].adfGeoTransform[GEOTRSFRM_NS_RES] >= 0)
{
fprintf( stderr, "GDAL Driver does not support positive NS resolution. Skipping %s\n",
dsFileName);
GDALClose(hDS);
continue;
}
psDatasetProperties[i].nRasterXSize = GDALGetRasterXSize(hDS);
psDatasetProperties[i].nRasterYSize = GDALGetRasterYSize(hDS);
double product_minX = psDatasetProperties[i].adfGeoTransform[GEOTRSFRM_TOPLEFT_X];
double product_maxY = psDatasetProperties[i].adfGeoTransform[GEOTRSFRM_TOPLEFT_Y];
double product_maxX = product_minX +
GDALGetRasterXSize(hDS) * psDatasetProperties[i].adfGeoTransform[GEOTRSFRM_WE_RES];
double product_minY = product_maxY +
GDALGetRasterYSize(hDS) * psDatasetProperties[i].adfGeoTransform[GEOTRSFRM_NS_RES];
GDALGetBlockSize(GDALGetRasterBand( hDS, 1 ),
&psDatasetProperties[i].nBlockXSize,
&psDatasetProperties[i].nBlockYSize);
if (bFirst)
{
if (proj)
projectionRef = CPLStrdup(proj);
minX = product_minX;
minY = product_minY;
maxX = product_maxX;
maxY = product_maxY;
nBands = GDALGetRasterCount(hDS);
bandProperties = (BandProperty*)CPLMalloc(nBands*sizeof(BandProperty));
for(j=0;j<nBands;j++)
{
GDALRasterBandH hRasterBand = GDALGetRasterBand( hDS, j+1 );
bandProperties[j].colorInterpretation = GDALGetRasterColorInterpretation(hRasterBand);
bandProperties[j].dataType = GDALGetRasterDataType(hRasterBand);
if (bandProperties[j].colorInterpretation == GCI_PaletteIndex)
{
bandProperties[j].colorTable = GDALGetRasterColorTable( hRasterBand );
if (bandProperties[j].colorTable)
{
bandProperties[j].colorTable = GDALCloneColorTable(bandProperties[j].colorTable);
}
}
else
bandProperties[j].colorTable = 0;
bandProperties[j].noDataValue = GDALGetRasterNoDataValue(hRasterBand, &bandProperties[j].bHasNoData);
}
}
else
{
if ((proj != NULL && projectionRef == NULL) ||
(proj == NULL && projectionRef != NULL) ||
(proj != NULL && projectionRef != NULL && EQUAL(proj, projectionRef) == FALSE))
{
fprintf( stderr, "gdalbuildvrt does not support heterogenous projection. Skipping %s\n",dsFileName);
GDALClose(hDS);
continue;
}
int _nBands = GDALGetRasterCount(hDS);
if (nBands != _nBands)
{
fprintf( stderr, "gdalbuildvrt does not support heterogenous band numbers. Skipping %s\n",
dsFileName);
GDALClose(hDS);
continue;
}
for(j=0;j<nBands;j++)
{
GDALRasterBandH hRasterBand = GDALGetRasterBand( hDS, j+1 );
if (bandProperties[j].colorInterpretation != GDALGetRasterColorInterpretation(hRasterBand) ||
bandProperties[j].dataType != GDALGetRasterDataType(hRasterBand))
{
fprintf( stderr, "gdalbuildvrt does not support heterogenous band characteristics. Skipping %s\n",
dsFileName);
GDALClose(hDS);
}
if (bandProperties[j].colorTable)
{
GDALColorTableH colorTable = GDALGetRasterColorTable( hRasterBand );
if (colorTable == NULL ||
GDALGetColorEntryCount(colorTable) != GDALGetColorEntryCount(bandProperties[j].colorTable))
{
fprintf( stderr, "gdalbuildvrt does not support heterogenous band characteristics. Skipping %s\n",
dsFileName);
GDALClose(hDS);
break;
}
/* We should check that the palette are the same too ! */
}
}
if (j != nBands)
continue;
if (product_minX < minX) minX = product_minX;
if (product_minY < minY) minY = product_minY;
if (product_maxX > maxX) maxX = product_maxX;
if (product_maxY > maxY) maxY = product_maxY;
}
if (resolutionStrategy == AVERAGE_RESOLUTION)
{
we_res += psDatasetProperties[i].adfGeoTransform[GEOTRSFRM_WE_RES];
ns_res += psDatasetProperties[i].adfGeoTransform[GEOTRSFRM_NS_RES];
}
else
{
if (bFirst)
{
we_res = psDatasetProperties[i].adfGeoTransform[GEOTRSFRM_WE_RES];
ns_res = psDatasetProperties[i].adfGeoTransform[GEOTRSFRM_NS_RES];
}
else if (resolutionStrategy == HIGHEST_RESOLUTION)
{
we_res = MIN(we_res, psDatasetProperties[i].adfGeoTransform[GEOTRSFRM_WE_RES]);
/* Yes : as ns_res is negative, the highest resolution is the max value */
ns_res = MAX(ns_res, psDatasetProperties[i].adfGeoTransform[GEOTRSFRM_NS_RES]);
}
else
{
we_res = MAX(we_res, psDatasetProperties[i].adfGeoTransform[GEOTRSFRM_WE_RES]);
/* Yes : as ns_res is negative, the lowest resolution is the min value */
ns_res = MIN(ns_res, psDatasetProperties[i].adfGeoTransform[GEOTRSFRM_NS_RES]);
}
}
psDatasetProperties[i].isFileOK = 1;
nCount ++;
bFirst = FALSE;
GDALClose(hDS);
}
else
{
fprintf( stderr, "Warning : can't open %s. Skipping it\n", dsFileName);
}
}
if (nCount == 0)
goto end;
if (resolutionStrategy == AVERAGE_RESOLUTION)
{
we_res /= nCount;
ns_res /= nCount;
}
rasterXSize = (int)(0.5 + (maxX - minX) / we_res);
rasterYSize = (int)(0.5 + (maxY - minY) / -ns_res);
hVRTDS = VRTCreate(rasterXSize, rasterYSize);
if (projectionRef)
{
//OE_NOTICE << "Setting projection to " << projectionRef << std::endl;
GDALSetProjection(hVRTDS, projectionRef);
}
double adfGeoTransform[6];
adfGeoTransform[GEOTRSFRM_TOPLEFT_X] = minX;
adfGeoTransform[GEOTRSFRM_WE_RES] = we_res;
adfGeoTransform[GEOTRSFRM_ROTATION_PARAM1] = 0;
adfGeoTransform[GEOTRSFRM_TOPLEFT_Y] = maxY;
adfGeoTransform[GEOTRSFRM_ROTATION_PARAM2] = 0;
adfGeoTransform[GEOTRSFRM_NS_RES] = ns_res;
GDALSetGeoTransform(hVRTDS, adfGeoTransform);
for(j=0;j<nBands;j++)
{
GDALRasterBandH hBand;
GDALAddBand(hVRTDS, bandProperties[j].dataType, NULL);
hBand = GDALGetRasterBand(hVRTDS, j+1);
GDALSetRasterColorInterpretation(hBand, bandProperties[j].colorInterpretation);
if (bandProperties[j].colorInterpretation == GCI_PaletteIndex)
{
GDALSetRasterColorTable(hBand, bandProperties[j].colorTable);
}
if (bandProperties[j].bHasNoData)
GDALSetRasterNoDataValue(hBand, bandProperties[j].noDataValue);
}
for(i=0;i<nInputFiles;i++)
{
if (psDatasetProperties[i].isFileOK == 0)
continue;
const char* dsFileName = files[i].c_str();
bool isProxy = true;
#if GDAL_VERSION_1_6_OR_NEWER
//Use a proxy dataset if possible. This helps with huge amount of files to keep the # of handles down
GDALProxyPoolDatasetH hDS =
GDALProxyPoolDatasetCreate(dsFileName,
psDatasetProperties[i].nRasterXSize,
psDatasetProperties[i].nRasterYSize,
GA_ReadOnly, TRUE, projectionRef,
psDatasetProperties[i].adfGeoTransform);
for(j=0;j<nBands;j++)
{
GDALProxyPoolDatasetAddSrcBandDescription(hDS,
bandProperties[j].dataType,
psDatasetProperties[i].nBlockXSize,
psDatasetProperties[i].nBlockYSize);
}
isProxy = true;
OE_DEBUG << LC << "Using GDALProxyPoolDatasetH" << std::endl;
#else // !GDAL_VERSION_1_6_OR_NEWER
OE_DEBUG << LC << "Using GDALDataset, no proxy support enabled" << std::endl;
//Just open the dataset
GDALDatasetH hDS = (GDALDatasetH)GDALOpen(dsFileName, GA_ReadOnly);
isProxy = false;
#endif
int xoffset = (int)
(0.5 + (psDatasetProperties[i].adfGeoTransform[GEOTRSFRM_TOPLEFT_X] - minX) / we_res);
int yoffset = (int)
(0.5 + (maxY - psDatasetProperties[i].adfGeoTransform[GEOTRSFRM_TOPLEFT_Y]) / -ns_res);
int dest_width = (int)
(0.5 + psDatasetProperties[i].nRasterXSize * psDatasetProperties[i].adfGeoTransform[GEOTRSFRM_WE_RES] / we_res);
int dest_height = (int)
(0.5 + psDatasetProperties[i].nRasterYSize * psDatasetProperties[i].adfGeoTransform[GEOTRSFRM_NS_RES] / ns_res);
for(j=0;j<nBands;j++)
{
VRTSourcedRasterBandH hVRTBand = (VRTSourcedRasterBandH)GDALGetRasterBand(hVRTDS, j + 1);
/* Place the raster band at the right position in the VRT */
VRTAddSimpleSource(hVRTBand, GDALGetRasterBand((GDALDatasetH)hDS, j + 1),
0, 0,
psDatasetProperties[i].nRasterXSize,
psDatasetProperties[i].nRasterYSize,
xoffset, yoffset,
dest_width, dest_height, "near",
VRT_NODATA_UNSET);
}
//Only dereference if it is a proxy dataset
if (isProxy)
{
GDALDereferenceDataset(hDS);
}
}
end:
CPLFree(psDatasetProperties);
for(j=0;j<nBands;j++)
{
GDALDestroyColorTable(bandProperties[j].colorTable);
}
CPLFree(bandProperties);
CPLFree(projectionRef);
return hVRTDS;
}
Here is the caller graph for this function:| static GDALDatasetH GDALAutoCreateWarpedVRTforPolarStereographic | ( | GDALDatasetH | hSrcDS, |
| const char * | pszSrcWKT, | ||
| const char * | pszDstWKT, | ||
| GDALResampleAlg | eResampleAlg, | ||
| double | dfMaxError, | ||
| const GDALWarpOptions * | psOptionsIn | ||
| ) | [static] |
Definition at line 495 of file ReaderWriterGDAL.cpp.
{
GDALWarpOptions *psWO;
int i;
VALIDATE_POINTER1( hSrcDS, "GDALAutoCreateWarpedVRTForPolarStereographic", NULL );
/* -------------------------------------------------------------------- */
/* Populate the warp options. */
/* -------------------------------------------------------------------- */
if( psOptionsIn != NULL )
psWO = GDALCloneWarpOptions( psOptionsIn );
else
psWO = GDALCreateWarpOptions();
psWO->eResampleAlg = eResampleAlg;
psWO->hSrcDS = hSrcDS;
psWO->nBandCount = GDALGetRasterCount( hSrcDS );
psWO->panSrcBands = (int *) CPLMalloc(sizeof(int) * psWO->nBandCount);
psWO->panDstBands = (int *) CPLMalloc(sizeof(int) * psWO->nBandCount);
for( i = 0; i < psWO->nBandCount; i++ )
{
psWO->panSrcBands[i] = i+1;
psWO->panDstBands[i] = i+1;
}
/* TODO: should fill in no data where available */
/* -------------------------------------------------------------------- */
/* Create the transformer. */
/* -------------------------------------------------------------------- */
psWO->pfnTransformer = GDALGenImgProjTransform;
psWO->pTransformerArg =
GDALCreateGenImgProjTransformer( psWO->hSrcDS, pszSrcWKT,
NULL, pszDstWKT,
TRUE, 1.0, 0 );
if( psWO->pTransformerArg == NULL )
{
GDALDestroyWarpOptions( psWO );
return NULL;
}
/* -------------------------------------------------------------------- */
/* Figure out the desired output bounds and resolution. */
/* -------------------------------------------------------------------- */
double adfDstGeoTransform[6];
int nDstPixels, nDstLines;
CPLErr eErr;
eErr =
GDALSuggestedWarpOutput( hSrcDS, psWO->pfnTransformer,
psWO->pTransformerArg,
adfDstGeoTransform, &nDstPixels, &nDstLines );
// override the suggestions:
nDstPixels = GDALGetRasterXSize( hSrcDS ) * 4;
nDstLines = GDALGetRasterYSize( hSrcDS ) / 2;
adfDstGeoTransform[0] = -180.0;
adfDstGeoTransform[1] = 360.0/(double)nDstPixels;
//adfDstGeoTransform[2] = 0.0;
//adfDstGeoTransform[4] = 0.0;
//adfDstGeoTransform[5] = (-90 -adfDstGeoTransform[3])/(double)nDstLines;
/* -------------------------------------------------------------------- */
/* Update the transformer to include an output geotransform */
/* back to pixel/line coordinates. */
/* */
/* -------------------------------------------------------------------- */
GDALSetGenImgProjTransformerDstGeoTransform(
psWO->pTransformerArg, adfDstGeoTransform );
/* -------------------------------------------------------------------- */
/* Do we want to apply an approximating transformation? */
/* -------------------------------------------------------------------- */
if( dfMaxError > 0.0 )
{
psWO->pTransformerArg =
GDALCreateApproxTransformer( psWO->pfnTransformer,
psWO->pTransformerArg,
dfMaxError );
psWO->pfnTransformer = GDALApproxTransform;
}
/* -------------------------------------------------------------------- */
/* Create the VRT file. */
/* -------------------------------------------------------------------- */
GDALDatasetH hDstDS;
hDstDS = GDALCreateWarpedVRT( hSrcDS, nDstPixels, nDstLines,
adfDstGeoTransform, psWO );
GDALDestroyWarpOptions( psWO );
if( pszDstWKT != NULL )
GDALSetProjection( hDstDS, pszDstWKT );
else if( pszSrcWKT != NULL )
GDALSetProjection( hDstDS, pszDstWKT );
else if( GDALGetGCPCount( hSrcDS ) > 0 )
GDALSetProjection( hDstDS, GDALGetGCPProjection( hSrcDS ) );
else
GDALSetProjection( hDstDS, GDALGetProjectionRef( hSrcDS ) );
return hDstDS;
}
Here is the caller graph for this function:| static void getFiles | ( | const std::string & | file, |
| const std::vector< std::string > & | exts, | ||
| std::vector< std::string > & | files | ||
| ) | [static] |
Definition at line 117 of file ReaderWriterGDAL.cpp.
{
if (osgDB::fileType(file) == osgDB::DIRECTORY)
{
osgDB::DirectoryContents contents = osgDB::getDirectoryContents(file);
for (osgDB::DirectoryContents::iterator itr = contents.begin(); itr != contents.end(); ++itr)
{
if (*itr == "." || *itr == "..") continue;
std::string f = osgDB::concatPaths(file, *itr);
getFiles(f, exts, files);
}
}
else
{
bool fileValid = false;
//If we have no _extensions specified, assume we should try everything
if (exts.size() == 0)
{
fileValid = true;
}
else
{
//Only accept files with the given _extensions
std::string ext = osgDB::getFileExtension(file);
for (unsigned int i = 0; i < exts.size(); ++i)
{
if (osgDB::equalCaseInsensitive(ext, exts[i]))
{
fileValid = true;
break;
}
}
}
if (fileValid)
{
files.push_back(osgDB::convertFileNameToNativeStyle(file));
}
}
}
Here is the caller graph for this function:| float Hue_2_RGB | ( | float | v1, |
| float | v2, | ||
| float | vH | ||
| ) |
Definition at line 45 of file ReaderWriterGDAL.cpp.
{
if ( vH < 0 ) vH += 1;
if ( vH > 1 ) vH -= 1;
if ( ( 6 * vH ) < 1 ) return ( v1 + ( v2 - v1 ) * 6 * vH );
if ( ( 2 * vH ) < 1 ) return ( v2 );
if ( ( 3 * vH ) < 2 ) return ( v1 + ( v2 - v1 ) * ( ( 2 / 3 ) - vH ) * 6 );
return ( v1 );
}
Here is the caller graph for this function:
1.7.3
