Raster Algebra

[jericks]

Hi all,

I am working on my raster branch again:

https://github.com/jericks/geoscript-groovy/commits/raster

One of the things I would like the geoscript raster module to be able to do is map algebra.  My Raster class now has add, minus, divide, multiply methods that use the wonderful org.geotools.coverage.processing.operation classes and jiffle.  I noticed that GeoTools didn't have Subtract and Divide operations, so I took a stab at implementing them:  

https://github.com/jericks/geoscript-groovy/blob/raster/src/main/groovy/geoscript/raster/Raster.groovy  (at the way bottom)

Divide works fine but Subtract gives me an exception:

java.lang.IllegalArgumentException: Illegal transform of type "ConstantTransform1D".

at org.geotools.coverage.Category.<init>(Category.java:494)

at org.geotools.coverage.Category.<init>(Category.java:429)

at org.geotools.coverage.Category.<init>(Category.java:398)

at org.geotools.coverage.processing.OperationJAI.deriveCategory(OperationJAI.java:841)

at org.geotools.coverage.processing.OperationJAI.deriveSampleDimension(OperationJAI.java:792)

at org.geotools.coverage.processing.OperationJAI.deriveGridCoverage(OperationJAI.java:572)

at org.geotools.coverage.processing.OperationJAI.doOperation(OperationJAI.java:302)

at org.geotools.coverage.processing.CoverageProcessor.doOperation(CoverageProcessor.java:542)

at org.geotools.coverage.processing.CoverageProcessor.doOperation(CoverageProcessor.java:564)

at org.geotools.coverage.processing.CoverageProcessor$doOperation.call(Unknown Source)

    ...

Caused by: org.opengis.referencing.operation.NoninvertibleTransformException: Transform is not invertible.

at org.geotools.referencing.operation.transform.AbstractMathTransform.inverse(AbstractMathTransform.java:585)

at org.geotools.referencing.operation.transform.LinearTransform1D.inverse(LinearTransform1D.java:167)

at org.geotools.coverage.Category.<init>(Category.java:530)

at org.geotools.coverage.GeophysicsCategory.<init>(GeophysicsCategory.java:58)

at org.geotools.coverage.Category.<init>(Category.java:485)

...

when called:

def processor = new CoverageProcessor()

def params = processor.getOperation("Subtract").parameters

params.parameter("Source0").value = this.coverage

params.parameter("Source1").value = other.coverage

def newCoverage = processor.doOperation(params)

new Raster(newCoverage)

so I switched to Jiffle and it works:

MapAlgebra mapAlgebra = new MapAlgebra()

List size = this.size

Raster output = mapAlgebra.calculate("dest = src1 - src2;", [src1: this, src2: other])

output

Sorry for the long email, but I was hoping someone could help.  Is using the org.geotools.coverage.processing.operation classes the right way to go, or should I just use Jiffle?

Thanks,

Jared

[andrea antonello]

Hi Jared,
a quick guess/comment.

Each time I had to deal with ConstantTransform1D errors it was due to
the fact that the resulting map had one single value available and
when it came to styling the map to render it broke.
Is it possible that the mapalgebra instead works because it doesn't
handle novalues on its own and handles the novalue as a value? I.e. to
make that work properly with novalues you should have a check like

if(isNull(src)){
dest = null
}else{

dest = src1 - src2
}

I am going by memory on this.

Cheers,
Andrea

> --
> --
> You received this message because you are subscribed to the GeoScript
> mailing list.
> To post to this group, send email to geos...@googlegroups.com
> To unsubscribe from this group, send email to
> geoscript+...@googlegroups.com
> Visit this group at http://groups.google.com/group/geoscript or see
> http://geoscript.org
>
> ---
> You received this message because you are subscribed to the Google Groups
> "GeoScript" group.
> To unsubscribe from this group and stop receiving emails from it, send an
> email to geoscript+...@googlegroups.com.
> For more options, visit https://groups.google.com/groups/opt_out.
>
>

[Justin Deoliveira]

Hey Jared,

Nice work. I more or less have the same thing going on my raster branch:

  https://github.com/jdeolive/geoscript-py/blob/raster/geoscript/layer/raster.py#L290

