C++ IRasterCoverage::size方法代码示例

int OperationHelperRaster::subdivideTasks(ExecutionContext *ctx,const IRasterCoverage& raster, const BoundingBox &bnds, std::vector<BoundingBox > &boxes)
{
    if ( !raster.isValid() || raster->size().isNull() || raster->size().ysize() == 0) {
        return ERROR1(ERR_NO_INITIALIZED_1, "Grid size");
        return iUNDEF;
    }

    int cores = std::min(QThread::idealThreadCount(),(int)raster->size().ysize());
    if (raster->size().linearSize() < 10000 || ctx->_threaded == false)
        cores = 1;

    boxes.clear();
    boxes.resize(cores);
    BoundingBox bounds = bnds;
    if ( bounds.isNull())
        bounds = BoundingBox(raster->size());
    int left = 0; //bounds.min_corner().x;
    int right = bounds.size().xsize();
    int top = bounds.size().ysize();
    int step = bounds.size().ysize() / cores;
    int currentY = 0;

    for(int i=0 ; i < cores; ++i){
        BoundingBox smallBox(Pixel(left, currentY,0), Pixel(right - 1, std::min(top - 1,currentY + step),bounds.zlength()) );
        boxes[i] = smallBox;
        currentY = currentY + step  ;
    }
    return cores;
}

BoundingBox OperationHelperRaster::initialize(const IRasterCoverage &inputRaster, IRasterCoverage &outputRaster, quint64 what)
{
    Resource resource(itRASTER);
    Size<> sz = inputRaster->size();
    BoundingBox box(sz);

    if ( what & itRASTERSIZE) {
        resource.addProperty("size", IVARIANT(sz.toString()));
    }
    if ( what & itENVELOPE) {
        if ( box.isNull() || !box.isValid()) {
            sz = inputRaster->size();
            box  = BoundingBox(sz);
        }
        Envelope bounds = inputRaster->georeference()->pixel2Coord(box);
        resource.addProperty("envelope", IVARIANT(bounds.toString()));
    }
    if ( what & itCOORDSYSTEM) {
        resource.addProperty("coordinatesystem", IVARIANT(inputRaster->coordinateSystem()->id()));
    }

    if ( what & itGEOREF) {
        if ( box.isNull() || !box.isValid()) {
            sz = inputRaster->size();
            box  = BoundingBox(sz);
        }
        if ( sz.xsize() == box.xlength() && sz.ysize() == box.ylength())
            resource.addProperty("georeference", IVARIANT(inputRaster->georeference()->id()));
    }
    if ( what & itDOMAIN) {
        resource.addProperty("domain", IVARIANT(inputRaster->datadef().domain<>()->id()));
    }
    resource.prepare();

    outputRaster.prepare(resource);
    if ( what & itTABLE) {
        if ( inputRaster->attributeTable().isValid())    {
            if ( inputRaster->datadef().domain<>() == outputRaster->datadef().domain<>()) {
                if ( outputRaster.isValid())
                    outputRaster->setAttributes(inputRaster->attributeTable());
            }
        }
    }
    if ( what & itDOMAIN){
        for(quint32 i = 0; i < outputRaster->size().zsize(); ++i){
            QString index = outputRaster->stackDefinition().index(i);
            outputRaster->setBandDefinition(index,DataDefinition(outputRaster->datadef().domain()));
        }
    }


    return box;
}

