geotools-gt2-users

[Geotools-gt2-users] masking a grid coverage

[Michael B. <mic...@gm...>]

Hi all,

In my app I have two single band grid coverages, one with double
values and one with integer values.  For each grid cell in the first
coverage I want to calculate the sum of values within a given
distance, using the second coverage as a mask (ignore pixels with
value <= 0).  Finally I want to binarize the result of the convolution
with a given threshold value.

I'm wondering what the most efficient way to do this is.  At present I
have code working that does the above steps without masking. It makes
use of JAI's ConvolveDescriptor class to sum neighbourhood values,
then the ImageWorker class to binarize the result.

To add masking to the process I guess I start out by creating a JAI
ROI object from the RenderedImage derived from the second coverage.
But after that I'm a bit stuck as to how to use the ROI in the
convolution stage, other than the naive approach of creating a new
input image with, for instance, TiledImage.set(img, roi).

I wonder if anyone can suggest an approach ?  Also, am I overlooking
any Geotools class methods that would do some or all of the above
without resorting to JAI directly ?

Any help much appreciated.

Michael

Re: [Geotools-gt2-users] masking a grid coverage

[Martin D. <mar...@ge...>]

Hello Michael

ROI is one aspect of JAI that I lack experience in. The approach that I'm
writting below may be less efficient than ROI. I posting it just in case, but it
may be worth to explore more on the ROI side.

A possible approach may be to binarize the mask in order to create a new image
containing only 0 and 1 values. Note: I would not recommand to use ImageWorker
if you can avoid it. I would rather suggest to use the JAI Binarize descriptor
directly:

http://download.java.net/media/jai/javadoc/1.1.3/jai-apidocs/javax/media/jai/operator/BinarizeDescriptor.html

Then you can multiply your numeric image by that mask, so the masked values get
0. You can then convolved as you did, those 0 values will change the sum. For
the last binarize step, I would suggest again the JAI operation rather than
ImageWorker.

However the above involves a lot of images (up 6) because of intermediate steps.
Maybe this is not an issue. But if it is, you could also create your own JAI
operations doing all this work in one step. It also give you more control on
what is computed exactly. If you wish to do so, the first step would be to
create a subclass of javax.media.jai.AreaOpImage.

	Martin

Re: [Geotools-gt2-users] masking a grid coverage

[Michael B. <mic...@gm...>]

Hi Martin,

Many thanks for your answer - I was actually hoping that you would respond !

I haven't studied ROI use in any detail yet, but at first glance it
doesn't seem to be as integrated into JAI operations as I thought it
would be.  But I'll look into it more.

I think I'll try your suggestion of extending AreaOpImage to combine
the steps in the masking and convolution process.  At some stage I'll
be working with large grid coverages so avoiding the overhead of
intermediate images would be very useful.  I'll let you know how I go
!

As an aside, I created a sub-class of GridCoverage2D to make it easier
to work with data that change frequently.  It caches new pixel values
until the coverage is queried, or the cache is full, at which time it
uses a TiledImage with the same backing data buffer as the grid
coverage and a JAI iterator to update the coverage.  This has been
working very well so far - although I saw a warning (from Andrea ?)
about using this approach a while ago.  I think it might be useful to
add some of the grid calculation methods that I'm using to this class
also, in which case it might be more generally useful for other
Geotools users.  And instead of a sub-class it could also be made into
a decorator which might be more convenient sometimes.

All the best,
Michael

