return raster->ipY;
}
+void
+rt_raster_get_phys_params(rt_raster rast,
+ double *i_mag, double *j_mag, double *theta_i, double *theta_ij)
+{
+ double o11, o12, o21, o22 ; /* geotransform coefficients */
+
+ if (rast == NULL) return ;
+ if ( (i_mag==NULL) || (j_mag==NULL) || (theta_i==NULL) || (theta_ij==NULL))
+ return ;
+
+ /* retrieve coefficients from raster */
+ o11 = rt_raster_get_x_scale(rast) ;
+ o12 = rt_raster_get_x_skew(rast) ;
+ o21 = rt_raster_get_y_skew(rast) ;
+ o22 = rt_raster_get_y_scale(rast) ;
+
+ rt_raster_calc_phys_params(o11, o12, o21, o22, i_mag, j_mag, theta_i, theta_ij);
+}
+
+void
+rt_raster_calc_phys_params(double xscale, double xskew, double yskew, double yscale,
+ double *i_mag, double *j_mag, double *theta_i, double *theta_ij)
+
+{
+ double theta_test ;
+
+ if ( (i_mag==NULL) || (j_mag==NULL) || (theta_i==NULL) || (theta_ij==NULL))
+ return ;
+
+ /* pixel size in the i direction */
+ *i_mag = sqrt(xscale*xscale + yskew*yskew) ;
+
+ /* pixel size in the j direction */
+ *j_mag = sqrt(xskew*xskew + yscale*yscale) ;
+
+ /* Rotation
+ * ========
+ * Two steps:
+ * 1] calculate the magnitude of the angle between the x axis and
+ * the i basis vector.
+ * 2] Calculate the sign of theta_i based on the angle between the y axis
+ * and the i basis vector.
+ */
+ *theta_i = acos(xscale/(*i_mag)) ; /* magnitude */
+ theta_test = acos(yskew/(*i_mag)) ; /* sign */
+ if (theta_test < M_PI_2){
+ *theta_i = -(*theta_i) ;
+ }
+
+
+ /* Angular separation of basis vectors
+ * ===================================
+ * Two steps:
+ * 1] calculate the magnitude of the angle between the j basis vector and
+ * the i basis vector.
+ * 2] Calculate the sign of theta_ij based on the angle between the
+ * perpendicular of the i basis vector and the j basis vector.
+ */
+ *theta_ij = acos(((xscale*xskew) + (yskew*yscale))/((*i_mag)*(*j_mag))) ;
+ theta_test = acos( ((-yskew*xskew)+(xscale*yscale)) /
+ ((*i_mag)*(*j_mag)));
+ if (theta_test > M_PI_2) {
+ *theta_ij = -(*theta_ij) ;
+ }
+}
+
+void
+rt_raster_set_phys_params(rt_raster rast,double i_mag, double j_mag, double theta_i, double theta_ij)
+{
+ double o11, o12, o21, o22 ; /* calculated geotransform coefficients */
+ int success ;
+
+ if (rast == NULL) return ;
+
+ success = rt_raster_calc_gt_coeff(i_mag, j_mag, theta_i, theta_ij,
+ &o11, &o12, &o21, &o22) ;
+
+ if (success) {
+ rt_raster_set_scale(rast, o11, o22) ;
+ rt_raster_set_skews(rast, o12, o21) ;
+ }
+}
+
+int
+rt_raster_calc_gt_coeff(double i_mag, double j_mag, double theta_i, double theta_ij,
+ double *xscale, double *xskew, double *yskew, double *yscale)
+{
+ double f ; /* reflection flag 1.0 or -1.0 */
+ double k_i ; /* shearing coefficient */
+ double s_i, s_j ; /* scaling coefficients */
+ double cos_theta_i, sin_theta_i ;
+
+ if ( (xscale==NULL) || (xskew==NULL) || (yskew==NULL) || (yscale==NULL)) {
+ return 0;
+ }
+
+ if ( (theta_ij == 0.0) || (theta_ij == M_PI)) {
+ return 0;
+ }
+
+ /* Reflection across the i axis */
+ f=1.0 ;
+ if (theta_ij < 0) {
+ f = -1.0;
+ }
+
+ /* scaling along i axis */
+ s_i = i_mag ;
+
+ /* shearing parallel to i axis */
+ k_i = tan(f*M_PI_2 - theta_ij) ;
+
+ /* scaling along j axis */
+ s_j = j_mag / (sqrt(k_i*k_i + 1)) ;
+
+ /* putting it altogether */
+ cos_theta_i = cos(theta_i) ;
+ sin_theta_i = sin(theta_i) ;
+ *xscale = s_i * cos_theta_i ;
+ *xskew = k_i * s_j * f * cos_theta_i + s_j * f * sin_theta_i ;
+ *yskew = -s_i * sin_theta_i ;
+ *yscale = -k_i * s_j * f * sin_theta_i + s_j * f * cos_theta_i ;
+ return 1;
+}
+
int32_t
rt_raster_get_srid(rt_raster raster) {
*/
double rt_raster_get_y_skew(rt_raster raster);
+/**
+* Calculates and returns the physically significant descriptors embodied
+* in the geotransform attached to the provided raster.
+*
+* @param rast the raster containing the geotransform of interest
+* @param i_mag size of a pixel along the transformed i axis
+* @param j_mag size of a pixel along the transformed j axis
+* @param theta_i angle by which the raster is rotated (radians positive clockwise)
+* @param theta_ij angle from transformed i axis to transformed j axis
+* (radians positive counterclockwise)
+*
+*/
+void rt_raster_get_phys_params(rt_raster rast,
+ double *i_mag, double *j_mag, double *theta_i, double *theta_ij) ;
+
+/**
+* Calculates the geotransform coefficients and applies them to the
+* supplied raster. The coefficients will not be applied if there was an
+* error during the calculation.
+*
+* This method affects only the scale and skew coefficients. The offset
+* parameters are not changed.
+*
+* @param rast the raster in which the geotransform will be stored.
+* @param i_mag size of a pixel along the transformed i axis
+* @param j_mag size of a pixel along the transformed j axis
+* @param theta_i angle by which the raster is rotated (radians positive clockwise)
+* @param theta_ij angle from transformed i axis to transformed j axis
+* (radians positive counterclockwise)
+*/
+void rt_raster_set_phys_params(rt_raster rast,
+ double i_mag, double j_mag, double theta_i, double theta_ij) ;
+
+
+/**
+* Calculates the physically significant descriptors embodied in an
+* arbitrary geotransform. Always succeeds unless one or more of the
+* output pointers is set to NULL.
+*
+* @param xscale geotransform coefficient o_11
+* @param xskew geotransform coefficient o_12
+* @param yskew geotransform coefficient o_21
+* @param yscale geotransform coefficient o_22
+* @param i_mag size of a pixel along the transformed i axis
+* @param j_mag size of a pixel along the transformed j axis
+* @param theta_i angle by which the raster is rotated (radians positive clockwise)
+* @param theta_ij angle from transformed i axis to transformed j axis
+* (radians positive counterclockwise)
+*/
+void rt_raster_calc_phys_params(double xscale,
+ double xskew, double yskew, double yscale,
+ double *i_mag, double *j_mag, double *theta_i, double *theta_ij) ;
+
+
+/**
+* Calculates the coefficients of a geotransform given the physically
+* significant parameters describing the transform. Will fail if any of the
+* result pointers is NULL, or if theta_ij has an illegal value (0 or PI).
+*
+* @param i_mag size of a pixel along the transformed i axis
+* @param j_mag size of a pixel along the transformed j axis
+* @param theta_i angle by which the raster is rotated (radians positive clockwise)
+* @param theta_ij angle from transformed i axis to transformed j axis
+* (radians positive counterclockwise)
+* @param xscale geotransform coefficient o_11
+* @param xskew geotransform coefficient o_12
+* @param yskew geotransform coefficient o_21
+* @param yscale geotransform coefficient o_22
+* @return 1 if the calculation succeeded, 0 if error.
+*/
+int rt_raster_calc_gt_coeff(double i_mag,
+ double j_mag, double theta_i, double theta_ij,
+ double *xscale, double *xskew, double *yskew, double *yscale) ;
+
/**
* Set raster's SRID
*
Datum RASTER_getYUpperLeft(PG_FUNCTION_ARGS);
Datum RASTER_getPixelWidth(PG_FUNCTION_ARGS);
Datum RASTER_getPixelHeight(PG_FUNCTION_ARGS);
-Datum RASTER_getRotation(PG_FUNCTION_ARGS);
+Datum RASTER_getGeotransform(PG_FUNCTION_ARGS);
/* Set all the properties of a raster */
Datum RASTER_setSRID(PG_FUNCTION_ARGS);
Datum RASTER_setSkewXY(PG_FUNCTION_ARGS);
Datum RASTER_setUpperLeftXY(PG_FUNCTION_ARGS);
Datum RASTER_setRotation(PG_FUNCTION_ARGS);
+Datum RASTER_setGeotransform(PG_FUNCTION_ARGS);
/* Get all the properties of a raster band */
Datum RASTER_getBandPixelType(PG_FUNCTION_ARGS);
POSTGIS_RT_DEBUGF(3, "call number %d", call_cntr);
nulls = palloc(sizeof(bool) * values_length);
- memset(nulls, FALSE, values_length);
+ memset(nulls, FALSE, values_length);
values[0] = CStringGetTextDatum(geomval2[call_cntr].geom);
values[1] = Float8GetDatum(geomval2[call_cntr].val);
}
/**
- * Return the raster rotation. The raster rotation is calculated from
- * the scale and skew values stored in the georeference. If the scale
- * and skew values indicate that the raster is not uniformly rotated
- * (the pixels are diamond-shaped), this function will return NaN.
+ * Set the geotransform of the supplied raster. Returns the raster.
*/
-PG_FUNCTION_INFO_V1(RASTER_getRotation);
-Datum RASTER_getRotation(PG_FUNCTION_ARGS)
+PG_FUNCTION_INFO_V1(RASTER_setGeotransform);
+Datum RASTER_setGeotransform(PG_FUNCTION_ARGS)
{
- rt_pgraster *pgraster;
- rt_raster raster;
- double xscale, xskew, yscale, yskew, xrot, yrot;
+ rt_pgraster *pgraster ;
+ rt_raster raster ;
+ float8 imag, jmag, theta_i, theta_ij, xoffset, yoffset ;
- if (PG_ARGISNULL(0)) PG_RETURN_NULL();
- pgraster = (rt_pgraster *)PG_DETOAST_DATUM_SLICE(PG_GETARG_DATUM(0), 0, sizeof(struct rt_raster_serialized_t));
+ if (PG_ARGISNULL(0) || PG_ARGISNULL(1) || PG_ARGISNULL(2) ||
+ PG_ARGISNULL(3) || PG_ARGISNULL(4) ||
+ PG_ARGISNULL(5) || PG_ARGISNULL(6))
+ PG_RETURN_NULL();
+
+ /* get the inputs */
+ pgraster = (rt_pgraster *)PG_DETOAST_DATUM_SLICE(PG_GETARG_DATUM(0),
+ 0, sizeof(struct rt_raster_serialized_t));
+ imag = PG_GETARG_FLOAT8(1) ;
+ jmag = PG_GETARG_FLOAT8(2) ;
+ theta_i = PG_GETARG_FLOAT8(3);
+ theta_ij = PG_GETARG_FLOAT8(4);
+ xoffset = PG_GETARG_FLOAT8(5);
+ yoffset = PG_GETARG_FLOAT8(6);
raster = rt_raster_deserialize(pgraster, TRUE);
if (!raster) {
- elog(ERROR, "RASTER_getRotation: Could not deserialize raster");
+ elog(ERROR, "RASTER_setGeotransform: Could not deserialize raster");
PG_RETURN_NULL();
}
- xscale = rt_raster_get_x_scale(raster);
- yscale = rt_raster_get_y_scale(raster);
+ /* store the new geotransform */
+ rt_raster_set_phys_params(raster, imag,jmag,theta_i,theta_ij) ;
+ rt_raster_set_offsets(raster, xoffset, yoffset) ;
- if (xscale == 0 || yscale == 0) {
- rt_raster_destroy(raster);
+ /* prep the return value */
+ pgraster = rt_raster_serialize(raster);
+ if ( ! pgraster ) PG_RETURN_NULL();
- /* cannot compute scale with a zero denominator */
- elog(NOTICE, "RASTER_getRotation: Could not divide by zero scale; cannot determine raster rotation.");
- PG_RETURN_FLOAT8(NAN);
- }
+ SET_VARSIZE(pgraster, pgraster->size);
- xskew = rt_raster_get_x_skew(raster);
- yskew = rt_raster_get_y_skew(raster);
+ rt_raster_destroy(raster);
+
+ PG_RETURN_POINTER(pgraster);
+}
- xrot = atan(yskew/xscale);
- yrot = atan(xskew/yscale);
- rt_raster_destroy(raster);
+/**
+ * Calculates the physically relevant parameters of the supplied raster's
+ * geotransform. Returns them as a set.
+ */
+PG_FUNCTION_INFO_V1(RASTER_getGeotransform);
+Datum RASTER_getGeotransform(PG_FUNCTION_ARGS)
+{
+ rt_pgraster *rast ;
+ rt_raster raster ;
+ float8 imag, jmag, theta_i, theta_ij, xoffset, yoffset ;
+ TupleDesc result_tuple ; /* for returning a composite */
+ HeapTuple heap_tuple ; /* instance of the tuple to return */
+ Oid result_oid ; /* internal code for the specific return type */
+ TypeFuncClass return_type ; /* is the return type a composite? */
+ Datum return_values[6] ;
+ bool nulls[6] ;
+
+ /* setup the return value infrastructure */
+ return_type = get_call_result_type(fcinfo, &result_oid, &result_tuple) ;
+ if (return_type != TYPEFUNC_COMPOSITE) {
+ rterror("RASTER_getGeotransform(): function returning record called in context that cannot accept type record") ;
+ PG_RETURN_NULL() ;
+ }
+ result_tuple = BlessTupleDesc(result_tuple) ;
+
+ /* get argument */
+ if (PG_ARGISNULL(0)) PG_RETURN_NULL();
+ rast = (rt_pgraster *)PG_DETOAST_DATUM_SLICE(PG_GETARG_DATUM(0),
+ 0, sizeof(struct rt_raster_serialized_t));
- if (FLT_EQ(xrot, yrot)) {
- PG_RETURN_FLOAT8(xrot);
+ raster = rt_raster_deserialize(rast, TRUE);
+ if (!raster) {
+ elog(ERROR, "RASTER_detGeotransform: Could not deserialize raster");
+ PG_RETURN_NULL();
}
- PG_RETURN_FLOAT8(NAN);
+ /* do the calculation */
+ rt_raster_calc_phys_params(
+ rt_raster_get_x_scale(raster),
+ rt_raster_get_x_skew(raster),
+ rt_raster_get_y_skew(raster),
+ rt_raster_get_y_scale(raster),
+ &imag, &jmag, &theta_i, &theta_ij) ;
+
+ rt_raster_destroy(raster);
+
+ /* prep the composite return value */
+ /* construct datum array */
+ return_values[0] = Float8GetDatum(imag) ;
+ return_values[1] = Float8GetDatum(jmag) ;
+ return_values[2] = Float8GetDatum(theta_i) ;
+ return_values[3] = Float8GetDatum(theta_ij) ;
+ return_values[4] = Float8GetDatum(rt_raster_get_x_offset(raster)) ;
+ return_values[5] = Float8GetDatum(rt_raster_get_y_offset(raster)) ;
+ memset(nulls, FALSE, 6);
+
+ /* stick em on the heap */
+ heap_tuple = heap_form_tuple(result_tuple, return_values, nulls) ;
+
+ /* return */
+ PG_RETURN_DATUM(HeapTupleGetDatum(heap_tuple)) ;
}
+
/**
* Set the rotation of the raster. This method will change the X Scale,
* Y Scale, X Skew and Y Skew properties all at once to keep the rotations
rt_pgraster *pgraster = NULL;
rt_raster raster;
double rotation = PG_GETARG_FLOAT8(1);
- double xscale, yscale, xskew, yskew, psize;
-
- if (PG_ARGISNULL(0)) PG_RETURN_NULL();
- pgraster = (rt_pgraster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
+ double imag, jmag, theta_i, theta_ij;
- /* no matter what, we don't rotate more than once around */
- if (rotation < 0) {
- rotation = (-2*M_PI) + fmod(rotation, (2*M_PI));
- }
- else {
- rotation = fmod(rotation, (2 * M_PI));
- }
+ if (PG_ARGISNULL(0)) PG_RETURN_NULL();
+ pgraster = (rt_pgraster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
raster = rt_raster_deserialize(pgraster, FALSE);
if (! raster ) {
PG_RETURN_NULL();
}
- xscale = rt_raster_get_x_scale(raster);
- yskew = rt_raster_get_y_skew(raster);
- psize = sqrt(xscale*xscale + yskew*yskew);
- xscale = psize * cos(rotation);
- yskew = psize * sin(rotation);
-
- yscale = rt_raster_get_y_scale(raster);
- xskew = rt_raster_get_x_skew(raster);
- psize = sqrt(yscale*yscale + xskew*xskew);
- yscale = psize * cos(rotation);
- xskew = psize * sin(rotation);
-
- rt_raster_set_scale(raster, xscale, yscale);
- rt_raster_set_skews(raster, xskew, yskew);
+ /* preserve all defining characteristics of the grid except for rotation */
+ rt_raster_get_phys_params(raster, &imag, &jmag, &theta_i, &theta_ij);
+ rt_raster_set_phys_params(raster, imag, jmag, rotation, theta_ij);
pgraster = rt_raster_serialize(raster);
rt_raster_destroy(raster);
AS 'MODULE_PATHNAME', 'RASTER_getPixelHeight'
LANGUAGE 'C' IMMUTABLE STRICT;
+CREATE TYPE geotransform AS (
+ imag double precision,
+ jmag double precision,
+ theta_i double precision,
+ theta_ij double precision,
+ xoffset double precision,
+ yoffset double precision);
+
+CREATE OR REPLACE FUNCTION st_geotransform(raster)
+ RETURNS geotransform
+ AS 'MODULE_PATHNAME', 'RASTER_getGeotransform'
+ LANGUAGE 'C' IMMUTABLE;
+
CREATE OR REPLACE FUNCTION st_rotation(raster)
RETURNS float8
- AS 'MODULE_PATHNAME','RASTER_getRotation'
- LANGUAGE 'C' IMMUTABLE STRICT;
+ AS $$ SELECT (ST_Geotransform($1)).theta_i $$
+ LANGUAGE 'sql' VOLATILE;
CREATE OR REPLACE FUNCTION st_metadata(
rast raster,
AS 'MODULE_PATHNAME','RASTER_setRotation'
LANGUAGE 'C' IMMUTABLE STRICT;
+CREATE OR REPLACE FUNCTION st_setgeotransform(rast raster,
+ imag double precision,
+ jmag double precision,
+ theta_i double precision,
+ theta_ij double precision,
+ xoffset double precision,
+ yoffset double precision)
+ RETURNS raster
+ AS 'MODULE_PATHNAME','RASTER_setGeotransform'
+ LANGUAGE 'C' IMMUTABLE;
+
+CREATE OR REPLACE FUNCTION st_setgeotransform(rast raster, gt geotransform)
+ RETURNS raster
+ AS $$ SELECT st_setgeotransform($1, ($2).imag, ($2).jmag,
+ ($2).theta_i, ($2).theta_ij,
+ ($2).xoffset, ($2).yoffset) $$
+ LANGUAGE 'sql' VOLATILE;
+
-----------------------------------------------------------------------
-- Raster Editors ST_SetGeoreference()
-----------------------------------------------------------------------
) = 0
FROM ST_TestRasterNgb(3, 3, 10) AS rast;
--- test st_stddev4ma, all different values
+-- test st_stddev4ma, values 1-9
SELECT
ST_Value(rast, 2, 2) = 5,
round(ST_Value(
-4|1x1, ip:7.5,2.5 scale:5,5 skew:1,1, srid:-1, width:1, height:1|0.197395559849881
-5|1x1, ip:7.5,2.5 scale:5,5 skew:3,7, srid:-1, width:1, height:1|NaN
-104|5.09901951359278|5.09901951359278|0|0|0
-105|8.60232526704263|5.8309518948453|0|0|0
-104|3.60555127546399|3.60555127546399|3.60555127546399|3.60555127546399|0.785398163397448
-105|6.08276253029822|4.12310562561766|4.12310562561766|6.08276253029822|0.785398163397448
-104|4.41588043316392|4.41588043316392|2.54950975679639|2.54950975679639|0.523598775598299
-105|7.44983221287567|5.04975246918104|2.91547594742265|4.30116263352131|0.523598775598299
+4|1x1, ip:7.5,2.5 scale:5,5 skew:1,1, srid:-1, width:1, height:1|-0.19739555984988
+5|1x1, ip:7.5,2.5 scale:5,5 skew:3,7, srid:-1, width:1, height:1|-0.950546840812075
+104|5.09901951359278|4.70678724331642|1.96116135138184|-0|0
+105|8.60232526704263|0.464990554975281|5.81238193719096|-0|0
+104|3.60555127546399|1.9414506867883|4.71495166791445|-3.60555127546399|0.785398163397448
+105|6.08276253029822|-3.78117670802321|4.43877265724465|-6.08276253029822|0.785398163397448
+104|4.41588043316392|3.09561664722963|4.0518091728751|-2.54950975679639|0.523598775598298
+105|7.44983221287567|-2.50349733546706|5.26616569159282|-4.30116263352131|0.523598775598299
projects / postgis / commitdiff