Ilwis::OperationImplementation::State PercentileFilterStretch::prepare(ExecutionContext *, const SymbolTable & )
{
    QString raster = _expression.parm(0).value();
    QString outputName = _expression.parm(0,false).value();

    // open the input map list
    if (!_inputObj.prepare(raster, itRASTER)) {
        ERROR2(ERR_COULD_NOT_LOAD_2, raster, "");
        return sPREPAREFAILED;
    }

    // Open the zonal map
    QString zones = _expression.parm(1).value();
    if (!_inputZones.prepare(zones, itRASTER)) {
        ERROR2(ERR_COULD_NOT_LOAD_2, zones, "");
        return sPREPAREFAILED;
    }

    // Open both percentile tables
    QString lowPerc = _expression.parm(2).value();
    if (!_lowPercentile.prepare(lowPerc, {"mustexist", true})){
        ERROR2(ERR_COULD_NOT_LOAD_2, lowPerc, "");
        return sPREPAREFAILED;
    }

    QString highPerc = _expression.parm(3).value();
    if (!_highPercentile.prepare(highPerc, {"mustexist", true})){
        ERROR2(ERR_COULD_NOT_LOAD_2, highPerc, "");
        return sPREPAREFAILED;
    }

    QString startDate = _expression.parm(4).value();
    Time date(startDate);
    _startDekad = dekadFromDate(date);
    if (!date.isValid()) {
        ERROR1(TR("Invalid date: "), startDate);
        return sPREPAREFAILED;
    }

    _outputObj = OperationHelperRaster::initialize(_inputObj, itRASTER, itDOMAIN | itGEOREF | itCOORDSYSTEM | itRASTERSIZE | itBOUNDINGBOX | itENVELOPE);
    if ( !_outputObj.isValid()) {
        ERROR1(ERR_NO_INITIALIZED_1, "output rastercoverage");
        return sPREPAREFAILED;
    }

    IRasterCoverage inputRaster = _inputObj.as<RasterCoverage>();
    // initialize tranquilizer
    initialize(inputRaster->size().xsize() * inputRaster->size().ysize());

    return sPREPARED;
}

IIlwisObject OperationHelperRaster::initialize(const IIlwisObject &inputObject, IlwisTypes tp, quint64 what)
{
    Resource resource(tp);
    if (inputObject->ilwisType() & itCOVERAGE) {
        ICoverage cov = inputObject.as<Coverage>();
        if (inputObject->ilwisType() == itRASTER) {
            IRasterCoverage gcInput = inputObject.as<RasterCoverage>();
            if ( what & itRASTERSIZE) {
                Size<> sz = gcInput->size();
                BoundingBox box(sz);
                resource.addProperty("size", IVARIANT(box.size()));
            }
            if ( what & itGEOREF) {
                resource.addProperty("georeference", IVARIANT(gcInput->georeference()));
            }
            if ( what & itDOMAIN) {
                resource.addProperty("domain", IVARIANT(gcInput->datadef().domain()));
            }
        }
        if ( what & itCOORDSYSTEM) {
            resource.addProperty("coordinatesystem", IVARIANT(cov->coordinateSystem()));
        }

     }

    resource.prepare();
    IIlwisObject obj;
    obj.prepare(resource);
    if (inputObject->ilwisType() & itCOVERAGE) {
        OperationHelper::initialize(inputObject, obj, tp, what);
    }

    return obj;
}

bool PinDataSource::setStackDomain(const QString& id) {
	IRasterCoverage raster;
	raster.prepare(_objid);

	bool ok;
	_actives.clear();
	quint64 oid = id.toULongLong(&ok);
	if (!ok) {
		kernel()->issues()->log(TR("No valid object id used: ") + id);
		return false;
	}
	RasterStackDefinition&  stack = raster->stackDefinitionRef();

	IDomain dom;
	if (dom.prepare(oid) && dom->ilwisType() == itITEMDOMAIN) {
		_stackDomain = dom;
		IItemDomain itemdom = dom.as< ItemDomain<DomainItem>>();
		if (itemdom->count() == raster->size().zsize()) {
			std::vector<QString> items;
			for (auto item : itemdom) {
				items.push_back(item->name());
			}
			stack.setSubDefinition(itemdom, items);
		}
		else {
			kernel()->issues()->log(TR("Item domain must have same size as tjhe number of bands in the container"));
			return false;
		}
	}

	for (quint32 i = 0; i < stack.count(); ++i) {
		QString name = QString::number(i + 1); // stack.index(i+1);
		QVariantMap data;
		data["name"] = name;
		data["active"] = true;
		if (_stackDomain->valueType() != itNUMERICITEM) {
			data["minvalue"] = "?";
			data["maxvalue"] = "?";
		}
		else {
			IIntervalDomain idomain = _stackDomain.as<IntervalDomain>();
			if (idomain->count() == stack.count()) {
				auto item = idomain->item(i)->as<Interval>();
				data["name"] = item->name();
				auto range = item->range().as<NumericRange>();
				data["minvalue"] = QString::number(range->min());
				data["maxvalue"] = QString::number(range->max());
			}
		}
		_actives.push_back(data);
	}
	return true;
}

