[OSM-dev] Help needed: postgis database transform()

[Richard Duivenvoorde]

Hi Joerg,
apparently the schema you use is not properly 'spatial enabled'.

If you have a postgresql database, and installed postgis properly, you 
either have to 'spatial enable' a database by using the right 
'gis'-template, OR you run two big sql file which give your database:
- a lot of spatial functions
- a table called 'geometry_columns' (to 'registre' a geometry column: 
give it the right srid/proj etc)
- a table called 'spatial_ref_sys' (which gives postgis the proper 
parameters for the proj-functions to do reprojection etc).

So if your spatial_ref_sys is empty, you should run a script which is 
called spatial_ref_sys.sql (here with me it's in 
/usr/share/postgresql-8.2-postgis/spatial_ref_sys.sql)
The other script by the way is in the same directory and called 
'lwpostgis.sql'.

select * from spatial_ref_sys where srid = 4326
gives me something like:
4326;"EPSG";4326;"GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 
84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],TOWGS84[0,0,0,0,0,0,0],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.01745329251994328,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]]";"+proj=longlat 
+ellps=WGS84 +datum=WGS84 +no_defs "

So IF you want to transform to or from 4326, this record should be 
available....

But besides this, if the import of your osm-data went ok, there is also 
an srid in the geometry itself defined:
select srid(way) from planet_osm_point limit 5
you can also view this with:
select AsEWKT(way) from planet_osm_point limit 5

Now the clue is that to use spatial functions, both need a or the same srid.

By the way:
900913 is normally NOT in the spatial_ref_sys because it's not an 
'official epgs' one. But you can put it in there, see:
http://spatialreference.org/ref/user/6/

INSERT into spatial_ref_sys (srid, auth_name, auth_srid, proj4text, 
srtext) values ( 900913, 'spatialreference.org', 6, '+proj=merc 
+a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 
+units=m +nadgrids=@null +wktext  +no_defs', 
'PROJCS["unnamed",GEOGCS["unnamed 
ellipse",DATUM["unknown",SPHEROID["unnamed",6378137,0]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433]],PROJECTION["Mercator_2SP"],PARAMETER["standard_parallel_1",0],PARAMETER["central_meridian",0],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["Meter",1],EXTENSION["PROJ4","+proj=merc 
+a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 
+units=m +nadgrids=@null +wktext  +no_defs"]]');

now you can do:

select osm_id,
asEWKT(way),
asEWKT(transform(way,900913)),
asEWKT(transform(way,4326))
from planet_osm_point limit 1

Hope this helps,

Regards,

Richard Duivenvoorde
Joerg Ostertag (OSM Tettnang/Germany) wrote:
> Probably my problem boils down to :
> 	Which SRID can I use to get data into the coordinate system used inside the
> 	planet.osm?
> 
> But here in Detail:
> I need a little help in reading from the mapnik posgis database. I plan to 
> expand the osmtrackfilter.pl to compare against the local data already 
> existing in my mapnik-postgis-database. So I tried to get all street-segments 
> inside a bounding box. For this i tried to select from the database. But 
> there I have the problem that the coordintesystems are different and I have 
> no clue how to convert between them.
> I want to operate with coordinates compatible to the planet.osm File.
> Something like 
> 	lat="48.8046469" lon="9.0401476"
> 
> So i thought i'll start with
> 	echo "select osm_id,asText(way) from planet_osm_point limit 5;" | psql gis
> Result:
> 	  osm_id   |                  astext                   
> 	-----------+-------------------------------------------
> 	 293106148 | POINT(-2876147.15643226 4560658.25114169)
> 	 293105894 | POINT(-2874584.97662211 4563768.28055896)
> 	 293107003 | POINT(-2871694.69962705 4559831.77530905)
> 	 292977203 | POINT(-2871581.82166339 4556767.11400322)
> 	 292976398 | POINT(-2870606.29556972 4556208.15857941)
> 	(5 rows)
> 
> Well at least I get some data from the database. 
> Then the next step is to get transformed Data in the coordinate system I 
> desire ... I did remember WGS-84 was SRID 4326, but either I'm wrong or 
> whatever. Doing the following 
> 	echo "select osm_id,asText(transform(way,4326)) from planet_osm_point limit
> 	         5;" | psql gis
> Results is an:
> 	ERROR:  AddToPROJ4SRSCache: Cannot find SRID (4326) in spatial_ref_sys
> 
> 
> Well I thought I'll have a look into this in "spatial_ref_sys".... but ....
> 	echo "select * from spatial_ref_sys;" | psql gis
> Result:
> 	 srid | auth_name | auth_srid | srtext | proj4text 
> 	------+-----------+-----------+--------+-----------
> 	(0 rows)
> 
> Then I found SRID 900913 inside osm2pgsql and tried:
> 	echo "select osm_id,asText(transform(way,900913)) from planet_osm_point limit
> 	 5;" | psql gis
> Result:
> 	  osm_id   |                  astext                   
> 	-----------+-------------------------------------------
> 	 293106148 | POINT(-2876147.15643226 4560658.25114169)
> 	 293105894 | POINT(-2874584.97662211 4563768.28055896)
> 	 293107003 | POINT(-2871694.69962705 4559831.77530905)
> 	 292977203 | POINT(-2871581.82166339 4556767.11400322)
> 	 292976398 | POINT(-2870606.29556972 4556208.15857941)
> 	(5 rows)
> 
> Well, this seems to be the same SRID which is used inside the Database. 
> 
> So my Question is:
> 	Which SRID can I use to get into the coordinate system I want?
> 
> 
> Jörg
> 
> http://www.ostertag.name/