3

I'm trying to define a default point setting for OSM layers in my application and I want users to be able to define the default latitude and longitude in EPSG:4326. What I'm doing right now is to use GEOS Point.transform to turn the EPSG:4326 input into EPSG:3857:

_default_point = Point(
    getattr(settings, 'DEFAULT_LON', None) or 0,
    getattr(settings, 'DEFAULT_LAT', None) or 0,
    srid=4326
)
_default_point.transform(3857)
default_lat, default_lon = _default_point.get_coords()

Then, i.e. for (-2.41544723477347,43.2140883725249) I'm getting (-268966.0382412352, 4810570.312724996). That, as for what I can see in the map, is wrong.

Now, when trying over psql

SELECT ST_AsText(ST_Transform(ST_SetSRID(ST_Point(-2.41544723477347,43.2140883725249),4326),3857));
                 st_astext                 
-------------------------------------------
 POINT(-268886.356213004 5344615.26048709)

Which is right.

What is going on then? Why are the results different?

EDIT:

select srid, proj4text from spatial_ref_sys where srid = 3857 or srid = 900913 or srid = 4326 yields:

   4326 | +proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs 
   3857 | +proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +units=m +k=1.0 +nadgrids=@null +no_defs
 900913 | +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 +no_defs

2 Answers 2

1

Passing the arguments as (DEFAULT_LAT, DEFAULT_LON) made me realize that those settings where in the wrong order in the configuration file and had its values transposed. Also, after transform, I was using the wrong order in unpack, and that way I was transposing the wrong result.

The right code is

_default_point = Point(
    getattr(settings, 'DEFAULT_LON', None) or 0,
    getattr(settings, 'DEFAULT_LAT', None) or 0,
    srid=4326
)
_default_point.transform(3857)

default_lon, default_lat = _default_point.coords
0

I haven't tried using GEOS directly so not sure how it looks up the proj details. I think it does use proj4 as well to transform.

Things to check: 1) In spatial_ref_sys table of postgis look up the value of proj4text field for srid = 3857

2) If there is a geos function that allows you to pass that in then at elast you can more or less compare apples and apples (or if GEOS can spit out it's corresponding proj4text for 3857.

3) The other thing that affects these are the grid shift files. I forget if 3857 is one affected by it or not, but you'll want to verify GEOS is using the same proj data set (include shift files) as PostGIS.

4) The actual proj library used. There have been changes along the way .

SELECT postgis_full_version();

Will list the project library version being used. Not sure how you check this with GEOS. I don't even recall GEOS being compilable with a proj version.

1
  • It turned out to be a problem with the order of the arguments and the order of Point.get_coords and Point.coords being wrongly assumed. Thanks anyway!
    – sebasmagri
    Commented Nov 18, 2012 at 1:16

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.