BoundingBox SelectionBase::boundingBox(const IRasterCoverage& raster) const
{
    BoundingBox box;
    std::set<QString> bands;
    for(const auto& epart : _expressionparts){
        box += epart._box;
        if ( epart._bands.size() != 0){
            for(QString b : epart._bands)
                bands.insert(b) ;
        }
    }
    if ( bands.size() != 0){
        if ( !box.isValid())
            box = raster->size();
        auto &p = box.max_corner();
        p.z = bands.size();
    }
    return box;
}

void RasterCoverage::copyBinary(const IRasterCoverage& raster, quint32 inputIndex, quint32 outputIndex) {
    if ( isNumericalUndef(inputIndex) || isNumericalUndef(outputIndex)){
        ERROR2(ERR_ILLEGAL_VALUE_2,TR("layer index"), isNumericalUndef(inputIndex) ? "input" : "output");
        return;
    }
    if ( inputIndex >= size().zsize()){
       //ERROR2(ERR_ILLEGAL_VALUE_2,TR("layer index"), "input");
    }
    IRasterCoverage gcNew;
    gcNew.set(this);
    Size<> inputSize =  raster->size();
    Size<> sz(inputSize.xsize(),inputSize.ysize(), outputIndex + 1);
    gcNew->georeference()->size(sz);
    PixelIterator iterIn(raster, BoundingBox(Pixel(0,0,inputIndex), Pixel(inputSize.xsize(), inputSize.ysize(), inputIndex + 1)));
    PixelIterator iterOut(gcNew, BoundingBox(Pixel(0,0,outputIndex), Pixel(inputSize.xsize(), inputSize.ysize(), outputIndex + 1)));
    if ( raster->id() == id() && inputIndex == outputIndex){
        ERROR2(ERR_ILLEGALE_OPERATION2, TR("copy"),TR("identical layers in same raster"));
        return;
    }
    std::for_each(iterOut, iterOut.end(), [&](double& v){
         v = *iterIn;
        ++iterIn;
    });
}

