Computing radius of gyration in PostGIS

Question

The radius of gyration is a measure of the standard deviation of distances between points (and their center of mass), i.e., how frequently and how far points are moving. A formal definition is:

I have a large table of localization data (x-y coordinates, not a geometry field - might be transformed if faster?). What is the best way to compute the radius of gyration with PostGreSQL / PostGIS?

For the center of mass, ST_Centroid(ST_Union(ST_MakePoint(x,y)))? Then ST_Distance? How to aggregate all this? ST_Distance(MakePoint(x,y),ST_Centroid(ST_Union(ST_MakePoint(x,y)))) does not work (ERROR: column "x" must appear in the GROUP BY clause or be used in an aggregate function).

Sorry, maybe a stupid question, I'm not so familiar with PostGreSQL.

Answer 1

I think it is better to store your points as geometries because you can use a spatial index for speeding up queries.

The query:

WITH points AS
    (SELECT
        ST_SetSRID(ST_MakePoint(x,y), 4326) AS geom
    FROM your_xy_table
    ),
centroid AS
    (SELECT
        ST_Centroid(ST_Union(points.geom)) AS geom
    FROM points),
multiobject AS
    (SELECT
        geom,
        generate_series(1,600) AS n 
    FROM points),
objects AS
    (SELECT
        n, 
        ST_GeometryN(geom, n) AS geometries 
    FROM multiobject 
        WHERE n <= ST_NumGeometries(geom)),
distance AS
    (SELECT DISTINCT ON
        (ST_Distance(geometries, centroid.geom))
        ST_Distance(geometries, centroid.geom) as distance
    FROM objects, centroid)

    SELECT
        stddev(distance) AS gyration
    FROM distance

First subquery: creates the points

Second subquery:calculates the centroid of your points

Third subquery: generates a series for your point objects, you have to change the value, if your data is bigger than 1000 objects

Fourth subquery: selects the single geometries for the distance calculation

Fifth subquery: calculates the distances

The last step is the calculation of the standard deviation of the distance values.

UPDATE

A second approach is this shorter query. The first subquery creates the points with the xy data in your table (but it is better to store the points as geometries, like mentioned above). You have to put in the EPSG number of your projection (e.g. 4326) in the first query. Change it when you work with a different projection.

WITH points AS
    (SELECT
        ST_SetSRID(ST_MakePoint(x,y), 4326) AS geom
    FROM your_xy_table),
centroid AS
    (SELECT
        ST_Centroid(ST_Union(points.geom)) AS geom
    FROM points),
objects AS
    (SELECT
        (ST_Dump(points.geom)).geom AS geometries
    FROM points),
distance AS
    (SELECT DISTINCT ON
        (ST_Distance(geometries, centroid.geom))
        ST_Distance(geometries, centroid.geom) AS distance
    FROM objects, centroid)

SELECT stddev(distance) FROM distance

With the first query maybe you can make other detailed calculations because you handle with the single objects with their numbers (n).

Answer 2

The formula appears correct, but I don't think the implementation is correct. Taking the standard deviation of the distance between each point and the centroid using the stddev function in PostgreSQL or other statistical software calculates the difference between each distance measure and the mean of the distance measures. In other words, the (r_i - r_{cm}) term needs to represent the distance between each point and the centroid, rather than (dist - avg(dist)). See this reproducible example. All steps are broken down in common table expressions for the sake of transparency.

--how many random points?
with n as (select 50 as n)
--generate a set of points using random XY
, xy as (select generate_series(1, n) as gid, st_makepoint(random()*1000, random()*5000) as geom_point from n)
--generate the centroid
, c as (select st_centroid(st_collect(geom_point)) as geom_centroid from xy)
--distance between each point and centroid
, d as (select gid, st_distance(geom_point, geom_centroid) as dist from xy, c)
--standard deviation of distance
, d_sd as (select stddev(dist) from d)
--square the distance
, d2 as (select gid, dist, dist^2 as dist_squared from d)
--sum the squares and divide by n 
, d2sum as (select sum(dist_squared) as dist_squared_sum from d2)
--divide by n
, d2sum_normalize as (select dist_squared_sum / n as dist_norm from d2sum, n)
--square root
, rog as (select sqrt(dist_norm) as radius_of_gyration from d2sum_normalize)
--report radius of gyration and stddev(distance)
select * from rog, d_sd