[raster] Add manual height/width support to ST_Resample

[serveralex]

Either add support to specify these parameters manually or use these values directly from the reference raster (which makes the semantics far more clear). The "normal" user would expect that the ST_Resample result raster would have the exact same parameters as the reference raster, which isn't the case currently.

[Bborie Park]

The problem is that specifying the width/height is exclusive to that of the X and Y scales. So, for ST_Resample(raster, reference_raster), should the reference raster's scales be used or the dimensions? The assumption that you're making about the resulting raster having the exact same parameters as the reference raster is incorrect if the two rasters are different (such as where the extent of one is larger than that of the other).

Using the reference raster's scale makes sense for ST_Resample(raster, reference_raster) as the default instead of the width and height because the resampled raster has the same alignment with the reference raster. This permits followup processing using a two raster mapalgebra function.

I'll add the options of specify a resampled raster's width/height to the ST_Resample variants. ST_Resample for a raster and a reference raster may need an additional flag indicating whether or not to use scales or width/height.

So, instead of

CREATE OR REPLACE FUNCTION st_resample(
rast raster,
ref raster,
algorithm text DEFAULT 'NearestNeighbour',
maxerr double precision DEFAULT 0.125
)

add

CREATE OR REPLACE FUNCTION st_resample(
rast raster,
ref raster,
algorithm text DEFAULT 'NearestNeighbour',
usedimensions boolean DEFAULT FALSE,
maxerr double precision DEFAULT 0.125
)

[pracine]

Resampling a raster to the dimensions of another raster does not make much sense to me.

You might want to align a raster to another one in which case you need to borrow its scales, its skew and its upperleft corner to determine the new grid on which to align the new resampled raster. We offer the possibility to borrow them from another raster because they are generally cumbersome to write by hand. You may also want to explicitly pass those parameters.

You may want to explicitly pass a new width and a new height for the resampled raster but since a height and a width are generally easy to write (being integers) it seems useless to have the possibility to get them from another raster. This is also a dangerous operation to do on the many rasters of a single table as you might end with rasters having different pixel size. Always better to specify pixel size.

I would prefer the addition of a new signature accepting width and height than having one borrowing of dimensions from another raster.

[Bborie Park]

That's my current thought. Adding another signature is the safe approach and allows end-users to do what they want…

[bnordgren]

Replying to dustymugs:

The assumption that you're making about the resulting raster having the exact same parameters as the reference raster is incorrect if the two rasters are different (such as where the extent of one is larger than that of the other).

Is there a solid justification for this? What does it matter if the requested output grid is smaller or larger than the input? All you need to do is iterate over every cell in the requested grid and perform the resample operation. In any case, the output and input grid definitions will never be identical (why resample to the same grid? you already have it!)

Replying to pracine:

You might want to align a raster to another one in which case you need to borrow its scales, its skew and its upperleft corner to determine the new grid on which to align the new resampled raster. We offer the possibility to borrow them from another raster because they are generally cumbersome to write by hand. You may also want to explicitly pass those parameters.

I agree. The unit of borrowing should be the entire grid definition (width, height, and the affine transform to the geopoint). Note that this also requires that the target raster "borrow" the SRID, as the affine transform is tied to a particular geospatial CRS.

However, I don't see why it should be necessary to alter the signature in order to do this. The user can construct a raster with no bands, configure its grid definition to their heart's content, then require that the resampling operation respect the grid definition specified in the "ref" raster. If they want to resample many sources to a common grid, they just keep specifying the same "ref" raster over and over again.

I would prefer the addition of a new signature accepting width and height than having one borrowing of dimensions from another raster.

Again, width and height alone are not useful without all of the other contextual information provided in a grid definition. All or nothing!

[pracine]

Replying to bnordgren:

Replying to pracine:

You might want to align a raster to another one in which case you need to borrow its scales, its skew and its upperleft corner to determine the new grid on which to align the new resampled raster. We offer the possibility to borrow them from another raster because they are generally cumbersome to write by hand. You may also want to explicitly pass those parameters.

