[raster] ST_Resample

[dustymugs]

With the use of GDAL's warp API, a full set of different functions is needed to provide end-user functions. Therefore, the following function underlies all the other functions:

_ST_Resample(
	rast raster,
	srid integer DEFAULT NULL,
	scalex double precision DEFAULT 0,
	scaley double precision DEFAULT 0,
	upperleftx double precision DEFAULT NULL,
	upperlefty double precision DEFAULT NULL,
	skewx double precision DEFAULT NULL,
	skewy double precision DEFAULT NULL,
	algorithm text DEFAULT 'NearestNeighbour',
	maxerr double precision DEFAULT 0.125
) -> raster

User-facing functions will provide needed validation before passing parameters to the above function.

As can be seen above, there is a "srid" parameter. The expectation is that ST_Transform will be refactored to call _ST_Resample.

The most comprehensive user function using _ST_Resample is:

ST_Resample(
	rast raster,
	srid integer DEFAULT NULL,
	scalex double precision DEFAULT 0,
	scaley double precision DEFAULT 0,
	upperleftx double precision DEFAULT NULL,
	upperlefty double precision DEFAULT NULL,
	skewx double precision DEFAULT NULL,
	skewy double precision DEFAULT NULL,
	algorithm text DEFAULT 'NearestNeighbour',
	maxerr double precision DEFAULT 0.125
) -> raster

Parameters passed to ST_Resample will be validated to ensure that at least one of the srid, scale, upperleft and skew parameters are provided.

returns the resampled raster

srid: the SRID of the projection to use when reprojecting the raster

scalex: the resampled raster's scale in the X axis. default is 0 indicating that the user isn't specifying the resampled raster's scale

scaley: the resampled raster's scale in the Y axis. default is 0 indicating that the user isn't specifying the resampled raster's scale

upperleftx: the X component of the corner of the upper-left pixel of the resampled raster. default is NULL indicating that the user isn't realigning the raster

upperlefty: the Y component of the corner of the upper-left pixel of the resampled raster. default is NULL indicating that the user isn't realigning the raster

skewx: the X component of the resampled raster's skew. default is NULL indicating that the skew of the original raster will be applied to the the resampled raster

skewy: the Y component of the resampled raster's skew. default is NULL indicating that the skew of the original raster will be applied to the the resampled raster

algorithm: the resampling algorithm to use when resampling the raster. default is 'NearestNeighbour'. possible algorithms are:

NearestNeighbour (default.  fastest performance but worst interpolation)

NearestNeighbor (for those wanting to use the American spelling)

Bilinear

Cubic

CubicSpline

Lanczos

maxerr: the threshold for transformation approximation by the resampling algorithm (in pixel units). default is 0.125, which is the same value used in GDAL gdalwarp utility. if set to zero, no approximation takes place.

The other variant of ST_Resample is:

ST_Resample(
	rast raster,
	ref raster,
	algorithm text DEFAULT 'NearestNeighbour',
	maxerr double precision DEFAULT 0.125
) -> raster

ref: the reference raster whose SRID, scale, upperleft coordinate and skew will be applied to the "rast" raster

Additional functions will be documented in following entries for this ticket.

[dustymugs]

For setting the raster's scale only:

1. ST_Rescale(
	rast raster,
	scalex double precision,
	scaley double precision,
	algorithm text DEFAULT 'NearestNeighbour',
	maxerr double precision DEFAULT 0.125
) -> raster

scalex: the resampled raster's scale in the X axis. default is 0 indicating that the user isn't specifying the resampled raster's scale

scaley: the resampled raster's scale in the Y axis. default is 0 indicating that the user isn't specifying the resampled raster's scale

2. ST_Rescale(
	rast raster,
	scalexy double precision,
	algorithm text DEFAULT 'NearestNeighbour',
	maxerr double precision DEFAULT 0.125
) -> raster

scalexy: the reprojected raster's scale in the X and Y axes. default is 0 indicating that the user isn't specifying the reprojected raster's scale

[dustymugs]

For changing the alignment of a raster:

1. ST_Realign(
	rast raster,
	upperleftx double precision,
	upperlefty double precision,
	skewx double precision DEFAULT NULL,
	skewy double precision DEFAULT NULL,
	algorithm text DEFAULT 'NearestNeighbour',
	maxerr double precision DEFAULT 0.125
) -> raster

upperleftx: the X component of the corner of the upper-left pixel of the resampled raster. default is NULL indicating that the user isn't realigning the raster

upperlefty: the Y component of the corner of the upper-left pixel of the resampled raster. default is NULL indicating that the user isn't realigning the raster

