[GIS] Coordinate values are out of range [-180 -90, 180 90] for GEOGRAPHY type

osm2pgsqlpostgissrid

After import some data from OSM file into Postgre, the coordinates numbers are too big and having trouble to cast to Geography:

Coordinate values are out of range [-180 -90, 180 90] for GEOGRAPHY type

In my table:

select ST_AsText( the_geom )  from routes limit 4

"LINESTRING(-4902130.86 -2265840.91,-4902083.04 -2265907.16)"
"LINESTRING(-4903296.22 -2266295.18,-4903363.69 -2266333.39)"
"LINESTRING(-5705594.67 -3518944.1,-5705523.22 -3518937.37)"
"LINESTRING(-5484056.84 -1884629.16,-5484093.47 -1884510.24)"

I think the numbers should be between (-90 / 90, -180 / 180) …

Info:

1: select ST_SRID( the_geom ) returns 900913 for all lines.

2: select UpdateGeometrySRID('public', 'routes', 'the_geom', 4326) does not converts the numbers, but SRID now shows 4326.

Best Answer

From @Jakub Kania tip :

ALTER TABLE routes ADD temp geometry;
select UpdateGeometrySRID('public', 'routes', 'temp', 4326) ;
update routes set temp = ST_Transform(the_geom,4326);
select UpdateGeometrySRID('public', 'routes', 'the_geom', 4326) ;
update routes set the_geom = temp;
alter table routes drop column temp;

Related Solutions

[GIS] Unions of polygons (Postgis ST_UNION for Geography type)

You can cast to geometry and do the calculation. If your data goes over the poles or date-line that could be problematic. Also if you have long edges, in which the great-circle interpolation line differs from a simple cartesian interpolation. In those cases you might have to create your own aggregate function that re-projects to something reasonable before running the operation.

SELECT ST_Transform(
  ST_Union(ST_Transform(geog::geometry, 26910)),
  4326)::geography 
FROM the_table

[GIS] Getting ST_ClosestPoint for Geography type

Sadly there isn't a geographic version of ST_ClosestPoint, so you will have to write your own function. There are two ways of calculating the nearest point of a great circle: spherical trigonometry, or 3D vector algebra. Luckily for you I have just written such a function for the latter method; I've not attempted the former because my spherical trig is pretty poor.

I've written a PL/Python function that you can add to your database, although you will need to install the Python libraries and PL/Python wrappers, which if you're using Linux should be simply to install the postgresql-python-x.x package (where x.x is the version of PostgreSQL you're using).

This function isn't necessarily the most efficient, it includes no bounds or error checking, and it only uses a sphere, but it certainly works well enough with the data I've tried it with. Feel free to modify accordingly.

First I define my own simple long/lat type:

CREATE TYPE lon_lat AS (lon float, lat float);

The function looks like this:

CREATE OR REPLACE FUNCTION geog_nearest_point(lp1 lon_lat, lp2 lon_lat, p lon_lat) RETURNS lon_lat AS
$$
import math
def sph2cart(lon, lat, r):
    return (r * math.cos(lat) * math.cos(lon),
            r * math.cos(lat) * math.sin(lon),
            r * math.sin(lat))

def cart2sph(x, y, z):
    return (math.atan2(y,x),
            math.atan2(z, math.sqrt(x * x + y * y)),
            math.sqrt(x * x + y * y + z * z))

def cross_prod(v1, v2):
    return (v1[1] * v2[2] - v1[2] * v2[1],
            v1[2] * v2[0] - v1[0] * v2[2],
            v1[0] * v2[1] - v1[1] * v2[0])

lv1 = sph2cart(math.radians(lp1['lon']), math.radians(lp1['lat']), 1.0)
lv2 = sph2cart(math.radians(lp2['lon']), math.radians(lp2['lat']), 1.0)
pv = sph2cart(math.radians(p['lon']), math.radians(p['lat']), 1.0)

f = cross_prod(lv1, lv2)
g = cross_prod(pv, f)
h = cross_prod(f, g)

nearest_point = cart2sph(h[0], h[1], h[2])
return (math.degrees(nearest_point[0]), math.degrees(nearest_point[1]))
$$
LANGUAGE 'plpythonu' IMMUTABLE;

And giving it the data you provided gives me this:

geog_db=# SELECT geog_nearest_point((41, 15)::lon_lat, (40, 5)::lon_lat, (41, 12)::lon_lat);
      geog_nearest_point
-------------------------------
 (40.6960040707,12.0295700818)
(1 row)

It works by finding the point f on the sphere which is normal to the plane that describes the great circle. This point is unique for each great circle.
Next, it uses f and the input point p to define a plane perpendicular to the input great circle, and determines its unique normal point g. The vectors from the centre of the sphere to f and g define a plane that is perpendicular to the previous two planes.
Finding the normal of this third plane gives us the point at which the two great circles intersect.
Finally, this Cartesian coordinate is converted back to latitude and longitude.

If anyone has a spherical trig method, which I think will be slightly more efficient, then I'd love to see it.

Best Answer

Related Solutions

[GIS] Unions of polygons (Postgis ST_UNION for Geography type)

[GIS] Getting ST_ClosestPoint for Geography type

Related Question