I agree. The unit of borrowing should be the entire grid definition (width, height, and the affine transform to the geopoint). Note that this also requires that the target raster "borrow" the SRID, as the affine transform is tied to a particular geospatial CRS.

I don't see the point of providing the width and the height at the same time as the affine transform… You are then stretching the raster to the footprint of the other one which is not a resampling operation but more like a wrapping operation. We decided (long ago) that wrapping would be a different function (ST_Wrap) yet to be designed…

Again the same argument about on-the-fly reprojection… We should not reproject on the fly to a new SRID before PostGIS do it itself. ST_Resample(raster, ref_raster) should fail with a warning if both SRID are different.

However, I don't see why it should be necessary to alter the signature in order to do this. The user can construct a raster with no bands, configure its grid definition to their heart's content, then require that the resampling operation respect the grid definition specified in the "ref" raster. If they want to resample many sources to a common grid, they just keep specifying the same "ref" raster over and over again.

Sure but this is tedious. We provide functions variants to make the life of our users easy so they choose our software!

[bnordgren]

Replying to pracine:

Replying to bnordgren:

I agree. The unit of borrowing should be the entire grid definition (width, height, and the affine transform to the geopoint). Note that this also requires that the target raster "borrow" the SRID, as the affine transform is tied to a particular geospatial CRS.

I don't see the point of providing the width and the height at the same time as the affine transform… You are then stretching the raster to the footprint of the other one which is not a resampling operation but more like a wrapping operation. We decided (long ago) that wrapping would be a different function (ST_Wrap) yet to be designed…

???!!!?? You lost me.

"Stretching" is performed only by the affine transform, not by the width/height values.

"Wrapping", I can't find anywhere on the trac site except for #454, which does not appear to be what you're talking about.

As to "the point" of providing width/height values at the same time as the affine transform: you're already requiring that the user provide width and height, as these are mandatory parameters of the raster type; you're just ignoring the values they provide in favor of other values which are computed by some unspecified means over which the user has no control. This caused user confusion that their provided values were ignored, hence this ticket. "The point" of respecting these values would be to produce the result requested by the user.

To put it another way, if the user specifically asked for a particular result, there should be a very good reason to give them something different. Is there a good reason?

Again the same argument about on-the-fly reprojection… We should not reproject on the fly to a new SRID before PostGIS do it itself. ST_Resample(raster, ref_raster) should fail with a warning if both SRID are different.

I was not making a statement about whether we should project on the fly or not. The affine transform is already tied to a particular SRID. If you blindly apply the transform without paying attention to which SRID it applies to, you'll successfully return when you should have raised an exception. Not an argument about on-the-fly projection at all.

However, I don't see why it should be necessary to alter the signature in order to do this. The user can construct a raster with no bands, configure its grid definition to their heart's content, then require that the resampling operation respect the grid definition specified in the "ref" raster. If they want to resample many sources to a common grid, they just keep specifying the same "ref" raster over and over again.

Sure but this is tedious. We provide functions variants to make the life of our users easy so they choose our software!

Let me get this straight: It's not tedious to ask them to supply a "reference raster" where the function ignores the required "height" and "width" parameters, but it is tedious to ask them to supply a "reference raster" when those parameters are respected? ???!?!!!!!

[serveralex]

To clearify the use case I'll explain why I need to set height/weight manually. Basically I have a very big raster table with square HxW raster blocks in WGS84 and one arcsec pixel scale that covers a big landscape. so basically i'm doing geostatistical interpolation of measuring samples with block kriging in realtime. that means i collect measuring samples and when a specific "new" limit is exceeded I'm interpolating the values for all raster blocks in that sector. The kriging has to be done in projected mode, so I transform the edge coordinates of the raster block to UTM projection and build a slightly bigger grid with square 25m pixels. Then I resample the resulting interpolated data to a WGS84 raster with the exact same extents as the original ones. But that won't work as the UTM raster extents are slightly bigger than the WGS84 ones which means the resulting raster won't be HxW in size any longer. So basically I have to set the height/weight parameters manually.