2009/1/7 Martin Desruisseaux <mar...@ge...>:
> Hello Michael
>
> ROI is one aspect of JAI that I lack experience in. The approach that I'm
> writting below may be less efficient than ROI. I posting it just in case, but it
> may be worth to explore more on the ROI side.
>
> A possible approach may be to binarize the mask in order to create a new image
> containing only 0 and 1 values. Note: I would not recommand to use ImageWorker
> if you can avoid it. I would rather suggest to use the JAI Binarize descriptor
> directly:
>
> http://download.java.net/media/jai/javadoc/1.1.3/jai-apidocs/javax/media/jai/operator/BinarizeDescriptor.html
>
> Then you can multiply your numeric image by that mask, so the masked values get
> 0. You can then convolved as you did, those 0 values will change the sum. For
> the last binarize step, I would suggest again the JAI operation rather than
> ImageWorker.
>
> However the above involves a lot of images (up 6) because of intermediate steps.
> Maybe this is not an issue. But if it is, you could also create your own JAI
> operations doing all this work in one step. It also give you more control on
> what is computed exactly. If you wish to do so, the first step would be to
> create a subclass of javax.media.jai.AreaOpImage.
>
>        Martin
>
>

Re: [Geotools-gt2-users] masking a grid coverage

[Michael B. <mic...@gm...>]

2009/1/7 Martin Desruisseaux wrote:
> However the above involves a lot of images (up 6) because of intermediate steps.
> Maybe this is not an issue. But if it is, you could also create your own JAI
> operations doing all this work in one step. It also give you more control on
> what is computed exactly. If you wish to do so, the first step would be to
> create a subclass of javax.media.jai.AreaOpImage.
>

Hi Martin,

I created a new operator, MaskedConvolveOpImage, together with
associated factory, descriptor and spi classes.  It is really just a
hack of JAI's ConvolveOpImage where I've added an ROI to the
constructor which is then queried (roi.contains method) when moving
the kernel and processing kernel cells.  It works very well :-)

Along the way I learned more about JAI's use of ROI but it does seem
to be quite limited.  Basically, the ROI only affects statistical
operators.  An ROI can be attached as a property to an image but
doesn't act as a mask for non-statistical image operations - rather it
is just kept in synch with the image if possible (e.g. cropped,
sub-sampled).  This puzzles me a little because I would have thought
there were many applications that would need selective processing of
pixels based on threshold value, juxtaposition etc.

Would the MaskedConvolve operator be a useful addition to the
gt-coverage module ?

Thanks again for your suggestion.
Michael

Re: [Geotools-gt2-users] masking a grid coverage

[Martin D. <mar...@ge...>]

Michael Bedward a écrit :
> Would the MaskedConvolve operator be a useful addition to the
> gt-coverage module ?

It could for images backed by integer values. For floating point values, I tend
to favor a different approach but I don't know if it is applicable in your case.

In my case, image backed by floating point values are measurements of some
geophysical quantity, for example an elevation in metres on a temperature in
degrees celcius. The areas to masks are areas where there is no data available.
I represent those data by Float.NaN value.

In some cases, we want to distinguish between different kind of "no data"
values. For example we may have an image of "Sea Surface Temperature" (SST)
measurements, with no data in some area because there is clouds, and no data in
other areas because we are on land. Some argue that we can not use Float.NaN for
"no data" value because they still want to distinguish the reason why the data
is missing (clouds, land, etc.).

This is not a well known fact, but actually there is many NaN values available
in IEEE 754. The Float.NaN constant is just one of them. We can build the other
ones with Float.intBitsToFloat(int) and a few knowledge of the bits pattern in a
IEEE 754 floating point number. So it is possible to represent "no data" values
by NaN and still be able to distinguish between "clouds, lands" or any other
category of your choice. The inconvenient is that we need to be careful with
portbility (we may have subtle difference in IEEE 754 interpretation between
different kind of processors).

If we use NaN values for "no data" instead of some special value like -9999,
then all JAI operations have a reasonable default behavior. In the case of the
convolve operator, if there is any NaN value under the convolution kernel, then
the result for that particular pixel is NaN. We don't get a totally irrelevant
value because of the mix of real data with special value like -9999; we get only
sensible result.

With the "NaN for nodata value" approach, there is no need for mask.

The inconvenient of the default implementation of the Convolve operator is that
if we use a kernel of said 4x4 and if there is only 1 NaN value among the 16
pixels under that kernel, then we get NaN. So we can loose a lot of area that
way. We could create a variant of the Convolve operator where NaN values are
ignored instead than summed.