bool SelectionRaster::execute(ExecutionContext *ctx, SymbolTable& symTable)
{
    if (_prepState == sNOTPREPARED)
        if((_prepState = prepare(ctx, symTable)) != sPREPARED)
            return false;
    IRasterCoverage outputRaster = _outputObj.as<RasterCoverage>();
    IRasterCoverage inputRaster = _inputObj.as<RasterCoverage>();

    std::map<Raw, quint32> raw2record;
    int keyColumn = _inputAttributeTable.isValid() ? _inputAttributeTable->columnIndex(inputRaster->primaryKey()) : iUNDEF;
    if (keyColumn != iUNDEF){
        std::vector<QVariant> values = _inputAttributeTable->column(keyColumn);
        for(quint32 rec=0; rec < values.size(); ++rec){
            Raw r = values[rec].toDouble();
            if ( !isNumericalUndef(r)){
                raw2record[r] = rec;
            }
        }
    }

    std::vector<int> extraAtrrib = organizeAttributes();


    std::vector<QString> selectionBands = bands(inputRaster);
    initialize(outputRaster->size().linearSize() * selectionBands.size());

    PixelIterator iterOut(outputRaster);
    int count = 0;
    bool numeric = outputRaster->datadef().domain()->ilwisType() == itNUMERICDOMAIN;
    for(QString band : selectionBands){
        PixelIterator iterIn = inputRaster->band(band, _box);

        PixelIterator iterEnd = iterIn.end();
        while(iterIn != iterEnd) {
            bool ok = true;
            double pixValue = *iterIn;
            double matchValue = pixValue;

            for(const auto& epart : _expressionparts){
                bool partOk = epart.match(iterIn.position(), matchValue,this);
                if ( epart._andor != loNONE)
                    ok =  epart._andor == loAND ? ok && partOk : ok || partOk;
                else
                    ok &= partOk;
                if (epart._type == ExpressionPart::ptATTRIBUTE && extraAtrrib.size() == 1){
                    if ( pixValue < 0 || pixValue >= _inputAttributeTable->recordCount()){
                        ok = false;
                        continue;
                    }
                    // pixValue == ID; ID < zero means undef, ID's start at zero.
                    if (pixValue >= 0) {
                        if (keyColumn != iUNDEF){
                            auto iter = raw2record.find(pixValue);
                            if ( iter != raw2record.end()){
                                quint32 rec = iter->second;
                                const Record& record = _inputAttributeTable->record(rec);
                                pixValue = record.cell(extraAtrrib[0]).toDouble();
                            }else
                                pixValue = rUNDEF;
                        }

                    }
                    else
                        pixValue = rUNDEF;
                }
            }
            if ( ok){
                *iterOut = pixValue;
            }else
                *iterOut = rUNDEF;

            ++iterIn;
            ++iterOut;
            updateTranquilizer(++count, 100);
        }
        // if there is an attribute table we must copy the correct attributes and records
        if ( keyColumn != iUNDEF && _attTable.isValid()){
            for(int recIndex=0; recIndex < _inputAttributeTable->recordCount(); ++recIndex){
                const Record& rec = _inputAttributeTable->record(recIndex);
                for(int i=0; i < extraAtrrib.size(); ++i){
                    _attTable->setCell(i, recIndex, rec.cell(extraAtrrib[i]));
                }
            }
        }
    }
    if ( numeric)
        outputRaster->statistics(ContainerStatistics<double>::pBASIC);

    outputRaster->setAttributes(_attTable);
    QVariant value;
    value.setValue<IRasterCoverage>(outputRaster);
    logOperation(outputRaster, _expression);
    ctx->setOutput(symTable, value, outputRaster->name(), itRASTER,outputRaster->resource());
    return true;


}

Ilwis::OperationImplementation::State SelectionRaster::prepare(ExecutionContext *ctx, const SymbolTable &st)
{
    OperationImplementation::prepare(ctx,st);
    if ( _expression.parameterCount() != 2) {
        ERROR3(ERR_ILLEGAL_NUM_PARM3,"rasvalue","1",QString::number(_expression.parameterCount()));
        return sPREPAREFAILED;
    }
    IlwisTypes inputType = itRASTER;
    QString raster = _expression.parm(0).value();
    if (!_inputObj.prepare(raster, inputType)) {
        ERROR2(ERR_COULD_NOT_LOAD_2,raster,"");
        return sPREPAREFAILED;
    }
    IRasterCoverage inputRaster = _inputObj.as<RasterCoverage>();
    _inputAttributeTable = inputRaster->attributeTable();
    quint64 copylist = itCOORDSYSTEM ;


    QString selector = _expression.parm(1).value();
    parseSelector(selector, inputRaster);

    std::vector<QString> selectionBands = bands(inputRaster);
     _box = boundingBox(_inputObj.as<RasterCoverage>());
    bool useOldGrf ;
    if ( _box.isNull() || !_box.isValid()){
        _box = inputRaster->size();
        useOldGrf = true;
    } else
        useOldGrf = false;

    if ( useOldGrf){
        copylist |= itGEOREF | itRASTERSIZE | itENVELOPE;
    }


    int selectedAttributes = attributeNames().size();
    if (selectedAttributes != 1)
        copylist |= itDOMAIN;

     _outputObj = OperationHelperRaster::initialize(_inputObj,inputType, copylist);
     if ( !_outputObj.isValid()) {
         ERROR1(ERR_NO_INITIALIZED_1, "output coverage");
         return sPREPAREFAILED;
     }
     IRasterCoverage outputRaster = _outputObj.as<RasterCoverage>();

     QString outputName = _expression.parm(0,false).value();
     if ( outputName != sUNDEF)
         _outputObj->name(outputName);

     if ( selectedAttributes > 1) {
         QString url = INTERNAL_CATALOG + "/" + outputName;
         Resource resource(url, itFLATTABLE);
         _attTable.prepare(resource);
     }
     if ( selectedAttributes == 1 && _inputAttributeTable.isValid()){
        QStringList names = attributeNames();
        //outputRaster->datadefRef().domain(_inputAttributeTable->columndefinition(names[0]).datadef().domain());
        outputRaster->setDataDefintions(_inputAttributeTable->columndefinition(names[0]).datadef().domain(), selectionBands, inputRaster->stackDefinition().domain());
     }else
         outputRaster->setDataDefintions(inputRaster->datadef().domain(), selectionBands, inputRaster->stackDefinition().domain());



     if ( (copylist & itGEOREF) == 0) {
         Resource resource(QUrl(INTERNAL_CATALOG + "/" + outputRaster->name() + "_grf_" + QString::number(outputRaster->id())),itGEOREF);
         resource.addProperty("size", IVARIANT(_box.size()));
         auto envelope = inputRaster->georeference()->pixel2Coord(_box);
         resource.addProperty("envelope", IVARIANT(envelope));
         resource.addProperty("coordinatesystem", IVARIANT(inputRaster->coordinateSystem()));
         resource.addProperty("name", _outputObj->name());
         resource.addProperty("centerofpixel",inputRaster->georeference()->centerOfPixel());
         IGeoReference  grf;
         grf.prepare(resource);
         outputRaster->georeference(grf);
         outputRaster->envelope(envelope);
    }


    return sPREPARED;
}

