[raster] Physically relevant geotransform settings

[bnordgren]

I only have the forward calculations worked out and described on the wiki. I have an approach in mind for the reverse calculations.

Algorithm documentation: DevWikiAffineParameters

This could stand to be reviewed/edited (for clarity and correctness) by someone other than me. Code will follow feedback on the proposed algorithm.

This approach allows users to specify the geotransform in terms of the following parameters:

• The magnitude of the i_b basis vector (distance between pixels along i axis)
• The magnitude of the j_b basis vector (distance between pixels along j axis)
• The angle θ_i (the angle by which the raster's grid needs to be rotated, positive clockwise -- for compatibility with "heading/bearing")
• The angle θ_ij (the angle from the i basis vector to the j basis vector positive counterclockwise -- for consistency with right-handed coordinate system)
• The offsets t_x and t_y

This would likely result in a single function with a signature like:

ST_SetGeoTransform(raster raster, 
                   float8 i_magnitude,
                   float8 j_magnitude, 
                   float8 theta_i default 0,
                   float8 theta_ij default -90)

... or the equivalent if I pooched the SQL syntax.

We should also downplay the scale/skew accessors after this is implemented. They're still relevant for copying individual parameters from an external source, but it's almost never appropriate to change just one or two of them. ...and the names for the parameters are somewhat misleading.

[bnordgren]

The method to calculate these physically relevant parameters given an arbitrary transform is presented on DevWikiRealParameters.

This would likely result in a single function with a signature like:

ST_GetGeoTransform(raster rast)

returning a setof 4 values:

• i_magnitude
• j_magnitude
• theta_i
• theta_ij

After implementation, ST_GetRotation should be modified to return theta_i from this function. Currently, ST_GetRotation returns NaN if theta_ij is not +-90 degrees.

If other accessors for individual parameters are desired, they may be implemented on top of the two general functions presented here with little or no loss of performance.

[bnordgren]

Tentative implementation here:  https://github.com/bnordgren/postgis/blob/ri-gen2-svn/raster/rt_core/rt_api.c#L3837

CUnit tests here:  https://github.com/bnordgren/postgis/blob/ri-gen2-svn/raster/test/core/cu_geotransform.c

I have the "setGeotransform" function implemented in rt_pg.c. "getGeotransform" needs to return 1 row with four columns, but I'm not sure how to do that.

[bnordgren]

The implementation in my git repo is now exposed to the SQL interface. This includes:

• A new type (geotransform)
• A getter (ST_Getgeotransform)
• Two setter variants (ST_setgeotransform)

The type :

CREATE TYPE geotransform AS (
        imag     double precision,
        jmag     double precision,
        theta_i  double precision,
        theta_ij double precision,
        xoffset  double precision,
        yoffset  double precision) ;

The getter signature:

ST_GetGeotransform(raster) -> geotransform

The setter signatures:

ST_SetGeotransform(rast raster,
                imag double precision, jmag double precision,
                theta_i double precision, theta_ij double precision,
                xoffset double precision, yoffset double precision) -> raster

ST_SetGeotransform(rast raster, gt geotransform) -> raster

These terms are consistent with the wiki pages, but not necessarily with the other accessors.

• theta_i has been referred to as "rotation"
• imag has been referred to as "pixelwidth"
• jmag has been referred to as "pixelheight"

The current terms for the last two are only intuitive if theta_i=0; and they mean exactly the opposite of the intuitive sense if theta_i= +-pi/2 (e.g., width is height and vice versa. However, imag and jmag always represent the pixel length along the i and j axes, respectively.)

xoffset and yoffset specify the location of the origin (i,j)=(0,0). This is identical to the upper left corner if and only if theta_i=0 and theta_ij=-pi/2. Otherwise, it could be any of the other three corners of the raster.

I'm flexible on the terminology, but whatever is chosen should be meticulously docmented, probably in the "raster best practices" section suggested by robe. Reasonable defaults could also be provided.

 https://github.com/bnordgren/postgis/tree/ri-gen2-svn

[robe]

Bryce,

If the ST_GetGeoTransform is the only thing that will be using geotransform, please use OUT parameters instead of a type. Types are a real pain to add in upgrades or change in PostgreSQL. I already yelled at other raster folks for their over love of types.

RETURNS TABLE would be fine too (probably clearer) but all the other code uses OUT parameters so we should probably stick with that convention for now. also the fact I don't have provision in docs for documenting RETURNS TABLE.

[bnordgren]

Oh my. that is nicer. Forgive my inexperience with postgresql implementation concepts.

Done.

ST_GetGeotransform(rast raster,
        OUT imag     double precision,
        OUT jmag     double precision,
        OUT theta_i  double precision,
        OUT theta_ij double precision,
        OUT xoffset  double precision,
        OUT yoffset  double precision)

ST_SetGeotransform(rast raster,
                imag double precision, jmag double precision,
                theta_i double precision, theta_ij double precision,
                xoffset double precision, yoffset double precision) -> raster

[bnordgren]

I think that #1556 duplicates this ticket, but I haven't had a chance to actually look at what's in SVN yet.

[bnordgren]

Yup, looks like this has been covered in #1556