[raster] There should be an option for setting the nodata value in certain variants of ST_Intersection(raster, raster)

[pracine]

All the variants of ST_Intersection accepting a extenttype argument should also allow to define a nodata value in case the input raster do not have one defined.

Otherwise it leads to strange results like this:

SELECT ST_AsBinary((gv).geom) geom,  (gv).val
FROM (SELECT ST_PixelAsPolygons(
           ST_Intersection(ST_AddBand(ST_MakeEmptyRaster(48, 63, 0, 0, 0.001, -0.001, 0, 0, 4269), '32BF'::text, 0, NULL), 
                           ST_AddBand(ST_MakeEmptyRaster(48, 63, 0.1, 0.1, 0.001, -0.001, 0, 0, 4269), '32BF'::text, 0, NULL)::geometry,
                           'SECOND')) gv) foo;

in which:

-the first raster do not have a nodata value defined,

-the second parameter is just a dummy geometry to 'cut' the first raster,

-we want the result to cover the extent of the geometry ('SECOND')

so that we should get a lower left rectangle filled with zeros (0) and the rest of the extent filled with nodata.

Instead we get just nodata

If ST_Intersection(raster, geometry, extenttype) would allow me to pass a nodata value, then I could make sure a correct nodata value is chosen and hence not annihilating my data.

On the same track, why do the variants of ST_Intersection() taking a geometry allow me to define the resulting extent but not the variant taking two rasters? Actually why does the one taking a geometry allows me to define the resulting extent at all? This is an intersection operation, the result should always cover just the intersection. If I want something more sophisticated I should use ST_MapAlgebra(raster, raster).

I'm working on a ST_Tile() plpgsql prototype using ST_Clip() and ST_Intersection(raster, geometry)…

[pracine]

So I see just two options to solve this:

1) We add a nodata value parameter and this nodata value is used just if the requested band do not have a nodata value defined (if it does the passed nodata value is ignored). If this option is chosen I would also expect to be able to give an extent to the raster/raster variants.

2) We remove the possibility to get an other extent then the 'INTERSECTION'.

That said I like the fact that we have (in some cases) the possibility to get an extent other than 'INTERSECTION' taking for granted that only the intersection is always filled with data…

[pracine]

Another example involving the two raster version of ST_Intersection:

SELECT ST_BandNodataValue(ST_Intersection(ST_AddBand(ST_MakeEmptyRaster(48, 63, 0, 0, 0.001, -0.001, 0, 0, 4269), '8BSI'::text, 0, NULL), 
                                   ST_AddBand(ST_MakeEmptyRaster(48, 63, 0.01, -0.01, 0.001, -0.001, 0, 0, 4269), '8BSI'::text, 1, NULL), 'FIRST'));

the result of this is -128 meaning that the nodata value chosen must be equal to ST_MinPossibleValue('8BSI'). I have no way to change this other than by assigning a nodata value to the first raster before calling ST_Intersection(). Might be fair enough though…

[pracine]

I think a general rule should be: As soon as a function CAN produce nodata values (this is the case of ST_Intersection(raster, raster), ST_Intersection(raster, geometry) and ST_Clip(raster, geometry)) there should be a way to specify the selected nodata value in case the input raster do not have one defined.

Otherwise we have to prepend ST_SetBandNodataValue() every time (if we do not want strange results in some cases).

[dustymugs]

I wouldn't mind getting rid of the "extenttype" parameter from the ST_Intersection(raster, geometry) variants. Like you said, ST_Intersections should ONLY be returning the intersection.

As for being able to specify a NODATA, you should either use ST_SetBandNoDataValue before calling ST_Intersection or use 2-raster ST_MapAlgebra.

[pracine]

Actually setting the nodata value with ST_SetBandNoDataValue on the first raster do not solve the first problem:

SELECT ST_BandNodataValue(ST_Intersection(ST_AddBand(ST_MakeEmptyRaster(48, 63, 0, 0, 0.001, -0.001, 0, 0, 4269), '8BSI'::text, 0, 2), 
                           ST_AddBand(ST_MakeEmptyRaster(48, 63, 0.1, 0.1, 0.001, -0.001, 0, 0, 4269), '8BSI'::text, 0, NULL)::geometry,
                           'SECOND'))

do not set 2 as the nodata value… It select the default nodata value assigned to the rasterization of the geometry which is 0.

so this seems like a real bug to me. Might be related to #1537 though.

[dustymugs]

The bug in #1537 should now be resolved. I'm leaning towards removing the "extenttype" parameter for ST_Intersection(raster, geometry) as it does cause confusion.

[pracine]

#1537 fixed…

Even if you restraint the extent to the intersection you must do something for when the provided raster do not have a nodata value defined. If I intersect a raster having nodata value with a circle, for example, I want to have some control over the nodata value assigned to the resulting raster. Doing ST_SetBandNoDataValue() afterward is not sufficient as all the nodata values, in this case, have been set to 0. I have no other way than calling ST_MapAlgebraExpr another time to convert all of them to the nodata value I wish.