Ilwis::OperationImplementation::State Selection::prepare(ExecutionContext *, const SymbolTable &)
{
    if ( _expression.parameterCount() != 2) {
        ERROR3(ERR_ILLEGAL_NUM_PARM3,"rasvalue","1",QString::number(_expression.parameterCount()));
        return sPREPAREFAILED;
    }
    IlwisTypes inputType = itRASTER;
    QString raster = _expression.parm(0).value();
    if (!_inputObj.prepare(raster, inputType)) {
        ERROR2(ERR_COULD_NOT_LOAD_2,raster,"");
        return sPREPAREFAILED;
    }
    IRasterCoverage inputRaster = _inputObj.as<RasterCoverage>();
    quint64 copylist = itCOORDSYSTEM;


    QString selector = _expression.parm(1).value();
    selector = selector.remove('"');

    int index = selector.indexOf("box=");
    Envelope box;
    if ( index != -1) {
        QString crdlist = "box(" + selector.mid(index+4) + ")";
        _box = BoundingBox(crdlist);
        box = inputRaster->georeference()->pixel2Coord(_box);
        copylist |= itDOMAIN | itTABLE;
        std::vector<qint32> vec{_box.min_corner().x, _box.min_corner().y,_box.min_corner().z};
        _base = vec;

    }
    index = selector.indexOf("polygon=");
    if ( index != -1)
    {
        //TODO:
        copylist |= itDOMAIN | itTABLE;
    }
    index = selector.indexOf("attribute=");
    if ( index != -1 ) {
        if (! inputRaster->attributeTable().isValid()) {
            ERROR2(ERR_NO_FOUND2,"attribute-table", "coverage");
            return sPREPAREFAILED;
        }
        _attribColumn =  selector.mid(index+10);
        copylist |= itRASTERSIZE | itGEOREF | itENVELOPE;
    }
    int indexindex = selector.indexOf("index=");
    if ( indexindex != -1) {
        copylist |= itDOMAIN | itGEOREF | itENVELOPE | itTABLE;
        _box = BoundingBox(inputRaster->size());
        QString zvalues = selector.mid(6);
        bool ok;
        _zvalue = zvalues.toInt(&ok);
        if ( !ok || _zvalue < 0) {
            ERROR3(ERR_ILLEGAL_PARM_3, TR("layer index"), zvalues,"Selection");
            return sPREPAREFAILED;
        }
        _box.min_corner().z = _zvalue;
        _box.max_corner().z = _zvalue;
        std::vector<qint32> vec{_box.min_corner().x, _box.min_corner().y,_box.min_corner().z};
        _base = vec;
    }

     _outputObj = OperationHelperRaster::initialize(_inputObj,inputType, copylist);
     if ( !_outputObj.isValid()) {
         ERROR1(ERR_NO_INITIALIZED_1, "output coverage");
         return sPREPAREFAILED;
     }
     IRasterCoverage outputRaster = _outputObj.as<RasterCoverage>();
     if ( (copylist & itDOMAIN) == 0) {
         outputRaster->datadefRef() = _attribColumn != "" ? inputRaster->attributeTable()->columndefinition(_attribColumn).datadef()
                                                   : outputRaster->datadefRef() = inputRaster->datadef();
     }
     QString outputName = _expression.parm(0,false).value();
     if ( outputName != sUNDEF)
         _outputObj->name(outputName);
     if ( (copylist & itGEOREF) == 0) {
        Resource resource(QUrl("ilwis://internalcatalog/georeference"),itGEOREF);
        resource.addProperty("size", IVARIANT(_box.size()));
        resource.addProperty("envelope", IVARIANT(box));
        resource.addProperty("coordinatesystem", IVARIANT(inputRaster->coordinateSystem()));
        resource.addProperty("name", _outputObj->name());
        resource.addProperty("centerofpixel",inputRaster->georeference()->centerOfPixel());
        IGeoReference  grf;
        grf.prepare(resource);
        outputRaster->georeference(grf);
        outputRaster->envelope(box);
    }
     if(indexindex != -1)  {
         Size<> sz(outputRaster->size().xsize(),outputRaster->size().xsize(), 1);
         outputRaster->size(sz);
     }

    return sPREPARED;
}

