ST_Project does not project perfectly east-west

[petermj]

I am using ST_Project to project a point east-west and/or north-south.

If I project the point north/south the latitude remains the same, if I project east/west the longitude also changes.

Simple example :

SELECT 
    -- ST_Project(pt, distance, 0) move pt by distance north
    -- ST_Project(pt, distance, radians(90)) move pt by distance east
    -- ST_Project(pt, distance, radians(180) move pt by distance south
    -- ST_Project(pt, distance, radians(270)) move pt by distance west
count,
ST_AsEWKT( (( ST_Project( st_geogfromtext('SRID=4326; POINT(-4.5 63)'), count*1000.0, radians(0))))::geometry),
ST_AsEWKT( (( ST_Project( st_geogfromtext('SRID=4326; POINT(-4.5 63)'), count*1000.0, radians(180))))::geometry),
ST_AsEWKT( (( ST_Project( st_geogfromtext('SRID=4326; POINT(-4.5 63)'), count*1000.0, pi())))::geometry),
ST_AsEWKT( (( ST_Project( st_geogfromtext('SRID=4326; POINT(-4.5 63)'), count*1000.0, pi()/2)))::geometry),
ST_AsEWKT( (( ST_Project( st_geogfromtext('SRID=4326; POINT(-4.5 63)'), count*1000.0, 3*pi()/2)))::geometry),
ST_AsEWKT( (( ST_Project( st_geogfromtext('SRID=4326; POINT(-4.5 63)'), count*1000.0, radians(90))))::geometry),
ST_AsEWKT( (( ST_Project( st_geogfromtext('SRID=4326; POINT(-4.5 63)'), count*1000.0, radians(270))))::geometry)
    FROM generate_series(-5, 5) as count
    ORDER BY count
;

Produces:

-5,'SRID=4326;POINT(-4.5 62.9551413817306)','SRID=4326;POINT(-4.5 63.0448583314294)','SRID=4326;POINT(-4.5 63.0448583314294)','SRID=4326;POINT(-4.59867214178292 62.9999655834464)','SRID=4326;POINT(-4.40132785821708 62.9999655834464)','SRID=4326;POINT(-4.59867214178292 62.9999655834464)','SRID=4326;POINT(-4.40132785821708 62.9999655834464)'
-4,'SRID=4326;POINT(-4.5 62.9641131283474)','SRID=4326;POINT(-4.5 63.0358866880749)','SRID=4326;POINT(-4.5 63.0358866880749)','SRID=4326;POINT(-4.57893773572975 62.9999779734007)','SRID=4326;POINT(-4.42106226427025 62.9999779734007)','SRID=4326;POINT(-4.57893773572975 62.9999779734007)','SRID=4326;POINT(-4.42106226427025 62.9999779734007)'
-3,'SRID=4326;POINT(-4.5 62.9730848634802)','SRID=4326;POINT(-4.5 63.0269150332574)','SRID=4326;POINT(-4.5 63.0269150332574)','SRID=4326;POINT(-4.55920331480765 62.9999876100357)','SRID=4326;POINT(-4.44079668519235 62.9999876100357)','SRID=4326;POINT(-4.55920331480765 62.9999876100357)','SRID=4326;POINT(-4.44079668519235 62.9999876100357)'
-2,'SRID=4326;POINT(-4.5 62.9820565871314)','SRID=4326;POINT(-4.5 63.0179433669741)','SRID=4326;POINT(-4.5 63.0179433669741)','SRID=4326;POINT(-4.53946888273384 62.9999944933485)','SRID=4326;POINT(-4.46053111726616 62.9999944933485)','SRID=4326;POINT(-4.53946888273384 62.9999944933485)','SRID=4326;POINT(-4.46053111726616 62.9999944933485)'
-1,'SRID=4326;POINT(-4.5 62.9910282993038)','SRID=4326;POINT(-4.5 63.0089716892226)','SRID=4326;POINT(-4.5 63.0089716892226)','SRID=4326;POINT(-4.51973444322554 62.999998623337)','SRID=4326;POINT(-4.48026555677446 62.999998623337)','SRID=4326;POINT(-4.51973444322554 62.999998623337)','SRID=4326;POINT(-4.48026555677446 62.999998623337)'
0,'SRID=4326;POINT(-4.5 63)','SRID=4326;POINT(-4.5 63)','SRID=4326;POINT(-4.5 63)','SRID=4326;POINT(-4.5 63)','SRID=4326;POINT(-4.5 63)','SRID=4326;POINT(-4.5 63)','SRID=4326;POINT(-4.5 63)'
1,'SRID=4326;POINT(-4.5 63.0089716892226)','SRID=4326;POINT(-4.5 62.9910282993038)','SRID=4326;POINT(-4.5 62.9910282993038)','SRID=4326;POINT(-4.48026555677446 62.999998623337)','SRID=4326;POINT(-4.51973444322554 62.999998623337)','SRID=4326;POINT(-4.48026555677446 62.999998623337)','SRID=4326;POINT(-4.51973444322554 62.999998623337)'
2,'SRID=4326;POINT(-4.5 63.0179433669741)','SRID=4326;POINT(-4.5 62.9820565871314)','SRID=4326;POINT(-4.5 62.9820565871314)','SRID=4326;POINT(-4.46053111726616 62.9999944933485)','SRID=4326;POINT(-4.53946888273384 62.9999944933485)','SRID=4326;POINT(-4.46053111726616 62.9999944933485)','SRID=4326;POINT(-4.53946888273384 62.9999944933485)'
3,'SRID=4326;POINT(-4.5 63.0269150332574)','SRID=4326;POINT(-4.5 62.9730848634802)','SRID=4326;POINT(-4.5 62.9730848634802)','SRID=4326;POINT(-4.44079668519235 62.9999876100357)','SRID=4326;POINT(-4.55920331480765 62.9999876100357)','SRID=4326;POINT(-4.44079668519235 62.9999876100357)','SRID=4326;POINT(-4.55920331480765 62.9999876100357)'
4,'SRID=4326;POINT(-4.5 63.0358866880749)','SRID=4326;POINT(-4.5 62.9641131283474)','SRID=4326;POINT(-4.5 62.9641131283474)','SRID=4326;POINT(-4.42106226427025 62.9999779734007)','SRID=4326;POINT(-4.57893773572975 62.9999779734007)','SRID=4326;POINT(-4.42106226427025 62.9999779734007)','SRID=4326;POINT(-4.57893773572975 62.9999779734007)'
5,'SRID=4326;POINT(-4.5 63.0448583314294)','SRID=4326;POINT(-4.5 62.9551413817306)','SRID=4326;POINT(-4.5 62.9551413817306)','SRID=4326;POINT(-4.40132785821708 62.9999655834464)','SRID=4326;POINT(-4.59867214178292 62.9999655834464)','SRID=4326;POINT(-4.40132785821708 62.9999655834464)','SRID=4326;POINT(-4.59867214178292 62.9999655834464)'

This means that the order of projection then becomes an issue, because if I have a point and I project it east 10km, then north 10km, this will be different to if I project it north 10km, then east 10km.

Example:

SELECT 
    -- ST_Project(pt, distance, 0) move pt by distance north
    -- ST_Project(pt, distance, radians(90)) move pt by distance east
    -- ST_Project(pt, distance, radians(180) move pt by distance south
    -- ST_Project(pt, distance, radians(270)) move pt by distance west
count,
ST_AsEWKT( ST_Project( ST_Project( st_geogfromtext('SRID=4326; POINT(-4.5 63)'), count*10000.0, radians(0)), count*10000.0,  radians(90))),
ST_AsEWKT( ST_Project( ST_Project( st_geogfromtext('SRID=4326; POINT(-4.5 63)'), count*10000.0, radians(90)), count*10000.0, radians(0)))
FROM generate_series(-5, 5) as count
    ORDER BY count
;

Produces:

-5,'SRID=4326;POINT(-5.47176629509741 62.5480247248945)','SRID=4326;POINT(-5.48664476153847 62.5479591959442)'
-4,'SRID=4326;POINT(-5.27978234144984 62.6389539645604)','SRID=4326;POINT(-5.28933810708621 62.6389203106193)'
-3,'SRID=4326;POINT(-5.08662259121748 62.7296192016654)','SRID=4326;POINT(-5.09201658867479 62.7296049603157)'
-2,'SRID=4326;POINT(-4.89227819418138 62.8200173614923)','SRID=4326;POINT(-4.89468392069683 62.8200131288902)'
-1,'SRID=4326;POINT(-4.69674028795425 62.9101453392753)','SRID=4326;POINT(-4.6973438189138 62.9101448085789)'
0,'SRID=4326;POINT(-4.5 63)','SRID=4326;POINT(-4.5 63)'
1,'SRID=4326;POINT(-4.30204844972743 63.089578178207)','SRID=4326;POINT(-4.3026561810862 63.0895787121552)'
2,'SRID=4326;POINT(-4.10287675065965 63.1788766778)','SRID=4326;POINT(-4.10531607930317 63.1788809624297)'
3,'SRID=4326;POINT(-3.90247601268086 63.2678922718583)','SRID=4326;POINT(-3.90798341132521 63.2679067765911)'
4,'SRID=4326;POINT(-3.7008373443622 63.3566217024552)','SRID=4326;POINT(-3.71066189291379 63.3566561887871)'
5,'SRID=4326;POINT(-3.49795185536835 63.4450616804805)','SRID=4326;POINT(-3.51335523846153 63.4451292415414)'

Info string:

'POSTGIS="2.4.1 r16012" PGSQL="100" GEOS="3.6.2-CAPI-1.10.2 4d2925d" PROJ="Rel. 4.9.3, 15 August 2016" GDAL="GDAL 2.2.2, released 2017/09/15" LIBXML="2.7.8" LIBJSON="0.12" LIBPROTOBUF="1.2.1" RASTER'

'PostgreSQL 10.0, compiled by Visual C++ build 1800, 64-bit'

Windows 10, tested inside pgAdmin 4 v2.0

[pramsey]

If I am on a line of latitude in the northern hemisphere and I want to move to a point 10km east and on the same line of latitude, do I set my heading directly east? No, I do not. I set it slightly north, so when I start walking, without deviation, the great circle route lands me at the point I desired.

Similarly, if I am in the northern hemisphere and set my direction due east and start walking, without deviation, I will follow a great circle from that point, which will inevitably lead me southwards.

The only place where I can set my heading due east and start walking and stay at the same latitude is the equator.

[petermj]

I agree with you if this were a geometry calculation, but this is a geography calculation, where north/south should follow lines of latitude, and east/west should follow lines of longitude.

If I move X degrees east, I would not expect my latitude to vary, the ST_Project function projects my movement in metres, so my latitude should not vary.

—-

Alternatively, you may be saying that this should use a geometry calculation, and projecting east, aligned with the cartesian grid, should line up with the longitude axis, and only the latitude should change. But this also does not work:

PostGIS does not seem to calculate this differently for the geometry and geography types.

SELECT 
    -- ST_Project(pt, distance, 0) move pt by distance north
    -- ST_Project(pt, distance, radians(90)) move pt by distance east
    -- ST_Project(pt, distance, radians(180) move pt by distance south
    -- ST_Project(pt, distance, radians(270)) move pt by distance west
count,
ST_AsEWKT( ST_Project( ST_Project( st_geogfromtext('SRID=4326; POINT(-4.5 63)'), count*10000.0, radians(0)), count*10000.0,  radians(90))),
ST_AsEWKT( ST_Project( ST_Project( st_geogfromtext('SRID=4326; POINT(-4.5 63)'), count*10000.0, radians(90)), count*10000.0, radians(0))),
ST_AsEWKT( ST_Project( ST_Project( st_geomfromtext('SRID=4326; POINT(-4.5 63)'), count*10000.0, radians(0)), count*10000.0,  radians(90))),
ST_AsEWKT( ST_Project( ST_Project( st_geomfromtext('SRID=4326; POINT(-4.5 63)'), count*10000.0, radians(90)), count*10000.0, radians(0)))
FROM generate_series(-5, 5) as count
    ORDER BY count

I feel like one of these geo-systems should allow me to project a point purely along an east-west axis, and I am beginning to feel that it is the geometry type that should be treating the axes separately, and so the st_project function should move in only one axis for azimuth=n*pi/2.

[pramsey]

Use ST_Translate in geometry space. ST_Project(g,a,b) is only defined for geography, so if you feed it a geometry, the system will quietly cast to geography and then run it for you in geography space.

[petermj]

Thank you for your help, I am calculating the metres per degree, at my reference point, and then using that for translating solely along each axis separately. (used to make a grid of tiles)

ST_Distance(ST_Translate(reference_point::geometry, 1.0, 0)::geography, reference_point) as metres_per_degree,

and then using ST_Translate to perform my translation…

    ST_Project(ST_Translate(reference_point::geometry, distance_x_metres/metres_per_degree, 0)::geography, distance_y_metre, radians(180))::geometry,