[postgis-tickets] [PostGIS] #3924: ST_Project does not project perfectly east-west

PostGIS trac at osgeo.org
Wed Nov 8 02:07:25 PST 2017 

#3924: ST_Project does not project perfectly east-west
---------------------+---------------------------
 Reporter:  petermj  |      Owner:  pramsey
     Type:  defect   |     Status:  new
 Priority:  medium   |  Milestone:  PostGIS 2.4.2
Component:  postgis  |    Version:  2.4.x
 Keywords:           |
---------------------+---------------------------
 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

--
Ticket URL: <https://trac.osgeo.org/postgis/ticket/3924>
PostGIS <http://trac.osgeo.org/postgis/>
The PostGIS Trac is used for bug, enhancement & task tracking, a user and developer wiki, and a view into the subversion code repository of PostGIS project.

More information about the postgis-tickets mailing list