LayerModel *LayerManager::create(QStandardItem *parentLayer, const ICoverage &cov, LayerManager *lm, const IOOptions &options)
{
    if (cov->coordinateSystem()->isUnknown() && lm->rootLayer()->screenCsy().isValid()){
        QString mes = QString("coordinate system 'unknown' not compatible with coordinate system of the layerview");
        kernel()->issues()->log(mes, IssueObject::itWarning);
        return 0;
    }
    if (parentLayer == 0)
        parentLayer = lm->layerTree()->invisibleRootItem();

	QString type = TypeHelper::type2name(cov->ilwisType());
	if (options.contains("createtype")) {
		type = options["createtype"].toString();
	}
    auto iter = _createLayers.find(type);
    if ( iter != _createLayers.end()){
        auto createFunc   = (*iter).second;
        QString layername = options.contains("layername") ? options["layername"].toString(): cov->name();
        LayerModel *layer = createFunc(lm, parentLayer,layername, cov->description(),options);
		if (!lm->rootLayer()->screenCsy().isValid()) {// first real layer sets the csy
			lm->rootLayer()->screenCsy(cov->coordinateSystem());
            lm->rootLayer()->coverageEnvelope(cov->envelope());
            if (cov->ilwisType() == itRASTER) {
                IRasterCoverage rc = cov.as<RasterCoverage>();
                if (rc->georeference()->grfType<UndeterminedGeoReference>()) {
                    Envelope envUndetermned(Coordinate(0, 0), Coordinate(rc->size().xsize(), rc->size().ysize())); // none.grf bounds
                    lm->rootLayer()->coverageEnvelope(envUndetermned);
                }
            }
	  		lm->setLayerListName(layername);
		}
        else {
            // adjust rootlayer envelop to fit all layers
            Coordinate crd1_trans;
            Coordinate crd2_trans;
            Envelope envelop = cov->envelope();
            if (!cov->coordinateSystem()->isEqual(lm->rootLayer()->screenCsy().ptr())) {
                Coordinate crn1 = envelop.min_corner();
                Coordinate crn2 = envelop.max_corner();

                crd1_trans = lm->rootLayer()->screenCsy()->coord2coord(cov->coordinateSystem(), crn1);
                crd2_trans = lm->rootLayer()->screenCsy()->coord2coord(cov->coordinateSystem(), crn2);
            }
            else {
                crd1_trans = envelop.min_corner();
                crd2_trans = envelop.max_corner();
            }
            Envelope orgenv = lm->rootLayer()->coverageEnvelope();
            if (!(orgenv.contains(crd1_trans) && orgenv.contains(crd2_trans))) {
                orgenv += crd1_trans;
                orgenv += crd2_trans;

                lm->rootLayer()->coverageEnvelope(orgenv);
                if (cov->ilwisType() == itRASTER) {
                    // TODO: check georef none ?
                }
            }
        }
        qint32 lowernodeid = options.contains("lowernodeid") ? options["lowernodeid"].toInt() : iUNDEF;
		layer->nodeId(lm->nextId());
		layer->fillData();
        addLayer(parentLayer, layer, lm, lowernodeid);

		lm->lastAddedCoverageLayer(layer);
    
        return layer;
    }
    return 0;

}

