[postgis-users] Function - Standard deviation ellipses

Francis Dupont Francis.X.Dupont at USherbrooke.ca
Wed Jun 25 10:50:35 PDT 2008


I created a function that generate a standard deviation ellipse based on
crimestat function in chapter 4.  I thought the community might be interested
so this is why I post it!!  Maybe it's not the right place to post the
function, if so let me know where it should go.

Before executing the function, you need to insert your point data in the temp
table.  I only tried the function with Quebec Lambert projection.

Feel free (this is essential!) to contribute and optimized the function.


Francis Dupont
francis.x.dupont at usherbrooke.ca


You will need that table to use the function generate_ellipse

CREATE TABLE tmp_data_ellipse
  x numeric NOT NULL,
  y numeric NOT NULL,
  id serial NOT NULL,
  CONSTRAINT pkid_tmp_data_ellipse PRIMARY KEY (id)

Function: std_dist_ellipse

    Create a standard distance ellipse, with major and minor axes reflecting the
directional variation of the point pattern store in tmp_data_ellipse table.


    num_of_std :    Number of standard deviation to generate for each group of
                    e.g.:  A value of 2 mean that an ellipse of distance of 1
                    deviation and 2 standard deviation will be created for each

    fix_point :     0 egual no action otherwise 1 mean that you want to fix the
                    of point at the same coordinate.  If, for some reason, all
the points
                    are at the same coordinate the differences between the value
and the mean
                    will be 0.  This will give a division by zero error which
will be catch
                    by the exception handler.  To avoid this problem, the fix
point option
                    let you move the point that have the same coordinates.  Each
of those points
                    are move randomly on a circle of 1 meter radius.


    Recordset of the ellipses created with the corresponding standard deviation.