Does it fit your use case?

	Martin

Re: [Geotools-gt2-users] masking a grid coverage

[Michael B. <mic...@gm...>]

2009/1/10 Martin Desruisseaux <mar...@ge...>:
>
> The inconvenient of the default implementation of the Convolve operator is that
> if we use a kernel of said 4x4 and if there is only 1 NaN value among the 16
> pixels under that kernel, then we get NaN. So we can loose a lot of area that
> way. We could create a variant of the Convolve operator where NaN values are
> ignored instead than summed.
>
> Does it fit your use case?
>
>        Martin
>

Yes, I work with some data that are (monitoring woodland clearing and
thinning) very similar to your examples in having various categories
of nodata: clouds, inapplicable areas, and areas with high uncertainty
ascribed to the classification.  So if I am attempting to calculate
local losses over time I have the similar issues to you with how to
control the convolution process.

Then there is another project where I am working on simulations of
animal populations where the summed value of habitat resources
available to an animal is constrained by the presence of other animal
territories.

In the present code for the modified convolve operator there are two
'mask' options: one which controls the positioning of the kernel and
the other which controls which source pixels are considered by the
kernel.  At the moment the operator only uses one controlling image
(an ROI object) for both options, but this would be easy to modify.
You could have the operator accept multiple ROIs, e.g. one to exclude
pixels where convolution isn't required / valid and a second to filter
source pixels that should not contribute to the kernel sum.  In the
latter case it might also be useful to have the option of generating
another output image to record the number or proportion of non-zero
kernel position that were available (unmasked) for each destination
pixel.

One can also imagine making the operator (or probably a second
operator) use a set of rules prescribing how to respond to source
image values directly rather than using separate controlling images.
I imagine this would result in slow processing - but it would avoid
the need to construct ROIs or similar so might be preferable in some
cases such as your example of differently valued NaNs in a source
image.

What do you think ?

Michael

Re: [Geotools-gt2-users] masking a grid coverage

[Martin D. <mar...@ge...>]

Hello Mickael

Michael Bedward a écrit :
> Then there is another project where I am working on simulations of
> animal populations where the summed value of habitat resources
> available to an animal is constrained by the presence of other animal
> territories.

This is a nice example, but I wonder: using a mask in this case means applying
the contraints in a "all or nothing" fashion. This is quite convenient for
exploratory work. But as the model become more elaborate, I would intuitively (I
may be totally wrong - I work more with oceans than habitats) expect the "all or
nothing" behavior to be replaced by some sigmoid function like the one used in
artificial neural network (http://en.wikipedia.org/wiki/Sigmoid_function). So if
R is the habitat resources, T the other animal territories and S some sigmoid
function, I would be tempted to compute the convolution (without mask) of:

  R * S(T)

Using a mask could be seen as a special case of the above formula where S is a
function going from 0 to 1 at a threshold value, without intermediate values.
The advantage of the above formula instead than a mask is to allow modelizer to
experiment different S functions having intermediate values, which is probably
closer to the real world.


> One can also imagine making the operator (or probably a second
> operator) use a set of rules prescribing how to respond to source
> image values directly rather than using separate controlling images.
> I imagine this would result in slow processing - but it would avoid
> the need to construct ROIs or similar so might be preferable in some
> cases such as your example of differently valued NaNs in a source
> image.

I don't have a clear opinion on this topic since I believe that it depends on
the experience gained with pratical work, and this experience is likely to vary
in different applications. The rules approach is nice if the rules can be simple
and "universal" enough, otherwise controlling images are preferable. But my
notion of "universal" is biased toward what I have seen in oceanography
(continuous functions are more frequent than discontinuities, hence my
preference for a sigmoid function in my previous paragraph), and my position
sometime change with time... So I can not really emmit a reliable advice here.

	Regards,

		Martin

Re: [Geotools-gt2-users] masking a grid coverage

[Michael B. <mic...@gm...>]

Hello Martin,
> This is a nice example, but I wonder: using a mask in this case means applying
> the contraints in a "all or nothing" fashion. This is quite convenient for
> exploratory work. But as the model become more elaborate, I would intuitively (I
> may be totally wrong - I work more with oceans than habitats) expect the "all or
> nothing" behavior to be replaced by some sigmoid function

Yes - I think you are right.  I have used various forms of sigmoid
functions, esp. inverse logistic, in other ecological models.  In my
present project the model is directed towards species that establish
and defend territories, so the crude representation of a bounded,
exclusive area is (possibly) not unreasonable - even though in the
real world adjacent territory boundaries are fuzzy and resource
availability there would be better modelled by a sigmoid function.
But then modelling is mostly about the trade-off between complexity
and detail of representation on the one hand against tractability,
testability and being able to actually parameterize the model on the
other hand.  I imagine this is also the case in oceanography.  More of
an art than a science :-)