bool PercentileFilterStretch::execute(ExecutionContext *ctx, SymbolTable& symTable)
{
    if (_prepState == sNOTPREPARED)
        if((_prepState = prepare(ctx,symTable)) != sPREPARED)
            return false;

    IRasterCoverage outputRaster = _outputObj.as<RasterCoverage>();
    IRasterCoverage inputRaster = _inputObj.as<RasterCoverage>();
    IRasterCoverage zoneRaster = _inputZones.as<RasterCoverage>();

    ITable low = _lowPercentile.as<Table>();
    ITable high = _highPercentile.as<Table>();

    PixelIterator iterIn(inputRaster, BoundingBox(), PixelIterator::fZXY);
    PixelIterator iterZone(zoneRaster, BoundingBox(), PixelIterator::fXYZ); // only one layer so Z is irrelevant
    PixelIterator iterOut(outputRaster, BoundingBox(), PixelIterator::fZXY);
    PixelIterator inEnd = iterIn.end();

    _nb = inputRaster->size().zsize();
    int rows = inputRaster->size().xsize();
    int cols = inputRaster->size().ysize();
    std::vector<double> slice(_nb);
    std::vector<int> percentage(_nb);
    int totalRows = low->recordCount(); // same as high->recordCount()
    int totalCols = low->columnCount(); // same as high->columnCount()
    std::vector<double> lowtab(low->columnCount() * low->recordCount());
    std::vector<double> hightab(high->columnCount() * high->recordCount());
    for (int row = 0; row < totalRows; ++row)
        for (int col = 0; col < totalCols; ++col) {
            // Let the first column in the percentile table match the start of time series
            int actCol = (col + totalCols - _startDekad) % totalCols;
            lowtab[actCol + row * totalCols] = low->cell(col, row).toDouble();
            hightab[actCol + row * totalCols] = high->cell(col, row).toDouble();
        }

    // timeseries are assumed to be 10 day periods.
    int pixCount = 0;
    while (iterIn != inEnd) {
        trq()->update(pixCount++);

        // get the time slice at the current location
        std::copy(iterIn, iterIn + _nb, slice.begin());
        // get the zone at the current location
        double dzone = *iterZone;     // row index into low and high percentile tables
        if (dzone == rUNDEF) {
            // for out of area locations set to zero percent
            std::fill(percentage.begin(), percentage.end(), 0);
        }
        else {
            int zone = (long) dzone;

            std::vector<double>::const_iterator liter = lowtab.begin() + zone * totalCols;
            std::vector<double>::const_iterator hiter = hightab.begin() + zone * totalCols;
            std::vector<int>::iterator piter = percentage.begin();
            for (std::vector<double>::const_iterator siter = slice.begin(); siter != slice.end(); siter++) {
                *piter = std::max(0, std::min(100, int(100.0 * (*hiter - *siter) / (*hiter - *liter))));
                if (*piter <= 5) *piter = 0;
                else if (*piter <= 10) *piter = 10;
                piter++;
                liter++;
                hiter++;
            }
        }
        std::copy(percentage.begin(), percentage.end(), iterOut);
        iterIn += _nb;
        iterOut += _nb;
        iterZone += 1;
    }

    trq()->update(rows * cols);
    trq()->inform("\nWriting...\n");
    trq()->stop();

    bool resource = true;
    if ( resource && ctx != 0) {
        QVariant value;
        value.setValue<IRasterCoverage>(outputRaster);
        ctx->setOutput(symTable, value, outputRaster->name(), itRASTER, outputRaster->resource() );
    }
    return resource;
}

