Skip to content

Commit

Permalink
Allow statistics to be selected via option.
Browse files Browse the repository at this point in the history
  • Loading branch information
@abellgithub
abellgithub committed Sep 26, 2016
1 parent 5e8875e commit 5954f31
Show file tree
Hide file tree
Showing 2 changed files with 175 additions and 59 deletions.
62 changes: 54 additions & 8 deletions io/gdal/GDALWriter.cpp
View file
Expand Up @@ -50,7 +50,6 @@ static PluginInfo const s_info = PluginInfo(

CREATE_STATIC_PLUGIN(1, 0, GDALWriter, Writer, s_info);


std::string GDALWriter::getName() const
{
return s_info.name;
Expand All @@ -67,11 +66,43 @@ void GDALWriter::addArgs(ProgramArgs& args)
args.add("gdaldriver", "GDAL writer driver name", m_drivername, "GTiff");
args.add("gdalopts", "GDAL driver options (name=value,name=value...)",
m_options);
args.add("output_type", "Statistics produced ('min', 'max', 'mean', "
"'idw', 'count', 'stdev' or 'all')", m_outputTypeString);
}


void GDALWriter::initialize()
{
for (auto& ts : m_outputTypeString)
{
Utils::trim(ts);
if (ts == "all")
{
m_outputTypes = ~0;
break;
}
if (ts == "min")
m_outputTypes |= Grid::statMin;
else if (ts == "max")
m_outputTypes |= Grid::statMax;
else if (ts == "count")
m_outputTypes |= Grid::statCount;
else if (ts == "mean")
m_outputTypes |= Grid::statMean;
else if (ts == "idw")
m_outputTypes |= Grid::statIdw;
else if (ts == "stdev")
m_outputTypes |= Grid::statStdDev;
else
{
std::ostringstream oss;
oss << "Invalid writers.gdal output type: '" << ts << "'.";
throw pdal_error(oss.str());
}
}
if (m_outputTypes == 0)
m_outputTypes = ~0;

gdal::registerDrivers();
}

Expand All @@ -93,7 +124,8 @@ void GDALWriter::write(const PointViewPtr view)
view->calculateBounds(m_bounds);
size_t width = ceil((m_bounds.maxx - m_bounds.minx) / m_edgeLength) + 1;
size_t height = ceil((m_bounds.maxy - m_bounds.miny) / m_edgeLength) + 1;
m_grid.reset(new Grid(width, height, m_edgeLength, m_radius, -9999.0));
m_grid.reset(new Grid(width, height, m_edgeLength, m_radius, -9999.0,
m_outputTypes));

for (PointId idx = 0; idx < view->size(); ++idx)
{
Expand Down Expand Up @@ -127,12 +159,26 @@ void GDALWriter::done(PointTableRef table)
m_grid->numBands(), Dimension::Type::Double, m_grid->noData());
if (err != gdal::GDALError::None)
throw pdal_error(raster.errorMsg());
raster.writeBand(m_grid->data("min"), 1, "min");
raster.writeBand(m_grid->data("max"), 2, "max");
raster.writeBand(m_grid->data("mean"), 3, "mean");
raster.writeBand(m_grid->data("idw"), 4, "idw");
raster.writeBand(m_grid->data("count"), 5, "count");
raster.writeBand(m_grid->data("stdev"), 6, "stdev");
int bandNum = 1;
uint8_t *buf;
buf = m_grid->data("min");
if (buf)
raster.writeBand(buf, bandNum++, "min");
buf = m_grid->data("max");
if (buf)
raster.writeBand(buf, bandNum++, "max");
buf = m_grid->data("mean");
if (buf)
raster.writeBand(buf, bandNum++, "mean");
buf = m_grid->data("idw");
if (buf)
raster.writeBand(buf, bandNum++, "idw");
buf = m_grid->data("count");
if (buf)
raster.writeBand(buf, bandNum++, "count");
buf = m_grid->data("stdev");
if (buf)
raster.writeBand(buf, bandNum++, "stdev");
}

} // namespace pdal
Expand Down
172 changes: 121 additions & 51 deletions io/gdal/GDALWriter.hpp
View file
Expand Up @@ -45,16 +45,37 @@ namespace pdal
class Grid
{
public:
static const int statCount = 1;
static const int statMin = 2;
static const int statMax = 4;
static const int statMean = 8;
static const int statStdDev = 16;
static const int statIdw = 32;

Grid(size_t width, size_t height, double edgeLength, double radius,
double noData) :
double noData, int outputTypes) :
m_width(width), m_height(height), m_edgeLength(edgeLength),
m_radius(radius), m_noData(noData),
m_count(width * height),
m_min(width * height, std::numeric_limits<double>::max()),
m_max(width * height, std::numeric_limits<double>::lowest()),
m_mean(width * height), m_stdDev(width * height),
m_idw(width * height), m_idwDist(width * height)
{}
m_radius(radius), m_noData(noData), m_outputTypes(outputTypes)
{
size_t size(width * height);

m_count.reset(new DataVec(size));
if (m_outputTypes & statMin)
m_min.reset(new DataVec(size,
std::numeric_limits<double>::max()));
if (m_outputTypes & statMax)
m_max.reset(new DataVec(size,
std::numeric_limits<double>::lowest()));
if (m_outputTypes & statIdw)
{
m_idw.reset(new DataVec(size));
m_idwDist.reset(new DataVec(size));
}
if ((m_outputTypes & statMean) || (m_outputTypes & statStdDev))
m_mean.reset(new DataVec(size));
if (m_outputTypes & statStdDev)
m_stdDev.reset(new DataVec(size));
}

