|
osgEarth 2.1.1
|
|
osgEarth 2.1.1
|
Inheritance diagram for GDALTileSource: Collaboration diagram for GDALTileSource:Public Member Functions | |
| GDALTileSource (const TileSourceOptions &options) | |
| virtual | ~GDALTileSource () |
| void | initialize (const std::string &referenceURI, const Profile *overrideProfile) |
| void | pixelToGeo (double x, double y, double &geoX, double &geoY) |
| osg::Image * | createImage (const TileKey &key, ProgressCallback *progress) |
| bool | isValidValue (float v, GDALRasterBand *band) |
| float | getInterpolatedValue (GDALRasterBand *band, double x, double y, bool applyOffset=true) |
| osg::HeightField * | createHeightField (const TileKey &key, ProgressCallback *progress) |
| bool | intersects (const TileKey &key) |
Static Public Member Functions | |
| static GDALRasterBand * | findBand (GDALDataset *ds, GDALColorInterp colorInterp) |
| static void | getPalleteIndexColor (GDALRasterBand *band, int index, osg::Vec4ub &color) |
Private Attributes | |
| GDALDataset * | _srcDS |
| GDALDataset * | _warpedDS |
| double | _geotransform [6] |
| double | _invtransform [6] |
| osg::Vec2d | _extentsMin |
| osg::Vec2d | _extentsMax |
| const GDALOptions | _options |
| unsigned int | _maxDataLevel |
Definition at line 611 of file ReaderWriterGDAL.cpp.
| GDALTileSource::GDALTileSource | ( | const TileSourceOptions & | options | ) | [inline] |
Definition at line 614 of file ReaderWriterGDAL.cpp.
:
TileSource( options ),
_srcDS(NULL),
_warpedDS(NULL),
_options(options),
_maxDataLevel(30)
{
}
| virtual GDALTileSource::~GDALTileSource | ( | ) | [inline, virtual] |
Definition at line 623 of file ReaderWriterGDAL.cpp.
| osg::HeightField* GDALTileSource::createHeightField | ( | const TileKey & | key, |
| ProgressCallback * | progress | ||
| ) | [inline, virtual] |
Creates a heightfield for the given TileKey The returned object is new and is the responsibility of the caller.
Reimplemented from osgEarth::TileSource.
Definition at line 1459 of file ReaderWriterGDAL.cpp.
{
if (key.getLevelOfDetail() > _maxDataLevel)
{
//OE_NOTICE << "Reached maximum data resolution key=" << key.getLevelOfDetail() << " max=" << _maxDataLevel << std::endl;
return NULL;
}
GDAL_SCOPED_LOCK;
int tileSize = _options.tileSize().value();
//Allocate the heightfield
osg::ref_ptr<osg::HeightField> hf = new osg::HeightField;
hf->allocate(tileSize, tileSize);
if (intersects(key))
{
//Get the meter extents of the tile
double xmin, ymin, xmax, ymax;
key.getExtent().getBounds(xmin, ymin, xmax, ymax);
//Just read from the first band
GDALRasterBand* band = _warpedDS->GetRasterBand(1);
double dx = (xmax - xmin) / (tileSize-1);
double dy = (ymax - ymin) / (tileSize-1);
for (int c = 0; c < tileSize; ++c)
{
double geoX = xmin + (dx * (double)c);
for (int r = 0; r < tileSize; ++r)
{
double geoY = ymin + (dy * (double)r);
float h = getInterpolatedValue(band, geoX, geoY);
hf->setHeight(c, r, h);
}
}
}
else
{
for (unsigned int i = 0; i < hf->getHeightList().size(); ++i) hf->getHeightList()[i] = NO_DATA_VALUE;
}
return hf.release();
}
Here is the call graph for this function:| osg::Image* GDALTileSource::createImage | ( | const TileKey & | key, |
| ProgressCallback * | progress | ||
| ) | [inline, virtual] |
Creates an image for the given TileKey. The returned object is new and is the responsibility of the caller.
Implements osgEarth::TileSource.
Definition at line 1019 of file ReaderWriterGDAL.cpp.
{
if (key.getLevelOfDetail() > _maxDataLevel)
{
OE_DEBUG << LC << "" << getName() << ": Reached maximum data resolution key="
<< key.getLevelOfDetail() << " max=" << _maxDataLevel << std::endl;
return NULL;
}
GDAL_SCOPED_LOCK;
int tileSize = _options.tileSize().value();
osg::ref_ptr<osg::Image> image;
if (intersects(key)) //TODO: I think this test is OBE -gw
{
//Get the extents of the tile
double xmin, ymin, xmax, ymax;
key.getExtent().getBounds(xmin, ymin, xmax, ymax);
double dx = (xmax - xmin) / (tileSize-1);
double dy = (ymax - ymin) / (tileSize-1);
int target_width = tileSize;
int target_height = tileSize;
int tile_offset_left = 0;
int tile_offset_top = 0;
int off_x = int((xmin - _geotransform[0]) / _geotransform[1]);
int off_y = int((ymax - _geotransform[3]) / _geotransform[5]);
int width = int(((xmax - _geotransform[0]) / _geotransform[1]) - off_x);
int height = int(((ymin - _geotransform[3]) / _geotransform[5]) - off_y);
if (off_x + width > _warpedDS->GetRasterXSize())
{
int oversize_right = off_x + width - _warpedDS->GetRasterXSize();
target_width = target_width - int(float(oversize_right) / width * target_width);
width = _warpedDS->GetRasterXSize() - off_x;
}
if (off_x < 0)
{
int oversize_left = -off_x;
tile_offset_left = int(float(oversize_left) / width * target_width);
target_width = target_width - int(float(oversize_left) / width * target_width);
width = width + off_x;
off_x = 0;
}
if (off_y + height > _warpedDS->GetRasterYSize())
{
int oversize_bottom = off_y + height - _warpedDS->GetRasterYSize();
target_height = target_height - (int)osg::round(float(oversize_bottom) / height * target_height);
height = _warpedDS->GetRasterYSize() - off_y;
}
if (off_y < 0)
{
int oversize_top = -off_y;
tile_offset_top = int(float(oversize_top) / height * target_height);
target_height = target_height - int(float(oversize_top) / height * target_height);
height = height + off_y;
off_y = 0;
}
OE_DEBUG << LC << "ReadWindow " << width << "x" << height << " DestWindow " << target_width << "x" << target_height << std::endl;
//Return if parameters are out of range.
if (width <= 0 || height <= 0 || target_width <= 0 || target_height <= 0)
{
return 0;
}
GDALRasterBand* bandRed = findBand(_warpedDS, GCI_RedBand);
GDALRasterBand* bandGreen = findBand(_warpedDS, GCI_GreenBand);
GDALRasterBand* bandBlue = findBand(_warpedDS, GCI_BlueBand);
GDALRasterBand* bandAlpha = findBand(_warpedDS, GCI_AlphaBand);
GDALRasterBand* bandGray = findBand(_warpedDS, GCI_GrayIndex);
GDALRasterBand* bandPalette = findBand(_warpedDS, GCI_PaletteIndex);
//The pixel format is always RGBA to support transparency
GLenum pixelFormat = GL_RGBA;
if (bandRed && bandGreen && bandBlue)
{
unsigned char *red = new unsigned char[target_width * target_height];
unsigned char *green = new unsigned char[target_width * target_height];
unsigned char *blue = new unsigned char[target_width * target_height];
unsigned char *alpha = new unsigned char[target_width * target_height];
//Initialize the alpha values to 255.
memset(alpha, 255, target_width * target_height);
image = new osg::Image;
image->allocateImage(tileSize, tileSize, 1, pixelFormat, GL_UNSIGNED_BYTE);
memset(image->data(), 0, image->getImageSizeInBytes());
//Nearest interpolation just uses RasterIO to sample the imagery and should be very fast.
if (!*_options.interpolateImagery() || _options.interpolation() == INTERP_NEAREST)
{
bandRed->RasterIO(GF_Read, off_x, off_y, width, height, red, target_width, target_height, GDT_Byte, 0, 0);
bandGreen->RasterIO(GF_Read, off_x, off_y, width, height, green, target_width, target_height, GDT_Byte, 0, 0);
bandBlue->RasterIO(GF_Read, off_x, off_y, width, height, blue, target_width, target_height, GDT_Byte, 0, 0);
if (bandAlpha)
{
bandAlpha->RasterIO(GF_Read, off_x, off_y, width, height, alpha, target_width, target_height, GDT_Byte, 0, 0);
}
for (int src_row = 0, dst_row = tile_offset_top;
src_row < target_height;
src_row++, dst_row++)
{
for (int src_col = 0, dst_col = tile_offset_left;
src_col < target_width;
++src_col, ++dst_col)
{
*(image->data(dst_col, dst_row) + 0) = red[src_col + src_row * target_width];
*(image->data(dst_col, dst_row) + 1) = green[src_col + src_row * target_width];
*(image->data(dst_col, dst_row) + 2) = blue[src_col + src_row * target_width];
*(image->data(dst_col, dst_row) + 3) = alpha[src_col + src_row * target_width];
}
}
image->flipVertical();
}
else
{
//Sample each point exactly
for (unsigned int c = 0; c < tileSize; ++c)
{
double geoX = xmin + (dx * (double)c);
for (unsigned int r = 0; r < tileSize; ++r)
{
double geoY = ymin + (dy * (double)r);
*(image->data(c,r) + 0) = getInterpolatedValue(bandRed, geoX,geoY,false);
*(image->data(c,r) + 1) = getInterpolatedValue(bandGreen,geoX,geoY,false);
*(image->data(c,r) + 2) = getInterpolatedValue(bandBlue, geoX,geoY,false);
if (bandAlpha != NULL)
*(image->data(c,r) + 3) = getInterpolatedValue(bandAlpha,geoX, geoY, false);
else
*(image->data(c,r) + 3) = 255;
}
}
}
delete []red;
delete []green;
delete []blue;
delete []alpha;
}
else if (bandGray)
{
unsigned char *gray = new unsigned char[target_width * target_height];
unsigned char *alpha = new unsigned char[target_width * target_height];
//Initialize the alpha values to 255.
memset(alpha, 255, target_width * target_height);
image = new osg::Image;
image->allocateImage(tileSize, tileSize, 1, pixelFormat, GL_UNSIGNED_BYTE);
memset(image->data(), 0, image->getImageSizeInBytes());
if (!*_options.interpolateImagery() || _options.interpolation() == INTERP_NEAREST)
{
bandGray->RasterIO(GF_Read, off_x, off_y, width, height, gray, target_width, target_height, GDT_Byte, 0, 0);
if (bandAlpha)
{
bandAlpha->RasterIO(GF_Read, off_x, off_y, width, height, alpha, target_width, target_height, GDT_Byte, 0, 0);
}
for (int src_row = 0, dst_row = tile_offset_top;
src_row < target_height;
src_row++, dst_row++)
{
for (int src_col = 0, dst_col = tile_offset_left;
src_col < target_width;
++src_col, ++dst_col)
{
*(image->data(dst_col, dst_row) + 0) = gray[src_col + src_row * target_width];
*(image->data(dst_col, dst_row) + 1) = gray[src_col + src_row * target_width];
*(image->data(dst_col, dst_row) + 2) = gray[src_col + src_row * target_width];
*(image->data(dst_col, dst_row) + 3) = alpha[src_col + src_row * target_width];
}
}
image->flipVertical();
}
else
{
for (int c = 0; c < tileSize; ++c)
{
double geoX = xmin + (dx * (double)c);
for (int r = 0; r < tileSize; ++r)
{
double geoY = ymin + (dy * (double)r);
float color = getInterpolatedValue(bandGray,geoX,geoY,false);
*(image->data(c,r) + 0) = color;
*(image->data(c,r) + 1) = color;
*(image->data(c,r) + 2) = color;
if (bandAlpha != NULL)
*(image->data(c,r) + 3) = getInterpolatedValue(bandAlpha,geoX,geoY,false);
else
*(image->data(c,r) + 3) = 255;
}
}
}
delete []gray;
delete []alpha;
}
else if (bandPalette)
{
//Pallete indexed imagery doesn't support interpolation currently and only uses nearest
//b/c interpolating pallete indexes doesn't make sense.
unsigned char *palette = new unsigned char[target_width * target_height];
image = new osg::Image;
image->allocateImage(tileSize, tileSize, 1, pixelFormat, GL_UNSIGNED_BYTE);
memset(image->data(), 0, image->getImageSizeInBytes());
bandPalette->RasterIO(GF_Read, off_x, off_y, width, height, palette, target_width, target_height, GDT_Byte, 0, 0);
for (int src_row = 0, dst_row = tile_offset_top;
src_row < target_height;
src_row++, dst_row++)
{
for (int src_col = 0, dst_col = tile_offset_left;
src_col < target_width;
++src_col, ++dst_col)
{
osg::Vec4ub color;
getPalleteIndexColor( bandPalette, palette[src_col + src_row * target_width], color );
*(image->data(dst_col, dst_row) + 0) = color.r();
*(image->data(dst_col, dst_row) + 1) = color.g();
*(image->data(dst_col, dst_row) + 2) = color.b();
*(image->data(dst_col, dst_row) + 3) = color.a();
}
}
image->flipVertical();
delete [] palette;
}
else
{
OE_WARN
<< LC << "Could not find red, green and blue bands or gray bands in "
<< _options.url()->full()
<< ". Cannot create image. " << std::endl;
return NULL;
}
}
// Moved this logic up into ImageLayer::createImageWrapper.
//if (!image.valid())
//{
// //OE_WARN << LC << "Illegal state-- should not get here" << std::endl;
// return ImageUtils::createEmptyImage();
//}
return image.release();
}
Here is the call graph for this function:| static GDALRasterBand* GDALTileSource::findBand | ( | GDALDataset * | ds, |
| GDALColorInterp | colorInterp | ||
| ) | [inline, static] |
Finds a raster band based on color interpretation
Definition at line 924 of file ReaderWriterGDAL.cpp.
{
GDAL_SCOPED_LOCK;
for (int i = 1; i <= ds->GetRasterCount(); ++i)
{
if (ds->GetRasterBand(i)->GetColorInterpretation() == colorInterp) return ds->GetRasterBand(i);
}
return 0;
}
| float GDALTileSource::getInterpolatedValue | ( | GDALRasterBand * | band, |
| double | x, | ||
| double | y, | ||
| bool | applyOffset = true |
||
| ) | [inline] |
Definition at line 1328 of file ReaderWriterGDAL.cpp.
{
double r, c;
GDALApplyGeoTransform(_invtransform, x, y, &c, &r);
//Account for slight rounding errors. If we are right on the edge of the dataset, clamp to the edge
double eps = 0.0001;
if (osg::equivalent(c, 0, eps)) c = 0;
if (osg::equivalent(r, 0, eps)) r = 0;
if (osg::equivalent(c, (double)_warpedDS->GetRasterXSize(), eps)) c = _warpedDS->GetRasterXSize();
if (osg::equivalent(r, (double)_warpedDS->GetRasterYSize(), eps)) r = _warpedDS->GetRasterYSize();
if (applyOffset)
{
//Apply half pixel offset
r-= 0.5;
c-= 0.5;
//Account for the half pixel offset in the geotransform. If the pixel value is -0.5 we are still technically in the dataset
//since 0,0 is now the center of the pixel. So, if are within a half pixel above or a half pixel below the dataset just use
//the edge values
if (c < 0 && c >= -0.5)
{
c = 0;
}
else if (c > _warpedDS->GetRasterXSize()-1 && c <= _warpedDS->GetRasterXSize()-0.5)
{
c = _warpedDS->GetRasterXSize()-1;
}
if (r < 0 && r >= -0.5)
{
r = 0;
}
else if (r > _warpedDS->GetRasterYSize()-1 && r <= _warpedDS->GetRasterYSize()-0.5)
{
r = _warpedDS->GetRasterYSize()-1;
}
}
float result = 0.0f;
//If the location is outside of the pixel values of the dataset, just return 0
if (c < 0 || r < 0 || c > _warpedDS->GetRasterXSize()-1 || r > _warpedDS->GetRasterYSize()-1)
return NO_DATA_VALUE;
if ( _options.interpolation() == INTERP_NEAREST )
{
band->RasterIO(GF_Read, (int)osg::round(c), (int)osg::round(r), 1, 1, &result, 1, 1, GDT_Float32, 0, 0);
if (!isValidValue( result, band))
{
return NO_DATA_VALUE;
}
}
else
{
int rowMin = osg::maximum((int)floor(r), 0);
int rowMax = osg::maximum(osg::minimum((int)ceil(r), (int)(_warpedDS->GetRasterYSize()-1)), 0);
int colMin = osg::maximum((int)floor(c), 0);
int colMax = osg::maximum(osg::minimum((int)ceil(c), (int)(_warpedDS->GetRasterXSize()-1)), 0);
if (rowMin > rowMax) rowMin = rowMax;
if (colMin > colMax) colMin = colMax;
float urHeight, llHeight, ulHeight, lrHeight;
band->RasterIO(GF_Read, colMin, rowMin, 1, 1, &llHeight, 1, 1, GDT_Float32, 0, 0);
band->RasterIO(GF_Read, colMin, rowMax, 1, 1, &ulHeight, 1, 1, GDT_Float32, 0, 0);
band->RasterIO(GF_Read, colMax, rowMin, 1, 1, &lrHeight, 1, 1, GDT_Float32, 0, 0);
band->RasterIO(GF_Read, colMax, rowMax, 1, 1, &urHeight, 1, 1, GDT_Float32, 0, 0);
/*
if (!isValidValue(urHeight, band)) urHeight = 0.0f;
if (!isValidValue(llHeight, band)) llHeight = 0.0f;
if (!isValidValue(ulHeight, band)) ulHeight = 0.0f;
if (!isValidValue(lrHeight, band)) lrHeight = 0.0f;
*/
if (!isValidValue(urHeight, band) || (!isValidValue(llHeight, band)) ||(!isValidValue(ulHeight, band)) || (!isValidValue(lrHeight, band)))
{
return NO_DATA_VALUE;
}
if ( _options.interpolation() == INTERP_AVERAGE )
{
double x_rem = c - (int)c;
double y_rem = r - (int)r;
double w00 = (1.0 - y_rem) * (1.0 - x_rem) * (double)llHeight;
double w01 = (1.0 - y_rem) * x_rem * (double)lrHeight;
double w10 = y_rem * (1.0 - x_rem) * (double)ulHeight;
double w11 = y_rem * x_rem * (double)urHeight;
result = (float)(w00 + w01 + w10 + w11);
}
else if ( _options.interpolation() == INTERP_BILINEAR )
{
//Check for exact value
if ((colMax == colMin) && (rowMax == rowMin))
{
//OE_NOTICE << "Exact value" << std::endl;
result = llHeight;
}
else if (colMax == colMin)
{
//OE_NOTICE << "Vertically" << std::endl;
//Linear interpolate vertically
result = ((float)rowMax - r) * llHeight + (r - (float)rowMin) * ulHeight;
}
else if (rowMax == rowMin)
{
//OE_NOTICE << "Horizontally" << std::endl;
//Linear interpolate horizontally
result = ((float)colMax - c) * llHeight + (c - (float)colMin) * lrHeight;
}
else
{
//OE_NOTICE << "Bilinear" << std::endl;
//Bilinear interpolate
float r1 = ((float)colMax - c) * llHeight + (c - (float)colMin) * lrHeight;
float r2 = ((float)colMax - c) * ulHeight + (c - (float)colMin) * urHeight;
//OE_INFO << "r1, r2 = " << r1 << " , " << r2 << std::endl;
result = ((float)rowMax - r) * r1 + (r - (float)rowMin) * r2;
}
}
}
return result;
}
| static void GDALTileSource::getPalleteIndexColor | ( | GDALRasterBand * | band, |
| int | index, | ||
| osg::Vec4ub & | color | ||
| ) | [inline, static] |
Definition at line 935 of file ReaderWriterGDAL.cpp.
{
const GDALColorEntry *colorEntry = band->GetColorTable()->GetColorEntry( index );
GDALPaletteInterp interp = band->GetColorTable()->GetPaletteInterpretation();
if (!colorEntry)
{
//FIXME: What to do here?
//OE_INFO << "NO COLOR ENTRY FOR COLOR " << rawImageData[i] << std::endl;
color.r() = 255;
color.g() = 0;
color.b() = 0;
color.a() = 1;
}
else
{
if (interp == GPI_RGB)
{
color.r() = colorEntry->c1;
color.g() = colorEntry->c2;
color.b() = colorEntry->c3;
color.a() = colorEntry->c4;
}
else if (interp == GPI_CMYK)
{
// from wikipedia.org
short C = colorEntry->c1;
short M = colorEntry->c2;
short Y = colorEntry->c3;
short K = colorEntry->c4;
color.r() = 255 - C*(255 - K) - K;
color.g() = 255 - M*(255 - K) - K;
color.b() = 255 - Y*(255 - K) - K;
color.a() = 255;
}
else if (interp == GPI_HLS)
{
// from easyrgb.com
float H = colorEntry->c1;
float S = colorEntry->c3;
float L = colorEntry->c2;
float R, G, B;
if ( S == 0 ) //HSL values = 0 - 1
{
R = L; //RGB results = 0 - 1
G = L;
B = L;
}
else
{
float var_2, var_1;
if ( L < 0.5 )
var_2 = L * ( 1 + S );
else
var_2 = ( L + S ) - ( S * L );
var_1 = 2 * L - var_2;
R = Hue_2_RGB( var_1, var_2, H + ( 1 / 3 ) );
G = Hue_2_RGB( var_1, var_2, H );
B = Hue_2_RGB( var_1, var_2, H - ( 1 / 3 ) );
}
color.r() = static_cast<unsigned char>(R*255.0f);
color.g() = static_cast<unsigned char>(G*255.0f);
color.b() = static_cast<unsigned char>(B*255.0f);
color.a() = static_cast<unsigned char>(255.0f);
}
else if (interp == GPI_Gray)
{
color.r() = static_cast<unsigned char>(colorEntry->c1*255.0f);
color.g() = static_cast<unsigned char>(colorEntry->c1*255.0f);
color.b() = static_cast<unsigned char>(colorEntry->c1*255.0f);
color.a() = static_cast<unsigned char>(255.0f);
}
}
}
Here is the call graph for this function:| void GDALTileSource::initialize | ( | const std::string & | referenceURI, |
| const Profile * | overrideProfile | ||
| ) | [inline, virtual] |
Initialize the TileSource. The profile should be computed and set here using setProfile()
Implements osgEarth::TileSource.
Definition at line 640 of file ReaderWriterGDAL.cpp.
{
GDAL_SCOPED_LOCK;
if ( !_options.url().isSet() || _options.url()->empty() )
{
OE_WARN << LC << "No URL or directory specified " << std::endl;
return;
}
URI uri = _options.url().value();
//Find the full path to the URL
//If we have a relative path and the map file contains a server address, just concat the server path and the _url together
if (osgEarth::isRelativePath(uri.full()) && osgDB::containsServerAddress(referenceURI))
{
uri = URI(osgDB::getFilePath(referenceURI) + std::string("/") + uri.full());
}
//If the path doesn't contain a server address, get the full path to the file.
if (!osgDB::containsServerAddress(uri.full()))
{
uri = URI(uri.full(), referenceURI);
}
StringTokenizer izer( ";" );
StringVector exts;
izer.tokenize( *_options.extensions(), exts );
//std::vector<std::string> exts;
//tokenize( _options.extensions().value(), exts, ";");
for (unsigned int i = 0; i < exts.size(); ++i)
{
OE_DEBUG << LC << "Using Extension: " << exts[i] << std::endl;
}
std::vector<std::string> files;
getFiles(uri.full(), exts, files);
OE_INFO << LC << "Driver found " << files.size() << " files:" << std::endl;
for (unsigned int i = 0; i < files.size(); ++i)
{
OE_INFO << LC << "" << files[i] << std::endl;
}
if (files.empty())
{
OE_WARN << LC << "Could not find any valid files " << std::endl;
return;
}
//If we found more than one file, try to combine them into a single logical dataset
if (files.size() > 1)
{
_srcDS = (GDALDataset*)build_vrt(files, HIGHEST_RESOLUTION);
if (!_srcDS)
{
OE_WARN << "[osgEarth::GDAL] Failed to build VRT from input datasets" << std::endl;
return;
}
}
else
{
//If we couldn't build a VRT, just try opening the file directly
//Open the dataset
_srcDS = (GDALDataset*)GDALOpen( files[0].c_str(), GA_ReadOnly );
if ( !_srcDS )
{
OE_WARN << LC << "Failed to open dataset " << files[0] << std::endl;
return;
}
}
//Create a spatial reference for the source.
const char* srcProj = _srcDS->GetProjectionRef();
if ( srcProj != 0L && overrideProfile != 0L )
{
OE_WARN << LC << "WARNING, overriding profile of a source that already defines its own SRS ("
<< this->getName() << ")" << std::endl;
}
osg::ref_ptr<const SpatialReference> src_srs;
if ( overrideProfile )
{
src_srs = overrideProfile->getSRS();
}
else if ( srcProj )
{
src_srs = SpatialReference::create( srcProj );
}
// assert SRS is present
if ( !src_srs.valid() )
{
// not found in the dataset; try loading a .prj file
std::string prjLocation = osgDB::getNameLessExtension( uri.full() ) + std::string(".prj");
std::string wkt;
if ( HTTPClient::readString( prjLocation, wkt ) == HTTPClient::RESULT_OK )
{
src_srs = SpatialReference::create( wkt );
}
if ( !src_srs.valid() )
{
OE_WARN << LC << "Dataset has no spatial reference information: " << uri.full() << std::endl;
return;
}
}
const Profile* profile = NULL;
if ( overrideProfile )
{
profile = overrideProfile;
}
if ( !profile && src_srs->isGeographic() )
{
profile = osgEarth::Registry::instance()->getGlobalGeodeticProfile();
}
//Note: Can cause odd rendering artifacts if we have a dataset that is mercator that doesn't encompass the whole globe
// if we take on the global profile.
/*
if ( !profile && src_srs->isMercator() )
{
profile = osgEarth::Registry::instance()->getGlobalMercatorProfile();
}*/
std::string warpedSRSWKT;
if ( profile && !profile->getSRS()->isEquivalentTo( src_srs.get() ) )
{
if ( profile->getSRS()->isGeographic() && (src_srs->isNorthPolar() || src_srs->isSouthPolar()) )
{
_warpedDS = (GDALDataset*)GDALAutoCreateWarpedVRTforPolarStereographic(
_srcDS,
src_srs->getWKT().c_str(),
profile->getSRS()->getWKT().c_str(),
GRA_NearestNeighbour,
5.0,
NULL);
}
else
{
_warpedDS = (GDALDataset*)GDALAutoCreateWarpedVRT(
_srcDS,
src_srs->getWKT().c_str(),
profile->getSRS()->getWKT().c_str(),
GRA_NearestNeighbour,
5.0,
NULL);
}
if ( _warpedDS )
{
warpedSRSWKT = _warpedDS->GetProjectionRef();
}
//GDALAutoCreateWarpedVRT(srcDS, src_wkt.c_str(), t_srs.c_str(), GRA_NearestNeighbour, 5.0, NULL);
}
else
{
_warpedDS = _srcDS;
warpedSRSWKT = src_srs->getWKT();
}
//Get the _geotransform
if (overrideProfile)
{
_geotransform[0] = overrideProfile->getExtent().xMin(); //Top left x
_geotransform[1] = overrideProfile->getExtent().width() / (double)_warpedDS->GetRasterXSize();//pixel width
_geotransform[2] = 0;
_geotransform[3] = overrideProfile->getExtent().yMax(); //Top left y
_geotransform[4] = 0;
_geotransform[5] = -overrideProfile->getExtent().height() / (double)_warpedDS->GetRasterYSize();//pixel height
}
else
{
_warpedDS->GetGeoTransform(_geotransform);
}
GDALInvGeoTransform(_geotransform, _invtransform);
//Compute the extents
// polar needs a special case when combined with geographic
if ( profile && profile->getSRS()->isGeographic() && (src_srs->isNorthPolar() || src_srs->isSouthPolar()) )
{
double ll_lon, ll_lat, ul_lon, ul_lat, ur_lon, ur_lat, lr_lon, lr_lat;
pixelToGeo(0.0, 0.0, ul_lon, ul_lat );
pixelToGeo(0.0, _warpedDS->GetRasterYSize(), ll_lon, ll_lat);
pixelToGeo(_warpedDS->GetRasterXSize(), _warpedDS->GetRasterYSize(), lr_lon, lr_lat);
pixelToGeo(_warpedDS->GetRasterXSize(), 0.0, ur_lon, ur_lat);
_extentsMin.x() = osg::minimum( ll_lon, osg::minimum( ul_lon, osg::minimum( ur_lon, lr_lon ) ) );
_extentsMax.x() = osg::maximum( ll_lon, osg::maximum( ul_lon, osg::maximum( ur_lon, lr_lon ) ) );
if ( src_srs->isNorthPolar() )
{
_extentsMin.y() = osg::minimum( ll_lat, osg::minimum( ul_lat, osg::minimum( ur_lat, lr_lat ) ) );
_extentsMax.y() = 90.0;
}
else
{
_extentsMin.y() = -90.0;
_extentsMax.y() = osg::maximum( ll_lat, osg::maximum( ul_lat, osg::maximum( ur_lat, lr_lat ) ) );
}
}
else
{
pixelToGeo(0.0, _warpedDS->GetRasterYSize(), _extentsMin.x(), _extentsMin.y());
pixelToGeo(_warpedDS->GetRasterXSize(), 0.0, _extentsMax.x(), _extentsMax.y());
}
OE_INFO << LC << "Geo extents: " << _extentsMin.x() << ", " << _extentsMin.y() << " => " << _extentsMax.x() << ", " << _extentsMax.y() << std::endl;
if ( !profile )
{
profile = Profile::create(
warpedSRSWKT,
//_warpedDS->GetProjectionRef(),
_extentsMin.x(), _extentsMin.y(), _extentsMax.x(), _extentsMax.y() );
OE_INFO << LC << "" << uri.full() << " is projected, SRS = "
<< warpedSRSWKT << std::endl;
//<< _warpedDS->GetProjectionRef() << std::endl;
}
//Compute the min and max data levels
double resolutionX = (_extentsMax.x() - _extentsMin.x()) / (double)_warpedDS->GetRasterXSize();
double resolutionY = (_extentsMax.y() - _extentsMin.y()) / (double)_warpedDS->GetRasterYSize();
double maxResolution = osg::minimum(resolutionX, resolutionY);
OE_INFO << LC << "Resolution= " << resolutionX << "x" << resolutionY << " max=" << maxResolution << std::endl;
if (_options.maxDataLevel().isSet())
{
_maxDataLevel = _options.maxDataLevel().value();
OE_INFO << "Using override max data level " << _maxDataLevel << std::endl;
}
else
{
unsigned int max_level = 30;
for (unsigned int i = 0; i < max_level; ++i)
{
_maxDataLevel = i;
double w, h;
profile->getTileDimensions(i, w, h);
double resX = (w / (double)_options.tileSize().value() );
double resY = (h / (double)_options.tileSize().value() );
if (resX < maxResolution || resY < maxResolution)
{
break;
}
}
OE_INFO << LC << "Max Data Level: " << _maxDataLevel << std::endl;
}
// record the data extent in profile space:
GeoExtent local_extent(
SpatialReference::create( warpedSRSWKT ), //_warpedDS->GetProjectionRef() ),
_extentsMin.x(), _extentsMin.y(), _extentsMax.x(), _extentsMax.y() );
GeoExtent profile_extent = local_extent.transform( profile->getSRS() );
getDataExtents().push_back( DataExtent(profile_extent, 0, _maxDataLevel) );
OE_INFO << LC << "Data Extents: " << profile_extent.toString() << std::endl;
//Set the profile
setProfile( profile );
}
Here is the call graph for this function:| bool GDALTileSource::intersects | ( | const TileKey & | key | ) | [inline] |
Definition at line 1506 of file ReaderWriterGDAL.cpp.
{
//Get the native extents of the tile
double xmin, ymin, xmax, ymax;
key.getExtent().getBounds(xmin, ymin, xmax, ymax);
return ! ( xmin >= _extentsMax.x() || xmax <= _extentsMin.x() || ymin >= _extentsMax.y() || ymax <= _extentsMin.y() );
}
Here is the call graph for this function:| bool GDALTileSource::isValidValue | ( | float | v, |
| GDALRasterBand * | band | ||
| ) | [inline] |
Definition at line 1299 of file ReaderWriterGDAL.cpp.
{
GDAL_SCOPED_LOCK;
float bandNoData = -32767.0f;
int success;
float value = band->GetNoDataValue(&success);
if (success)
{
bandNoData = value;
}
//Check to see if the value is equal to the bands specified no data
if (bandNoData == v) return false;
//Check to see if the value is equal to the user specified nodata value
if (getNoDataValue() == v) return false;
//Check to see if the user specified a custom min/max
if (v < getNoDataMinValue()) return false;
if (v > getNoDataMaxValue()) return false;
//Check within a sensible range
if (v < -32000) return false;
if (v > 32000) return false;
return true;
}
| void GDALTileSource::pixelToGeo | ( | double | x, |
| double | y, | ||
| double & | geoX, | ||
| double & | geoY | ||
| ) | [inline] |
Definition at line 1013 of file ReaderWriterGDAL.cpp.
{
geoX = _geotransform[0] + _geotransform[1] * x + _geotransform[2] * y;
geoY = _geotransform[3] + _geotransform[4] * x + _geotransform[5] * y;
}
osg::Vec2d GDALTileSource::_extentsMax [private] |
Definition at line 1524 of file ReaderWriterGDAL.cpp.
osg::Vec2d GDALTileSource::_extentsMin [private] |
Definition at line 1523 of file ReaderWriterGDAL.cpp.
double GDALTileSource::_geotransform[6] [private] |
Definition at line 1520 of file ReaderWriterGDAL.cpp.
double GDALTileSource::_invtransform[6] [private] |
Definition at line 1521 of file ReaderWriterGDAL.cpp.
unsigned int GDALTileSource::_maxDataLevel [private] |
Definition at line 1534 of file ReaderWriterGDAL.cpp.
const GDALOptions GDALTileSource::_options [private] |
Reimplemented from osgEarth::TileSource.
Definition at line 1531 of file ReaderWriterGDAL.cpp.
GDALDataset* GDALTileSource::_srcDS [private] |
Definition at line 1518 of file ReaderWriterGDAL.cpp.
GDALDataset* GDALTileSource::_warpedDS [private] |
Definition at line 1519 of file ReaderWriterGDAL.cpp.
1.7.3