Ilwis::OperationImplementation::State AggregateRaster::prepare(ExecutionContext *, const SymbolTable & )
{
    QString raster = _expression.parm(0).value();
    QString outputName = _expression.parm(0,false).value();
    int copylist = itDOMAIN | itCOORDSYSTEM;

    if (!_inputObj.prepare(raster, itRASTER)) {
        ERROR2(ERR_COULD_NOT_LOAD_2,raster,"");
        return sPREPAREFAILED;
    }
    _method = toMethod(_expression.parm(1).value());
    if ( _method == NumericStatistics::pLAST) {
        ERROR2(ERR_ILLEGAL_VALUE_2, "parameter value", " aggregation method");
        return sPREPAREFAILED;
    }
    bool ok;
    quint32 groupSz = _expression.parm(2).value().toInt(&ok);
    if (!ok) {
        QString blist = _expression.parm(2).value();
        blist.remove("{");
        blist.remove("}");
        QStringList dims = blist.split(" ");
        if ( dims.size() > 0) {
            for(int i=0; i < dims.size(); ++i)     {
                quint32 val = dims[i].toInt(&ok);
                if ( ok) {
                    _groupSize[i] = val;
                }else
                    break;
            }
        }

    }else {
        _groupSize[0] = groupSz;
        _groupSize[1] = groupSz;

    }
    if ( !ok || groupSize() < 2) {
        ERROR2(ERR_ILLEGAL_VALUE_2, "aggregation group size", QString::number(groupSize()));
        return sPREPAREFAILED;
    }

    _grouped = _expression.parm(3).value().toLower() == "true";
    if ( !_grouped)
        copylist |= itGEOREF;

    _outputObj = OperationHelperRaster::initialize(_inputObj,itRASTER, copylist);
    if ( !_outputObj.isValid()) {
        ERROR1(ERR_NO_INITIALIZED_1, "output rastercoverage");
        return sPREPAREFAILED;
    }
    QString outputBaseName = outputName;
    int index = 0;
    if ( (index = outputName.lastIndexOf(".")) != -1) {
        outputBaseName = outputName.left(index);
    }
    IRasterCoverage inputRaster = _inputObj.as<RasterCoverage>();
    IRasterCoverage outputRaster = _outputObj.as<RasterCoverage>();
    if ( outputName != sUNDEF)
        _outputObj->name(outputName);
    outputRaster->coordinateSystem(inputRaster->coordinateSystem());

    BoundingBox box(inputRaster->size());
    if ( _grouped) {
        int xs = box.xlength();
        int ys = box.ylength();
        int zs = box.zlength();
        int newxs = xs / groupSize(0);
        int newys = ys / groupSize(1);
        int newzs = zs / groupSize(2);
        box = BoundingBox(Size<>(newxs, newys, newzs));

    }
    if ( _expression.parameterCount() == 5 || _grouped) {
        Envelope envlope = inputRaster->envelope();
        Resource resource(QUrl("ilwis://internalcatalog/georeference"),itGEOREF);
        resource.addProperty("size", IVARIANT(box.size()));
        resource.addProperty("envelope", IVARIANT(envlope));
        resource.addProperty("coordinatesystem", inputRaster->coordinateSystem()->id());
        resource.addProperty("name", outputBaseName);
        resource.addProperty("centerofpixel",inputRaster->georeference()->centerOfPixel());
        IGeoReference  grf;
        if (grf.prepare(resource)) {
            mastercatalog()->addItems({resource});
            outputRaster->georeference(grf);
        }
        outputRaster->envelope(envlope);
        outputRaster->size(box.size());
    }

    initialize(outputRaster->size().linearSize());

    return sPREPARED;
}