skewx: the X component of the resampled raster's skew. default is NULL indicating that the skew of the original raster will be applied to the the resampled raster

skewy: the Y component of the resampled raster's skew. default is NULL indicating that the skew of the original raster will be applied to the the resampled raster

2. ST_Realign(
	rast raster,
	upperleftx double precision,
	upperlefty double precision,
	algorithm text,
	maxerr double precision DEFAULT 0.125
) -> raster

[dustymugs]

For changing the skew of a raster:

1. ST_Reskew(
	rast raster,
	skewx double precision,
	skewy double precision,
	algorithm text DEFAULT 'NearestNeighbour',
	maxerr double precision DEFAULT 0.125
)

skewx: the X component of the resampled raster's skew. default is NULL indicating that the skew of the original raster will be applied to the the resampled raster

skewy: the Y component of the resampled raster's skew. default is NULL indicating that the skew of the original raster will be applied to the the resampled raster

[pracine]

What is the advantage of having four names (ST_Reskew, ST_Realign, ST_Rescale, ST_Resample) for what are all different variants of resample?

Two disadvantages I see:

- We have to maintain four different pages of documentation having a lot in common. - I'm scared that people would go directly to ST_Resample without seeing/using the other functions.

robe might have an opinion on that?

[dustymugs]

I think having four (technically five including ST_Transform as all these functions run against rt_raster_gdal_warp at the lowest level) works better due to cleaner function declaration and purpose. Having five different functions especially helps with being able to specify default parameters.

The disadvantage with having only ST_Resample (not counting ST_Transform due to it being a parallel to the geometry version) is that there is no easy way to maintain a consistent order for the function parameters if you want to make variations for just specific elements and still be able to set default parameter values.

I'm of the opinion that people will preferentially use the shortcut variations instead of the full ST_Resample. Say for example, the user wants to change the upper left corner:

ST_Resample(
	rast,
	NULL,
	0,
	0,
	-500000,
	400000
)

The above call is more involved than calling:

ST_Realign(
	rast,
	-500000,
	400000
)

I would expect that if a user needs to do multiple resampling operations that isn't addressed in a single variation (such as variation #1 of ST_Realign) and runs as efficiently as possible, they'd use ST_Resample.

[robe]

I agree with Bborie on this. To me this is similar in concept to ST_Relate.

Essentially (well as least in concept), you can do all the st_intersects, sT_Within, ST_crosses etc. with ST_Relate and &&, but(of course we do have speed enhancements and so forth behind the others which makes it more efficient to use those than using ST_Relate.

Then again you may find in future that you can optimize certain kinds of operations if you know that is the intent of people to do that specific thing. Those enhancements would become more obvious if they are separate functions.

I also hate having functions with the same name with argument orders swapped. Drives me a bit nuts, but that's just me.

[pracine]

ok

[dustymugs]

Pierre,

I've a question regarding how realigning a raster using the upper left pixel coordinate should happen.

From the image above, when a raster is being resampled with a new upper left corner, should it be like the upper or lower change? The red outline indicates the extent of the resampled raster.

I've written the resampling code to behave like the upper where shifting the corner preserves the dimensions of the raster.

Thoughts?

[pracine]

It is not explicit in the spec but you never want to do that... You would lose most of your raster. What you might want to do however it to realign a raster so it matches the grid of a second one. To do that you pass the upper left corner x and y of the second to ST_Resample and you just snap the upper left corner of the first to the grid defined by this x and y and the provided pixelsize.

In other words x and y are not really the upperleft corner of the second raster but just the upper left corner of any pixel of the grid you want to snap to.

[dustymugs]

I agree that the upper example can be detrimental to your raster. I'll tweak the code to behave like the lower example, which would preserve the raster data. If I'm incorrect about this, could you post a visual aid?

[dustymugs]

Have you had the chance to read my question regarding ST_Transform (and thus ST_Resample) in the devel mailing list? If not, I'll post the content here...

From what I can tell, the GDAL warp API underlying ST_Resample automatically de-skews the input raster. Is this what we want? Or do we want the input raster's skew maintained in the resampled raster?

[pracine]

Replying to dustymugs:

I agree that the upper example can be detrimental to your raster. I'll tweak the code to behave like the lower example, which would preserve the raster data. If I'm incorrect about this, could you post a visual aid?

I don't think so. That would unwisely change the shape of the raster footprint. I knew I described this better. It is in the ST_AsRaster() specs where I describe how to align the desired raster. I will try to describe it better here...

In the ST_Resample specs, the second series of variant to not define what should be the new upperleft corner of the raster. They define a grid on which the raster should snap to. So it is neither of your drawings... An alignment grid is defined by x, y and a pixelsize. The pixel size can be defined with just 'pixelsize', meaning that scalex = scaley and skewx = skewy = 0 of with some of those pixelsize parameters (variants 7 and 9). This is why I did not call them ulx and uly but just x and y. Those x and y are arbitrary upperleft corner of ANY pixel of the new grid to snap to. You have to compute the new ulx and uly of the resampled raster according/to snap to this grid.

Sorry if this was not explained...

[pracine]

Replying to dustymugs:

Have you had the chance to read my question regarding ST_Transform (and thus ST_Resample) in the devel mailing list? If not, I'll post the content here... From what I can tell, the GDAL warp API underlying ST_Resample automatically de-skews the input raster. Is this what we want? Or do we want the input raster's skew maintained in the resampled raster?

Shoul I understand here that it is not possible to resample to a rotated grid using GDAL?

I agree that it is not clear what should be the desired behaviour with the first series of variant when the raster is rotated... If the raster is not rotated only changing the pixel size make sense (variant 1 and 3). You want the upperleft corner to stay the same. But what if the raster is rotated? I would solve this like this:

1) Take for granted that when no skewx and skewy are provided they are 0 and hence we lose the rotation if there was one.