CREATE OR REPLACE FUNCTION std_dist_ellipse(num_of_std smallint, fix_point
    stats record;
    data record;
    sum_square_dif_avg_x numeric;   --SUM [(Xi - avg(X))^2]
    sum_square_dif_avg_y numeric;   --SUM [(Yi - avg(Y))^2]
    sum_dif_avg_x_y numeric;        --SUM [(Xi - avg(X) * (Yi - avg(Y)]
    sum_square_dif_avg_x_y numeric; --SUM [((Xi - avg(X) * (Yi - avg(Y))^2]
    C numeric;                      --Constant used (see standard ellipse
distance reference for more information)
    theta numeric;                  --Angle calculated in rad
    SDx numeric;                    --Standard deviation in X
    SDy numeric;                    --Standard deviation in Y
    SDx_sum_x_y_cos_sin_theta numeric;
    SDy_sum_x_y_sin_cos_theta numeric;
    stdev_rank smallint;            --Rank of the standard deviation being
    LastX numeric;                  --Value of the last X being verified (use in
fix_point option)
    LastY numeric;                  --Value of the last Y being verified (use in
fix_point option)
    r numeric;                      --radius of the circule (use in fix_point

        For more information on the formula used in this function, take a look
at the following websites:

        Spatial Analysis Online:   
        CrimeStat (chapter 4):      http://www.icpsr.umich.edu/CRIMESTAT/

        A somewhat different form of directional analysis is used within many
packages to summarise the distribution
        of point sets around a mean centre. The average variation in the
distance of points around the mean can be viewed as
        a circle or set of circles at a set of standard distances (like standard
deviations in univariate statistics).
        In addition, separate variations in x- and y-coordinate values may be
used to generate a standard distance
        ellipse, with major and minor axes reflecting the directional variation
of the point pattern.

    --Check if fix point option is activated
    IF (fix_point = 1) THEN

        --Reset value so there's no changed for the first point
        LastX := NULL;
        LastY := NULL;
        r     := 1;     --set a 1 meter radius circle

        FOR data IN SELECT id, x, y FROM tmp_data_ellipse ORDER BY x, y LOOP

            --Check if the last point have the same coordinate as the last one
            --Check null value to avoid calculation for the first point.
            IF ((LastX = data.x) AND (LastY = data.y)) AND ((LastX IS NOT NULL)

                theta := random() * 2 * pi();   --angle value in radian between
0 and 2pi

                --Move the point on a circle of r radius and update the point in
the temporary tmp_data_ellipse
                UPDATE tmp_data_ellipse
                SET x = data.x + ((cos(theta))*r),
                    y = data.y + ((sin(theta))*r)
                WHERE id = data.id;

                --Save the last point being processed
                LastX := data.x;
                LastY := data.y;
            END IF;

        END LOOP;

    END IF;

    --Calculate the average of the absolute values of X and Y, the X and Y value
of the centroid in the projection used and the number of points
    SELECT avg(abs(x)) as avg_x, avg(abs(y)) as avg_y, avg(x) as center_x,
avg(y) as center_y, COUNT(*) as n FROM tmp_data_ellipse INTO stats;

    --Reset all var values before calculation
    sum_square_dif_avg_x    := 0;
    sum_square_dif_avg_y    := 0;
    sum_dif_avg_x_y         := 0;
    sum_square_dif_avg_x_y  := 0;

    --Pre-calculate the values used in formula.  It is necessary to used the
absolute values in calculation to be in
    --a trigonometric circle otherwise you'll get the wrong angle when the angle
value calculated is beyond pi.
    --The negative values in X or Y depends on the cartographic projection you
are using.
    FOR data IN SELECT x, y FROM tmp_data_ellipse LOOP

        sum_square_dif_avg_x := sum_square_dif_avg_x + ((abs(data.x) -
        sum_square_dif_avg_y := sum_square_dif_avg_y + ((abs(data.y) -
        sum_dif_avg_x_y := sum_dif_avg_x_y + ((abs(data.x) - stats.avg_x) *
(abs(data.y) - stats.avg_y));
        sum_square_dif_avg_x_y := sum_square_dif_avg_x_y + (((abs(data.x) -
stats.avg_x) * (abs(data.y) - stats.avg_y))^2);


    --Calculate the square root of the C value (as in the formula)
    C := sqrt(((sum_square_dif_avg_x - sum_square_dif_avg_y)^2) + (4 *

    --Calulate the angle of the ellipse
    theta := atan((sum_square_dif_avg_x - sum_square_dif_avg_y + C) / (2 *

    --Reset values and calculate the standard deviation in X and Y for the
    SDx_sum_x_y_cos_sin_theta := 0;
    SDy_sum_x_y_sin_cos_theta := 0;

    --Use the var center_x and center_y because they are in projection
coordinate values.
    FOR data IN SELECT x, y FROM tmp_data_ellipse LOOP

        SDx_sum_x_y_cos_sin_theta := SDx_sum_x_y_cos_sin_theta + ((((data.x -
stats.center_x) * cos(theta))- ((data.y - stats.center_y) * sin(theta)))^2);
        SDy_sum_x_y_sin_cos_theta := SDy_sum_x_y_sin_cos_theta + ((((data.x -
stats.center_x) * sin(theta))- ((data.y - stats.center_y) * cos(theta)))^2);


    --Note, again, that 2 is subtracted from the
    --number of points in both denominators to produce an unbiased estimate of
    --the standard deviational ellipse since there are two constants from which
    --distance along each axis is measured (mean of X, mean of Y) (Crimestat,
    SDx := sqrt((2 * SDx_sum_x_y_cos_sin_theta) / (stats.n - 2));
    SDy := sqrt((2 * SDy_sum_x_y_sin_cos_theta) / (stats.n - 2));

    --Temporary insert for validating the result
    INSERT INTO tmp_data_ellipse(x, y)
    VALUES(sum_square_dif_avg_x, sum_square_dif_avg_y);
    INSERT INTO tmp_data_ellipse(x, y)
    VALUES(sum_dif_avg_x_y, sum_square_dif_avg_x_y);
    INSERT INTO tmp_data_ellipse(x, y)
    VALUES(C, C);
    INSERT INTO tmp_data_ellipse(x, y)
    VALUES(theta * 180 / pi(), theta);
    INSERT INTO tmp_data_ellipse(x, y)
    VALUES(SDx, SDy);

    --Produce an ellipse for each number of standard deviation.  Start at
num_of_std (usually at 2 standard deviation) downto 1
    stdev_rank := num_of_std;
    WHILE stdev_rank >= 1 LOOP
        SELECT ellipse(stats.center_x, stats.center_y, (stdev_rank * SDx),
(stdev_rank * SDy), theta, 40) as the_geom, stdev_rank INTO data;
        RETURN NEXT data;
        stdev_rank := stdev_rank - 1;   --go to next rank


    --Stop the process if a division by zero occured
    WHEN DIVISION_BY_ZERO THEN RAISE EXCEPTION 'error: Division by zero, check
the coordinates of your point dataset.';

More information about the postgis-users mailing list