> I don't have a clear opinion on this topic since I believe that it depends on
> the experience gained with pratical work, and this experience is likely to vary
> in different applications. The rules approach is nice if the rules can be simple
> and "universal" enough, otherwise controlling images are preferable. But my
> notion of "universal" is biased toward what I have seen in oceanography
> (continuous functions are more frequent than discontinuities, hence my
> preference for a sigmoid function in my previous paragraph), and my position
> sometime change with time... So I can not really emmit a reliable advice here.

Would it be possible / useful to add new operators to a new or
existing unsupported module ?  I need to do more such work for my
project anyway, and if I could add them to the geotools unsupported
modules it would be very helpful for me to get feedback on the
operators.

Michael

Re: [Geotools-gt2-users] masking a grid coverage

[Martin D. <mar...@ge...>]

Michael Bedward a écrit :
> Would it be possible / useful to add new operators to a new or
> existing unsupported module ?  I need to do more such work for my
> project anyway, and if I could add them to the geotools unsupported
> modules it would be very helpful for me to get feedback on the
> operators.

Yes we can add the operation in an unsupported module. Just for information, if
the operators are pure JAI, the following project may be worth a look. I have no
idea how good it is and if it still active however:

    https://jai-operators.dev.java.net/

If the operators involve interaction with GeoTools classes, then it is probably
better to keep them as an unsupported module in the geotools SVN.

	Martin

Re: [Geotools-gt2-users] masking a grid coverage

[Michael B. <mic...@gm...>]

2009/1/13 Martin Desruisseaux wrote:

> Yes we can add the operation in an unsupported module. Just for information, if
> the operators are pure JAI, the following project may be worth a look. I have no
> idea how good it is and if it still active however:
>
>    https://jai-operators.dev.java.net/

Thanks Martin.  There seems to be some recent activity there.  Yes, my
code is pure JAI so I'll ask about contributing it to jai-operators
and see what happens :)

Michael

Re: [Geotools-gt2-users] masking a grid coverage

[Michael B. <mic...@gm...>]

> 2009/1/13 Martin Desruisseaux wrote:
>
>> Yes we can add the operation in an unsupported module. Just for information, if
>> the operators are pure JAI, the following project may be worth a look. I have no
>> idea how good it is and if it still active however:
>>
>>    https://jai-operators.dev.java.net/
>

Hi Martin

I know you're not short of things to do or think about at the moment
:-)  But Ive been following up on the discussion that we had about jai
operators...

The jai-operators project seems to be hibernating at the moment.  The
owner is busy with other things.  So I've created a google code
project (Simone's suggestion) to host the new operators that I'm
writing plus some related work on a map expression interpreter (sort
of following the syntax of r.mapcalc in GRASS).

http://code.google.com/p/jai-tools/

There isn't any operator code there yet - just some preliminary bits
for the interpreter.  I'll upload the masked convolve operator that we
chatted about once I've tidied up the code.  Also Simone and others
might use the project to publish operators that they write.

Michael