[pracine]

The width or the weight?

[Bborie Park]

Added support for width/height in r8153.

The function signature changes are:

_st_resample(
	rast raster,
	algorithm text DEFAULT 'NearestNeighbour', maxerr double precision DEFAULT 0.125,
	srid integer DEFAULT NULL,
	scalex double precision DEFAULT 0, scaley double precision DEFAULT 0,
	gridx double precision DEFAULT NULL, gridy double precision DEFAULT NULL,
	skewx double precision DEFAULT 0, skewy double precision DEFAULT 0,
	width integer DEFAULT NULL, height integer DEFAULT NULL
)

width: the number of columns in the raster. width/height are mutually exclusive to scalex/scaley

height: the number of rows in the raster. width/height are mutually exclusive to scalex/scaley

st_resample(
	rast raster,
	ref raster,
	algorithm text DEFAULT 'NearestNeighbour',
	maxerr double precision DEFAULT 0.125,
	usescale boolean DEFAULT TRUE
)

usescale: if FALSE, the reference raster's width and height are used instead of the X and Y scales.

New function signatures are:

st_resample(
	rast raster,
	width integer, height integer,
	srid integer DEFAULT NULL,
	gridx double precision DEFAULT NULL, gridy double precision DEFAULT NULL,
	skewx double precision DEFAULT 0, skewy double precision DEFAULT 0,
	algorithm text DEFAULT 'NearestNeighbour', maxerr double precision DEFAULT 0.125
)

st_resample(
	rast raster,
	ref raster,
	usescale boolean,
	algorithm text DEFAULT 'NearestNeighbour',
	maxerr double precision DEFAULT 0.125
)

[bnordgren]

Replying to dustymugs:

Added support for width/height in r8153.

The function signature changes are:

…

width: the number of columns in the raster. width/height are mutually exclusive to scalex/scaley

height: the number of rows in the raster. width/height are mutually exclusive to scalex/scaley

…

usescale: if FALSE, the reference raster's width and height are used instead of the X and Y scales.

Ok, I think what we need here is an explanation of how the un-specified/ignored parameters are calculated. On top of that, we need a justification for why these parameters are treated as mutually exclusive when they are not. In particular, why are you disallowing the perfectly legitimate case where the user provides a reference raster and wants you to return something which is exactly on the same grid.

Optionally calculating missing parameters is a whole different ball game than imposing an external criteria (relationship between width/height and scalex/scaley) the user may not want.

[Bborie Park]

As ST_Resample behaves the same way as GDAL gdalwarp utility, this latest change maintains that behavior. gdalwarp does not permit specifying scale and width/height at the same time.

​http://gdal.org/gdalwarp.html

Look at switches -tr and -ts.

You can only have one or the other as they're dependent upon one another. Using GDAL's suggested extent for the warped raster, you can compute the dependent variable. If width/height is provided, the scales are computed for the suggested extent. If scale is provided, the width and height are also computed.

[bnordgren]

Using gdalwarp as a convenience in the short term is a good thing to get us up and running quickly. But it really doesn't add the functionality requested in this ticket.

height and width are not at all dependent on scalex and scaley unless you're imposing some external constraint (gdalwarp's suggested size). In this case, the external constraint computed the wrong height and width for severalex, so it will also compute the wrong scalex/scaley when provided the correct height and width. It is the external constraint which is incorrect for this application.

I'm reopening this for postgis future, because it's not solved. We can close it as "wontfix" if we never plan on adding this functionality, or we can leave it open to remind future us of something that would be nice. Maybe we can even change the title to something reflecting the need to eliminate the assumptions imposed on the operation by gdalwarp. But we can't close it as "fixed."

[robe]

Milestone renamed