2) Add a variant to the first series taking skewx ans skewy (like in the second series) allowing to specify the new rotation if we want to keep one.

Do you want me to update the specs to reflect that?

[dustymugs]

Replying to pracine:

I don't think so. That would unwisely change the shape of the raster footprint. I knew I described this better. It is in the ST_AsRaster() specs where I describe how to align the desired raster. I will try to describe it better here... In the ST_Resample specs, the second series of variant to not define what should be the new upperleft corner of the raster. They define a grid on which the raster should snap to. So it is neither of your drawings... An alignment grid is defined by x, y and a pixelsize. The pixel size can be defined with just 'pixelsize', meaning that scalex = scaley and skewx = skewy = 0 of with some of those pixelsize parameters (variants 7 and 9). This is why I did not call them ulx and uly but just x and y. Those x and y are arbitrary upperleft corner of ANY pixel of the new grid to snap to. You have to compute the new ulx and uly of the resampled raster according/to snap to this grid. Sorry if this was not explained...

I think I understand what you're describing. I'm going to have to let this sink in for a bit to understand the calculation involved.

[dustymugs]

Replying to pracine: \> Shoul I understand here that it is not possible to resample to a rotated grid using GDAL?

I agree that it is not clear what should be the desired behaviour with the first series of variant when the raster is rotated... If the raster is not rotated only changing the pixel size make sense (variant 1 and 3). You want the upperleft corner to stay the same. But what if the raster is rotated? I would solve this like this: 1) Take for granted that when no skewx and skewy are provided they are 0 and hence we lose the rotation if there was one. 2) Add a variant to the first series taking skewx ans skewy (like in the second series) allowing to specify the new rotation if we want to keep one. Do you want me to update the specs to reflect that?

No. You can resample to a rotated grid as you can specify the output raster's geotransform, which gets used in the third stage of the GDALCreateGenImgProjTransformer2 function.

 http://gdal.org/gdal__alg_8h.html#94cd172f78dbc41d6f407d662914f2e3

I think the desired behavior is to maintain the skew of the input raster. Only by specifying skew parameters should the resampled raster's skew change.

Is it possible for you (or anyone else for that matter) to send me a small set of skewed rasters? I could use landsat imagery but I'm looking for something more manageable.

[pracine]

Replying to dustymugs:

Replying to pracine:

I don't think so. That would unwisely change the shape of the raster footprint. I knew I described this better. It is in the ST_AsRaster() specs where I describe how to align the desired raster. I will try to describe it better here... In the ST_Resample specs, the second series of variant to not define what should be the new upperleft corner of the raster. They define a grid on which the raster should snap to. So it is neither of your drawings... An alignment grid is defined by x, y and a pixelsize. The pixel size can be defined with just 'pixelsize', meaning that scalex = scaley and skewx = skewy = 0 of with some of those pixelsize parameters (variants 7 and 9). This is why I did not call them ulx and uly but just x and y. Those x and y are arbitrary upperleft corner of ANY pixel of the new grid to snap to. You have to compute the new ulx and uly of the resampled raster according/to snap to this grid. Sorry if this was not explained...

I think I understand what you're describing. I'm going to have to let this sink in for a bit to understand the calculation involved.

This is basically the same behaviour as ST_SnapToGrid... We should even provide a set of alias to the second series of ST_Resample variants named SnapToGrid?. The first set of variant do not really 'snap to grid'. They only change the pixel size without touching the upperleft corner.