Indeed I think using jiffle makes a whole lot of sense, i am just not sure in what context. If it makes sense for simple operations like this or more for full blown scripts that involve more complex operations. 

One thing I have been hoping to do is analyze the performance of using jiffle vs going through the CoverageProcessor.

-Justin

--
Justin Deoliveira

OpenGeo - http://opengeo.org

Enterprise support for open source geospatial.

[Andrea Aime]

On Wed, Feb 27, 2013 at 4:09 PM, Justin Deoliveira <jdeo...@opengeo.org> wrote:

Hey Jared,

Nice work. I more or less have the same thing going on my raster branch:

  https://github.com/jdeolive/geoscript-py/blob/raster/geoscript/layer/raster.py#L290

Indeed I think using jiffle makes a whole lot of sense, i am just not sure in what context. If it makes sense for simple operations like this or more for full blown scripts that involve more complex operations. 

One thing I have been hoping to do is analyze the performance of using jiffle vs going through the CoverageProcessor.

My gut feeling about it:

- if you are doing only a few operations CoverageProcessor should be faster, as it can rely on native JAI if installed

- if you have a complex expression or a full script, I'd expect Jiffle to be faster instead, since every JAI operation

creates a tile set that ends up filling the tile cache, when you have too many the tile cache will start to suffer and

the computation slow down. With Jiffle there is really just input and output, no intermediate tiles

Cheers

Andrea

 

--

==

Our support, Your Success! Visit http://opensdi.geo-solutions.it for more information.

==

Ing. Andrea Aime 

@geowolf

Technical Lead

GeoSolutions S.A.S.

Via Poggio alle Viti 1187

55054  Massarosa (LU)

Italy

phone: +39 0584 962313

fax: +39 0584 1660272

mob: +39  339 8844549

http://www.geo-solutions.it

http://twitter.com/geosolutions_it

-------------------------------------------------------

[Justin Deoliveira]

On Wed, Feb 27, 2013 at 8:13 AM, Andrea Aime <andre...@geo-solutions.it> wrote:

On Wed, Feb 27, 2013 at 4:09 PM, Justin Deoliveira <jdeo...@opengeo.org> wrote:

Hey Jared,

Nice work. I more or less have the same thing going on my raster branch:

  https://github.com/jdeolive/geoscript-py/blob/raster/geoscript/layer/raster.py#L290

Indeed I think using jiffle makes a whole lot of sense, i am just not sure in what context. If it makes sense for simple operations like this or more for full blown scripts that involve more complex operations. 

One thing I have been hoping to do is analyze the performance of using jiffle vs going through the CoverageProcessor.

My gut feeling about it:

- if you are doing only a few operations CoverageProcessor should be faster, as it can rely on native JAI if installed

- if you have a complex expression or a full script, I'd expect Jiffle to be faster instead, since every JAI operation

creates a tile set that ends up filling the tile cache, when you have too many the tile cache will start to suffer and

the computation slow down. With Jiffle there is really just input and output, no intermediate tiles

Nice! Thanks Andrea. That saves me hours of time :) 

Cheers

Andrea

 

--

==

Our support, Your Success! Visit http://opensdi.geo-solutions.it for more information.

==

Ing. Andrea Aime 

@geowolf

Technical Lead

GeoSolutions S.A.S.

Via Poggio alle Viti 1187

55054  Massarosa (LU)

Italy

phone: +39 0584 962313

fax: +39 0584 1660272

mob: +39  339 8844549

http://www.geo-solutions.it

http://twitter.com/geosolutions_it

-------------------------------------------------------

--

--
You received this message because you are subscribed to the GeoScript mailing list.
To post to this group, send email to geos...@googlegroups.com
To unsubscribe from this group, send email to geoscript+...@googlegroups.com
Visit this group at http://groups.google.com/group/geoscript or see http://geoscript.org
 
---
You received this message because you are subscribed to the Google Groups "GeoScript" group.
To unsubscribe from this group and stop receiving emails from it, send an email to geoscript+...@googlegroups.com.
For more options, visit https://groups.google.com/groups/opt_out.
 
 

[Jared Erickson]

Hi Andrea,

You were right on. My unit test that was throwing the ConstantTransform1D error had a raster/coverage that only had a single value. Thanks for the pointer, it saved me hours.

Jared