sij/mapnik
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity Actions

Merge pull request #2686 from rouault/gdal_input_dataset_rasterio

Browse source 
Gdal.input dataset rasterio improvements
This commit is contained in:
Dane Springmeyer 2015-02-22 22:14:52 -08:00
commit 8919ae2f35
3 changed files with 85 additions and 52 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

64
plugins/input/gdal/gdal_datasource.cpp
View file

@ -45,37 +45,9 @@ using mapnik::layer_descriptor;
using mapnik::datasource_exception;
/*
* Opens a GDALDataset and returns a pointer to it.
* Caller is responsible for calling GDALClose on it
*/
inline GDALDataset* gdal_datasource::open_dataset() const
{
MAPNIK_LOG_DEBUG(gdal) << "gdal_datasource: Opening " << dataset_name_;
GDALDataset *dataset;
#if GDAL_VERSION_NUM >= 1600
if (shared_dataset_)
{
dataset = reinterpret_cast<GDALDataset*>(GDALOpenShared((dataset_name_).c_str(), GA_ReadOnly));
}
else
#endif
{
dataset = reinterpret_cast<GDALDataset*>(GDALOpen((dataset_name_).c_str(), GA_ReadOnly));
}
if (! dataset)
{
throw datasource_exception(CPLGetLastErrorMsg());
}
return dataset;
}
gdal_datasource::gdal_datasource(parameters const& params)
: datasource(params),
dataset_(nullptr),
desc_(gdal_datasource::name(), "utf-8"),
nodata_value_(params.get<double>("nodata")),
nodata_tolerance_(*params.get<double>("nodata_tolerance",1e-12))
@ -104,11 +76,27 @@ gdal_datasource::gdal_datasource(parameters const& params)
shared_dataset_ = *params.get<mapnik::boolean_type>("shared", false);
band_ = *params.get<mapnik::value_integer>("band", -1);
GDALDataset *dataset = open_dataset();
#if GDAL_VERSION_NUM >= 1600
if (shared_dataset_)
{
dataset_ = reinterpret_cast<GDALDataset*>(GDALOpenShared((dataset_name_).c_str(), GA_ReadOnly));
}
else
#endif
{
dataset_ = reinterpret_cast<GDALDataset*>(GDALOpen((dataset_name_).c_str(), GA_ReadOnly));
}
nbands_ = dataset->GetRasterCount();
width_ = dataset->GetRasterXSize();
height_ = dataset->GetRasterYSize();
if (! dataset_)
{
throw datasource_exception(CPLGetLastErrorMsg());
}
MAPNIK_LOG_DEBUG(gdal) << "gdal_featureset: opened Dataset=" << dataset_;
nbands_ = dataset_->GetRasterCount();
width_ = dataset_->GetRasterXSize();
height_ = dataset_->GetRasterYSize();
desc_.add_descriptor(mapnik::attribute_descriptor("nodata", mapnik::Double));
double tr[6];
@ -140,7 +128,7 @@ gdal_datasource::gdal_datasource(parameters const& params)
}
else
{
if (dataset->GetGeoTransform(tr) != CPLE_None)
if (dataset_->GetGeoTransform(tr) != CPLE_None)
{
MAPNIK_LOG_DEBUG(gdal) << "gdal_datasource GetGeotransform failure gives="
<< tr[0] << "," << tr[1] << ","
@ -187,8 +175,6 @@ gdal_datasource::gdal_datasource(parameters const& params)
extent_.init(x0, y0, x1, y1);
}
GDALClose(dataset);
MAPNIK_LOG_DEBUG(gdal) << "gdal_datasource: Raster Size=" << width_ << "," << height_;
MAPNIK_LOG_DEBUG(gdal) << "gdal_datasource: Raster Extent=" << extent_;
@ -196,6 +182,8 @@ gdal_datasource::gdal_datasource(parameters const& params)
gdal_datasource::~gdal_datasource()
{
MAPNIK_LOG_DEBUG(gdal) << "gdal_featureset: Closing Dataset=" << dataset_;
GDALClose(dataset_);
}
datasource::datasource_t gdal_datasource::type() const
@ -232,7 +220,7 @@ featureset_ptr gdal_datasource::features(query const& q) const
gdal_query gq = q;
// TODO - move to std::make_shared, but must reduce # of args to <= 9
return featureset_ptr(new gdal_featureset(*open_dataset(),
return featureset_ptr(new gdal_featureset(*dataset_,
band_,
gq,
extent_,
@ -254,7 +242,7 @@ featureset_ptr gdal_datasource::features_at_point(coord2d const& pt, double tol)
gdal_query gq = pt;
// TODO - move to std::make_shared, but must reduce # of args to <= 9
return featureset_ptr(new gdal_featureset(*open_dataset(),
return featureset_ptr(new gdal_featureset(*dataset_,
band_,
gq,
extent_,

1
plugins/input/gdal/gdal_datasource.hpp
View file

@ -56,6 +56,7 @@ public:
mapnik::layer_descriptor get_descriptor() const;
private:
GDALDataset* open_dataset() const;
GDALDataset* dataset_;
mapnik::box2d<double> extent_;
std::string dataset_name_;
int band_;

72
plugins/input/gdal/gdal_featureset.cpp
View file

@ -80,7 +80,6 @@ gdal_featureset::~gdal_featureset()
{
MAPNIK_LOG_DEBUG(gdal) << "gdal_featureset: Closing Dataset=" << &dataset_;
GDALClose(&dataset_);
}
feature_ptr gdal_featureset::next()
@ -376,7 +375,10 @@ feature_ptr gdal_featureset::get_feature(mapnik::query const& q)
raster_nodata = red->GetNoDataValue(&raster_has_nodata);
GDALColorTable *color_table = red->GetColorTable();
bool has_nodata = nodata_value_ || raster_has_nodata;
if (has_nodata && !color_table)
// we can deduce the alpha channel from nodata in the Byte case
// by reusing the reading of R,G,B bands directly
if (has_nodata && !color_table && red->GetRasterDataType() != GDT_Byte)
{
double apply_nodata = nodata_value_ ? *nodata_value_ : raster_nodata;
// read the data in and create an alpha channel from the nodata values
@ -402,20 +404,62 @@ feature_ptr gdal_featureset::get_feature(mapnik::query const& q)
}
}
}
raster_io_error = red->RasterIO(GF_Read, x_off, y_off, width, height, image.getBytes() + 0,
image.width(), image.height(), GDT_Byte, 4, 4 * image.width());
if (raster_io_error == CE_Failure) {
throw datasource_exception(CPLGetLastErrorMsg());
/* Use dataset RasterIO in priority in 99.9% of the cases */
if( red->GetBand() == 1 && green->GetBand() == 2 && blue->GetBand() == 3 )
{
int nBandsToRead = 3;
if( alpha != NULL && alpha->GetBand() == 4 && !raster_has_nodata )
{
nBandsToRead = 4;
alpha = NULL; // to avoid reading it again afterwards
}
raster_io_error = dataset_.RasterIO(GF_Read, x_off, y_off, width, height,
image.getBytes(),
image.width(), image.height(), GDT_Byte,
nBandsToRead, NULL,
4, 4 * image.width(), 1);
if (raster_io_error == CE_Failure) {
throw datasource_exception(CPLGetLastErrorMsg());
}
}
raster_io_error = green->RasterIO(GF_Read, x_off, y_off, width, height, image.getBytes() + 1,
image.width(), image.height(), GDT_Byte, 4, 4 * image.width());
if (raster_io_error == CE_Failure) {
throw datasource_exception(CPLGetLastErrorMsg());
else
{
raster_io_error = red->RasterIO(GF_Read, x_off, y_off, width, height, image.getBytes() + 0,
image.width(), image.height(), GDT_Byte, 4, 4 * image.width());
if (raster_io_error == CE_Failure) {
throw datasource_exception(CPLGetLastErrorMsg());
}
raster_io_error = green->RasterIO(GF_Read, x_off, y_off, width, height, image.getBytes() + 1,
image.width(), image.height(), GDT_Byte, 4, 4 * image.width());
if (raster_io_error == CE_Failure) {
throw datasource_exception(CPLGetLastErrorMsg());
}
raster_io_error = blue->RasterIO(GF_Read, x_off, y_off, width, height, image.getBytes() + 2,
image.width(), image.height(), GDT_Byte, 4, 4 * image.width());
if (raster_io_error == CE_Failure) {
throw datasource_exception(CPLGetLastErrorMsg());
}
}
raster_io_error = blue->RasterIO(GF_Read, x_off, y_off, width, height, image.getBytes() + 2,
image.width(), image.height(), GDT_Byte, 4, 4 * image.width());
if (raster_io_error == CE_Failure) {
throw datasource_exception(CPLGetLastErrorMsg());
// In the case we skipped initializing the alpha channel
if (has_nodata && !color_table && red->GetRasterDataType() == GDT_Byte)
{
double apply_nodata = nodata_value_ ? *nodata_value_ : raster_nodata;
if( apply_nodata >= 0 && apply_nodata <= 255 )
{
int len = image.width() * image.height();
GByte* pabyBytes = (GByte*) image.getBytes();
for (int i = 0; i < len; ++i)
{
// TODO - we assume here the nodata value for the red band applies to all bands
// more details about this at http://trac.osgeo.org/gdal/ticket/2734
if (std::fabs(apply_nodata - pabyBytes[4*i]) < nodata_tolerance_)
pabyBytes[4*i + 3] = 0;
else
pabyBytes[4*i + 3] = 255;
}
}
}
}
else if (grey)