[pracine]

Replying to dustymugs:

Replying to pracine: \> Shoul I understand here that it is not possible to resample to a rotated grid using GDAL?

I agree that it is not clear what should be the desired behaviour with the first series of variant when the raster is rotated... If the raster is not rotated only changing the pixel size make sense (variant 1 and 3). You want the upperleft corner to stay the same. But what if the raster is rotated? I would solve this like this: 1) Take for granted that when no skewx and skewy are provided they are 0 and hence we lose the rotation if there was one. 2) Add a variant to the first series taking skewx ans skewy (like in the second series) allowing to specify the new rotation if we want to keep one. Do you want me to update the specs to reflect that?

No. You can resample to a rotated grid as you can specify the output raster's geotransform, which gets used in the third stage of the GDALCreateGenImgProjTransformer2 function.  http://gdal.org/gdal__alg_8h.html#94cd172f78dbc41d6f407d662914f2e3 I think the desired behavior is to maintain the skew of the input raster. Only by specifying skew parameters should the resampled raster's skew change.

Actually in all the ST_Resample signatures, pixelsize, pixelsizex and pixelsizey should be renamed scale, scalex and scaley. For me if you don't specify skewx and skewy it mean you want them to be set to 0.

Is it possible for you (or anyone else for that matter) to send me a small set of skewed rasters? I could use landsat imagery but I'm looking for something more manageable.

Can't you just do:

CREATE TABLE skewedrasttable AS SELECT ST_SetSkew(rast, x, y) FROM yourtable

Much easier than in ArcGIS...

[dustymugs]

I'll implement what we've discussed for both snapping to a grid and only setting a skew when specified by a user, thus defaulting to de-skewing a raster.

I believe no parameters I proposed are in the "pixelsize" style, just scale.

[dustymugs]

The functions ST_Resample, ST_Rescale, ST_Reskew and ST_SnapToGrid have been committed in r7633. ST_Transform now points to _ST_Resample which underlies all functions of the ST_Resample family.

Only nagging concerns in my head are ST_Reskew (whether or not the raster outputted is correct) and ST_SnapToGrid (whether or not the lower right corner of the aligned raster contains the lower right corner of the original raster). In my testing and the regression test, everything looks right but only real testing will tell.

[pracine]

I get this when making check:

--- rt_resample_expected	2011-07-21 10:14:23 -0400
+++ /tmp/pgis_reg_5796/test_41_out	2011-07-21 10:21:48 -0400
@@ -8,7 +8,7 @@
 1.12|992163|10|10|1|1000.000|-1000.000|0.000|0.000|-500001.000|600000.000|t|t|t
 1.13|992163|10|10|1|1000.000|-1000.000|0.000|0.000|-500000.000|600009.000|t|t|t
 1.14|992163|10|11|1|1000.000|-1000.000|0.000|0.000|-500100.000|600950.000|t|t|t
-1.15|992163|200|200|1|50.000|-50.000|0.000|0.000|-500010.000|600007.000|t|t|t
+1.15|992163|201|200|1|50.000|-50.000|0.000|0.000|-500040.000|600007.000|t|t|t
 1.16|992163|83|83|1|121.000|-121.000|0.000|0.000|-500093.000|600039.000|t|t|t
 1.17|993310|243|243|1|50.000|-50.000|0.000|0.000|950710.000|1409307.000|t|t|t
 1.18|993309|242|243|1|50.000|-50.000|0.000|0.000|950760.000|1409107.000|t|t|t
@@ -72,7 +72,7 @@
 5.24|992163|83|83|1|121.000|-121.000|0.000|0.000|-500000.000|600040.000|t|t|t
 5.25|992163|83|83|1|121.000|-121.000|0.000|0.000|-500000.000|600048.000|t|t|t
 5.26|992163|83|83|1|121.000|-121.000|0.000|0.000|-500098.000|600040.000|t|t|t
-5.27|992163|83|83|1|121.000|-121.000|0.000|0.000|-500037.000|600030.000|t|t|t
+5.27|992163|83|83|1|121.000|-121.000|0.000|0.000|-500084.000|600030.000|t|t|t
 5.3|992163|10|10|1|1000.000|-1000.000|0.000|0.000|-500001.000|600000.000|t|t|t
 5.4|992163|10|10|1|1000.000|-1000.000|0.000|0.000|-500001.000|600000.000|t|t|t
 5.5|992163|11|10|1|1000.000|-1000.000|0.000|0.000|-500999.000|600000.000|t|t|t

[dustymugs]

Fixed in r7657. Was a floating point issue causing my results to be incorrect while your results are correct.