I have two tables in EPSG:4326, one for home locations and one for site locations. I need to be able to generate the XY offset in meters between the home location, and the site locations. I then need to be able to accurately regenerate the site locations, based on the XY offset from the home location. Currently when I use ST_Project to recreate a site location by using the specified XY offset from the home location. The results are slightly off in the Y direction. Generating lines using the same logic used to calculate the offset, is accurate. My environment is PostgreSQL 9.4.1 PostGIS 2.1.
The SQL script is my test data for this issue. Any help in getting me on the right track would be greatly appreciated. The home and sites test tables are the two tables stored in EPSG:4326. The sites table contains a delta_x, and delta_y field representing the X & Y offset in meters. The offset_lines view contains lines generated using similar logic as was used to calculate the offset. These lines appear accurate and terminate at the exact location as the site. The points_using_offsets view is a view of the site locations generated using the offsets from the home location and the ST_Project command. These locations are all slightly shifted in the Y direction.
CREATE TABLE home ( gid integer, geog geography(Point,4326) ); CREATE TABLE sites ( gid integer, geog geography(Point,4326), delta_x float, delta_y float ); insert into home (gid, geog) values (1, ST_MakePoint(-78.9570, 41.6241)::geography); insert into sites (gid, geog) values (1,ST_MakePoint(-79.1133, 41.7690)::geography); insert into sites (gid, geog) values (2,ST_MakePoint(-78.7726, 41.7801)::geography); insert into sites (gid, geog) values (3,ST_MakePoint(-78.7753, 41.4041)::geography); insert into sites (gid, geog) values (4,ST_MakePoint(-79.4226, 41.5032)::geography); UPDATE sites AS st set delta_x = foo.delta_x, delta_y = foo.delta_y FROM ( SELECT CASE WHEN ST_X(s.geog::geometry) < ST_X(h.geog::geometry) THEN ST_Distance_Spheroid(ST_MakePoint(ST_X(h.geog::geometry), ST_Y(h.geog::geometry) ),ST_MakePoint(ST_X(s.geog::geometry), ST_Y(h.geog::geometry)),'SPHEROID["WGS 84",6378137,298.257223563]')*-1 ELSE ST_Distance_Spheroid(ST_MakePoint(ST_X(h.geog::geometry), ST_Y(h.geog::geometry)),ST_MakePoint(ST_X(s.geog::geometry), ST_Y(h.geog::geometry)),'SPHEROID["WGS 84",6378137,298.257223563]') END AS delta_x, CASE WHEN ST_Y(s.geog::geometry) < ST_Y(h.geog::geometry) THEN ST_Distance_Spheroid(ST_MakePoint(ST_X(s.geog::geometry), ST_Y(s.geog::geometry)),ST_MakePoint(ST_X(s.geog::geometry), ST_Y(h.geog::geometry)),'SPHEROID["WGS 84",6378137,298.257223563]')*-1 ELSE ST_Distance_Spheroid(ST_MakePoint(ST_X(s.geog::geometry), ST_Y(s.geog::geometry)),ST_MakePoint(ST_X(s.geog::geometry), ST_Y(h.geog::geometry)),'SPHEROID["WGS 84",6378137,298.257223563]') END AS delta_y, s.gid FROM home as h, sites as s) as foo where st.gid = foo.gid; CREATE OR REPLACE VIEW offset_lines AS SELECT s.gid as gid, ST_SetSRID(ST_MakeLine(ARRAY[ST_MakePoint(ST_X(h.geog::geometry), ST_Y(h.geog::geometry)),ST_MakePoint(ST_X(s.geog::geometry), ST_Y(h.geog::geometry)),ST_MakePoint(ST_X(s.geog::geometry), ST_Y(s.geog::geometry))]),4326) as geom FROM home as h, sites as s; CREATE OR REPLACE VIEW points_using_offsets AS SELECT CASE WHEN delta_x > 0 and delta_y > 0 THEN ST_Project(ST_Project(h.geog, abs(s.delta_x), radians(90)),abs(s.delta_y),radians(0)) WHEN delta_x > 0 and delta_y < 0 THEN ST_Project(ST_Project(h.geog, abs(s.delta_x), radians(90)),abs(s.delta_y),radians(180)) WHEN delta_x < 0 and delta_y < 0 THEN ST_Project(ST_Project(h.geog, abs(s.delta_x), radians(270)),abs(s.delta_y),radians(180)) WHEN delta_x < 0 and delta_y > 0 THEN ST_Project(ST_Project(h.geog, abs(s.delta_x), radians(270)),abs(s.delta_y),radians(0)) END as geom, s.gid as gid FROM home as h, sites as s;
أكثر...
The SQL script is my test data for this issue. Any help in getting me on the right track would be greatly appreciated. The home and sites test tables are the two tables stored in EPSG:4326. The sites table contains a delta_x, and delta_y field representing the X & Y offset in meters. The offset_lines view contains lines generated using similar logic as was used to calculate the offset. These lines appear accurate and terminate at the exact location as the site. The points_using_offsets view is a view of the site locations generated using the offsets from the home location and the ST_Project command. These locations are all slightly shifted in the Y direction.
CREATE TABLE home ( gid integer, geog geography(Point,4326) ); CREATE TABLE sites ( gid integer, geog geography(Point,4326), delta_x float, delta_y float ); insert into home (gid, geog) values (1, ST_MakePoint(-78.9570, 41.6241)::geography); insert into sites (gid, geog) values (1,ST_MakePoint(-79.1133, 41.7690)::geography); insert into sites (gid, geog) values (2,ST_MakePoint(-78.7726, 41.7801)::geography); insert into sites (gid, geog) values (3,ST_MakePoint(-78.7753, 41.4041)::geography); insert into sites (gid, geog) values (4,ST_MakePoint(-79.4226, 41.5032)::geography); UPDATE sites AS st set delta_x = foo.delta_x, delta_y = foo.delta_y FROM ( SELECT CASE WHEN ST_X(s.geog::geometry) < ST_X(h.geog::geometry) THEN ST_Distance_Spheroid(ST_MakePoint(ST_X(h.geog::geometry), ST_Y(h.geog::geometry) ),ST_MakePoint(ST_X(s.geog::geometry), ST_Y(h.geog::geometry)),'SPHEROID["WGS 84",6378137,298.257223563]')*-1 ELSE ST_Distance_Spheroid(ST_MakePoint(ST_X(h.geog::geometry), ST_Y(h.geog::geometry)),ST_MakePoint(ST_X(s.geog::geometry), ST_Y(h.geog::geometry)),'SPHEROID["WGS 84",6378137,298.257223563]') END AS delta_x, CASE WHEN ST_Y(s.geog::geometry) < ST_Y(h.geog::geometry) THEN ST_Distance_Spheroid(ST_MakePoint(ST_X(s.geog::geometry), ST_Y(s.geog::geometry)),ST_MakePoint(ST_X(s.geog::geometry), ST_Y(h.geog::geometry)),'SPHEROID["WGS 84",6378137,298.257223563]')*-1 ELSE ST_Distance_Spheroid(ST_MakePoint(ST_X(s.geog::geometry), ST_Y(s.geog::geometry)),ST_MakePoint(ST_X(s.geog::geometry), ST_Y(h.geog::geometry)),'SPHEROID["WGS 84",6378137,298.257223563]') END AS delta_y, s.gid FROM home as h, sites as s) as foo where st.gid = foo.gid; CREATE OR REPLACE VIEW offset_lines AS SELECT s.gid as gid, ST_SetSRID(ST_MakeLine(ARRAY[ST_MakePoint(ST_X(h.geog::geometry), ST_Y(h.geog::geometry)),ST_MakePoint(ST_X(s.geog::geometry), ST_Y(h.geog::geometry)),ST_MakePoint(ST_X(s.geog::geometry), ST_Y(s.geog::geometry))]),4326) as geom FROM home as h, sites as s; CREATE OR REPLACE VIEW points_using_offsets AS SELECT CASE WHEN delta_x > 0 and delta_y > 0 THEN ST_Project(ST_Project(h.geog, abs(s.delta_x), radians(90)),abs(s.delta_y),radians(0)) WHEN delta_x > 0 and delta_y < 0 THEN ST_Project(ST_Project(h.geog, abs(s.delta_x), radians(90)),abs(s.delta_y),radians(180)) WHEN delta_x < 0 and delta_y < 0 THEN ST_Project(ST_Project(h.geog, abs(s.delta_x), radians(270)),abs(s.delta_y),radians(180)) WHEN delta_x < 0 and delta_y > 0 THEN ST_Project(ST_Project(h.geog, abs(s.delta_x), radians(270)),abs(s.delta_y),radians(0)) END as geom, s.gid as gid FROM home as h, sites as s;
أكثر...