As for ST_Clip() ST_Intersection(raster, geometry) should let us provide a nodata value.

The same logic applies to ST_Intersection(raster, raster). What if I want the two bands, the second raster has -10 as nodata value, many of then in the intersecting area and the first raster has no nodata value defined (and worse many pixel having the value -10)? What should be the nodata value assigned to the first raster?

I think we should have a parameter for each raster to define the resulting nodata value. If I don't provide it then ST_MinPossibleValue() should be used (at the risk of transforming some withdata values into nodata values).

I think ST_MapAlgebra provides all the tools for ST_Intersection to be able to define or redefine the nodata value. You have to set those values using the alternate expressions and then redefine the nodata value with ST_SetBandNodataValue().

You might say then there should also be a nodata value paramater to ST_MapAlgebraExpr but no: such a parameter could be conflictual with the value set with the alternate expressions. It is the role of the caller to set the nodata value to the same defined with the alternate expressions.

What do you think?

[dustymugs]

I don't believe you need to have a nodata parameter as the intersection isn't doing anything to the intersecting rasters' bands unless you use 'OTHER' for the "returnband" parameter.

returnband: can be FIRST, SECOND, BOTH, OTHER

FIRST returns the band of rast1 in the intersecting extent. returning raster will have one band. SECOND returns the band of rast2 in the intersecting extent. returning raster will have one band. BOTH returns the bands of rast1 and rast2 in the intersection extent. returning raster will have two bands. OTHER returns the computed band based upon rast1 and rast2 in the intersecting extent. returning raster will have one band. If OTHER, must provide a regprocedure to otherfunc

So, specifying a NODATA isn't useful as that attribute would be copied over from the source band's metadata. If the source band's metadata doesn't have a NODATA, the output band shouldn't have a NODATA.

I did just notice that the user doc for ST_Intersection is not complete. The section in the Working Specs is complete and correct though.

[pracine]

I would expect, if I intersect a raster representing a circle surrounded with nodata values with another raster, that the area surrounding the circle would be filled with nodata values in BOTH bands…

Otherwise we can not say that ST_Intersection is taking nodata value into account. It is just an extent intersection.

Actually I guess this is what is happening based on the way you are calling ST_MapAlgebraExpr: since you are not providing alternate expressions, any time one pixel from any band is nodata both resulting bands gets nodata. This should be the expected behavior of ST_MapAlgebra. I check…

[dustymugs]

If I'm intersecting two rasters using ST_Intersection, I'm expecting that the returning raster would have the pixels of one or both bands in the intersecting area.

In your example with a circle raster and another raster with the function returning the band of the other raster, I'd expect the band returned to be constrained spatially to just the intersecting area and the pixels outside the circle to have a NODATA value.

In the situation where the other raster doesn't have a NODATA, it would be nice to have a NODATA parameter.

[pracine]

Replying to dustymugs:

If I'm intersecting two rasters using ST_Intersection, I'm expecting that the returning raster would have the pixels of one or both bands in the intersecting area.

Agree. The question is not about which band to return. It's about taking nodata value into account or not. If I intersect a nodata pixel with a withdata pixel, I expect to get nodata. Intersect means "there is a value in the first raster AND a value in the second raster". It is both an extent intersection AND a value intersection. Actually the nodata values define the extent the rasters. They are implicit masks.

That's the logic behind ST_Intersects(raster, geometry): It is not just an extent intersection. It is also a value intersection. If the geometry intersects only with nodata pixel, ST_Intersects returns false.

In your example with a circle raster and another raster with the function returning the band of the other raster, I'd expect the band returned to be constrained spatially to just the intersecting area and the pixels outside the circle to have a NODATA value.

In the situation where the other raster doesn't have a NODATA, it would be nice to have a NODATA parameter.

Correct.

[dustymugs]

So, I'll add a NODATA parameter for the ST_Intersection(raster, geometry) variants. If NODATA parameter is specified, it should supersede whatever NODATA value is specified for the raster band.

Adding a NODATA parameter for ST_Intersection(raster, raster) is harder as we'd need two parameters as the two bands used could have different pixel types.

Agreed?

[dustymugs]

Just to be certain before I make the changes.

Yes to adding NODATA parameter for the ST_Intersection(raster, geometry) variants.

No to adding NODATA parameters for the ST_Intersection(raster, raster) variants.

[pracine]

Replying to dustymugs:

Just to be certain before I make the changes.

Yes to adding NODATA parameter for the ST_Intersection(raster, geometry) variants.

Yes

No to adding NODATA parameters for the ST_Intersection(raster, raster) variants.

No, but we will have to do it later. In 2.1? I would then leave the ticket open but postpone it.

You agree with the logic?

[dustymugs]