int horizontalIndex(double x)
{ return x / m_edgeLength; }
Expand All @@ -76,23 +97,45 @@ class Grid
}

int numBands() const
{ return 6; }
{
int num = 0;

if (m_outputTypes & statCount)
num++;
if (m_outputTypes & statMin)
num++;
if (m_outputTypes & statMax)
num++;
if (m_outputTypes & statMean)
num++;
if (m_outputTypes & statIdw)
num++;
if (m_outputTypes & statStdDev)
num++;
return num;
}

uint8_t *data(const std::string& name)
{
if (name == "count")
return (uint8_t *)m_count.data();
return (m_outputTypes & statCount ?
(uint8_t *)m_count->data() : nullptr);
if (name == "min")
return (uint8_t *)m_min.data();
return (m_outputTypes & statMin ?
(uint8_t *)m_min->data() : nullptr);
if (name == "max")
return (uint8_t *)m_max.data();
return (m_outputTypes & statMax ?
(uint8_t *)m_max->data() : nullptr);
if (name == "mean")
return (uint8_t *)m_mean.data();
return (m_outputTypes & statMean ?
(uint8_t *)m_mean->data() : nullptr);
if (name == "idw")
return (uint8_t *)m_idw.data();
return (m_outputTypes & statIdw ?
(uint8_t *)m_idw->data() : nullptr);
if (name == "stdev")
return (uint8_t *)m_stdDev.data();
throw pdal_error("Requested invalid grid data '" + name + "'.");
return (m_outputTypes & statStdDev ?
(uint8_t *)m_stdDev->data() : nullptr);
return nullptr;
}

void addPoint(double x, double y, double z)
Expand Down Expand Up @@ -195,48 +238,68 @@ class Grid

size_t offset = (y * m_width) + x;

double& count = m_count[offset];
double& count = (*m_count)[offset];
count++;

double& min = m_min[offset];
min = std::min(val, min);

double& max = m_max[offset];
max = std::max(val, max);
if (m_min)
{
double& min = (*m_min)[offset];
min = std::min(val, min);
}

double& mean = m_mean[offset];
double delta = val - mean;
if (m_max)
{
double& max = (*m_max)[offset];
max = std::max(val, max);
}

mean += delta / count;
if (m_mean)
{
double& mean = (*m_mean)[offset];
double delta = val - mean;
if (m_stdDev)
{
mean += delta / count;

double& stdDev = m_stdDev[offset];
stdDev += delta * (val - mean);
double& stdDev = (*m_stdDev)[offset];
stdDev += delta * (val - mean);
}
}

double& idw = m_idw[offset];
idw += val / dist;
if (m_idw)
{
double& idw = (*m_idw)[offset];
idw += val / dist;

double& idwDist = m_idwDist[offset];
idwDist += 1 / dist;
double& idwDist = (*m_idwDist)[offset];
idwDist += 1 / dist;
}
}

void finalize()
{
// See
// https://en.wikipedia.org/wiki/Algorithms_for_calculating_variance
// https://en.wikipedia.org/wiki/Inverse_distance_weighting
for (size_t i = 0; i < m_count.size(); ++i)
{
if (m_count[i])
{
m_stdDev[i] = sqrt(m_stdDev[i] / m_count[i]);
m_idw[i] /= m_idwDist[i];
}
else
{
m_min[i] = m_noData;
m_max[i] = m_noData;
}
}
if (m_stdDev)
for (size_t i = 0; i < m_count->size(); ++i)
if ((*m_count)[i])
(*m_stdDev)[i] = sqrt((*m_stdDev)[i] / (*m_count)[i]);

if (m_idw)
for (size_t i = 0; i < m_count->size(); ++i)
if ((*m_count)[i])
(*m_idw)[i] /= (*m_idwDist)[i];

if (m_min)
for (size_t i = 0; i < m_count->size(); ++i)
if (!(*m_count)[i])
(*m_min)[i] = m_noData;

if (m_max)
for (size_t i = 0; i < m_count->size(); ++i)
if (!(*m_count)[i])
(*m_max)[i] = m_noData;
}

size_t width() const
Expand All @@ -252,13 +315,18 @@ class Grid
double m_edgeLength;
double m_radius;
double m_noData;
std::vector<double> m_count;
std::vector<double> m_min;
std::vector<double> m_max;
std::vector<double> m_mean;
std::vector<double> m_stdDev;
std::vector<double> m_idw;
std::vector<double> m_idwDist;

typedef std::vector<double> DataVec;
typedef std::unique_ptr<DataVec> DataPtr;
DataPtr m_count;
DataPtr m_min;
DataPtr m_max;
DataPtr m_mean;
DataPtr m_stdDev;
DataPtr m_idw;
DataPtr m_idwDist;

int m_outputTypes;
};
typedef std::unique_ptr<Grid> GridPtr;

Expand All @@ -281,6 +349,8 @@ class PDAL_DLL GDALWriter : public Writer
double m_edgeLength;
double m_radius;
StringList m_options;
StringList m_outputTypeString;
int m_outputTypes;
GridPtr m_grid;
};

Expand Down

0 comments on commit 5954f31

Please sign in to comment.