I'll add NODATA parameters to both ST_Intersection(raster, raster) and ST_Intersection(raster, geometry) today. There's no point in not doing so when I'll be tweaking the code.

[pracine]

Also: At first I was thinking that those nodata value parameters should be used only when the requested band do not have a nodata value defined but I think this will become confusing.

I think in the absence of such an argument, if the raster do not have a nodata value defined the ST_MinPossibleValue() should be used.

If the argument is provided and the provided raster DO have a nodata value defined then the provided value should replace the nodata value of the raster (actually it does not replace it as we have to produce a new raster anyway).

You can achieve this by passing the provided nodata value as the alternate expressions. Here is an example I use in ST_Tile():

ST_SetBandNodataValue(ST_MapAlgebraExpr(rast, 1, geomrast, 1, 'RAST1', newpixtype, 'SECOND', newnodata::text, newnodata::text, newnodata), newnodata);

As I said previously it is the responsability of the caller of ST_MapAlgebraExpr() to determine and set the new nodata value and the new pixel type.

[dustymugs]

So, only one NODATA parameter for ST_Intersection(raster, raster)?

[pracine]

Replying to dustymugs:

So, only one NODATA parameter for ST_Intersection(raster, raster)?

No. Two. Actually it depends on what band are requested. If 'BOTH' band are requested and no raster has a nodata value then we need two nodata values. I think the parameter should be an array so it can hold on or two values.

[dustymugs]

So for various returnband:

FIRST

ST_MapAlgebraExpr(rast1, band1, rast2, band2, '[rast1.val]', ST_BandPixelType(rast1, band1), extenttype, nodata2, nodata1, nodata1);

SECOND

ST_MapAlgebraExpr(rast1, band1, rast2, band2, '[rast2.val]', ST_BandPixelType(rast2, band2), extenttype, nodata2, nodata1, nodata2);

BOTH is FIRST and then SECOND

[pracine]

I would say:

FIRST

ST_MapAlgebraExpr(rast1, band1, rast2, band2, '[rast1.val]', ST_BandPixelType(rast1, band1), 'INTERSECTION', nodata1, nodata1, nodata1); SECOND

ST_MapAlgebraExpr(rast1, band1, rast2, band2, '[rast2.val]', ST_BandPixelType(rast2, band2), 'INTERSECTION', nodata2, nodata2, nodata2);

BOTH is FIRST and then SECOND

The extent should always be 'INTERSECTION' unless we want to support more possibilities and nodata1 (not nodata2) should replace the existing nodata for 'FIRST'.

[dustymugs]

So you also want to get rid of the extenttype parameter?

Do you want to make the changes? At this point I'm torn between confusion and frustration.

[pracine]

So here are the list of changes I'm planning:

-Removal of the extenttype parameter. ST_Intersection should always return the geometrical intersection of both raster. If you want something else use ST_MapAlgebra() or the future ST_SymDifference() or ST_Difference().

-A new set of ST_Intersection() variant accepting an array of nodata values. The number of values in the array should match the number of band requested (FIRST and SECOND: 1, BOTH: 2). I'm worry about the third case: OTHER… What should we do in this case? What is the need for this option anyway?

-When one of the input raster corresponding to a requested band do not have a nodata value defined but the other one does AND no nodata value array is provided, we will set the nodata value of the one missing it to ST_MinPossibleValue() with a warning similar to the one planned in ST_Clip() (see #1581): "The FIRST (or SECOND) band did not have a nodata value defined and no nodata value were provided. Setting the FIRST (or SECOND) band nodata value to XXX. Use a variant accepting a nodata array to define it yourself." Very long notice…

[pracine]

extenttype parameter removed in r9229.

[pracine]

I redesigned the variants if I add a double precision[] nodataval parameter. I end up with: (= means this parameter has a default value. I shortened everything for ease of reading)

— raster/geometry
rast, geom, reg
rast, geom, nodata=, reg=
rast, band, geom, reg
rast, band, geom, nodata=, reg=

— raster/raster no band number
ras1, rast2, return=, nodaval[]=, reg=
ras1, rast2, return=, nodaval=, reg=
ras1, rast2, nodaval[], reg= (opt)
ras1, rast2, nodaval, reg= (opt)
ras1, rast2, reg (opt)

— raster/raster with band number
ras1, b1, rast2, b2, return=, nodaval[]=, reg=
ras1, b1, rast2, b2, return=, nodaval=, reg=
ras1, b1, rast2, b2, nodaval[], reg= (opt)
ras1, b1, rast2, b2, nodaval, reg= (opt)
ras1, b1, rast2, b2, reg (opt)

Maybe we can simplify by removing the three last one (the opt ones) of the two last groups and force users to pass something or a casted NULL in one of the two first one?

Could we, in the doc, when we have parameter accepting a single value or an array, like here the new nodataval parameter, not list the variants accepting an array and simply say in the doc "Every of those variant also accept an array of value"?

[pracine]

Fixed in r9336.