900913->4326 Transform

[pramsey]

Originally reported on postgis-users against 1.3.6, strange results are also available on 1.5:

geog=# select postgis_full_version();
                                                    postgis_full_version                                                    
----------------------------------------------------------------------------------------------------------------------------
 POSTGIS="1.5.0SVN" GEOS="3.2.0-CAPI-1.6.0" PROJ="Rel. 4.6.1, 21 August 2008" USE_STATS (procs from 1.5 r4836 need upgrade)
(1 row)

geog=# SELECT ASTEXT(ST_Transform(ST_SetSRID(ST_Point(-633510.090428,7506727.67383),900913),4326));
                  astext                   
-------------------------------------------
 POINT(-5.69091796875412 55.7394816986905)
(1 row)

geog=# SELECT ASTEXT(ST_Transform(ST_SetSRID(ST_Point(-305748.113141,7834489.65112),900913),4326));
ERROR:  transform: couldn't project point (-305748 7.83449e+06 0): failed to load NAD27-83 correction file (-38)
HINT:  PostGIS was unable to transform the point because either no grid shift files were found, or the point does not lie within the range for which the grid shift is defined. Refer to the ST_Transform() section of the PostGIS manual for details on how to configure PostGIS to alter this behaviour.
geog=# SELECT ASTEXT(ST_Transform(ST_SetSRID(ST_Point(-305748.113141,7834489.65112),900913),4269));
                  astext                   
-------------------------------------------
 POINT(-2.74658203125265 57.3620898624176)
(1 row)

geog=# SELECT ASTEXT(ST_Transform(ST_SetSRID(ST_Point(-305748.113141,7834489.65112),900913),4267));
                  astext                   
-------------------------------------------
 POINT(-2.74658203125265 57.3620898615614)
(1 row)

Note that for some reason we get NAD27/83 grid failure in the 900913->4326 case but *not* in the transforms into NAD83 and NAD27!

[pramsey]

Actually, note in the example above the first call to transform works, the second fails. This is also exhibited in the postgis-users case, though the failure in his case is wrong answers.

​http://postgis.refractions.net/pipermail/postgis-users/2009-November/025234.html

[mcayland]

Hmmm I can't recreate this on SVN trunk or 1.4 here (Debian Lenny, amd64):

pg83@zeno:~$ psql -d trunk Welcome to psql 8.3.7, the PostgreSQL interactive terminal.

Type: \copyright for distribution terms

\h for help with SQL commands \? for help with psql commands \g or terminate with semicolon to execute query \q to quit

trunk=# SELECT ASTEXT(ST_Transform(ST_SetSRID(ST_Point(-633510.090428,7506727.67383),900913),4326));

astext

---

POINT(-5.69091796875412 55.7394816986905)

(1 row)

trunk=# SELECT ASTEXT(ST_Transform(ST_SetSRID(ST_Point(-633510.090428,7506727.67383),900913),4326));

astext

---

POINT(-5.69091796875412 55.7394816986905)

(1 row)

trunk=# SELECT ASTEXT(ST_Transform(ST_SetSRID(ST_Point(-633510.090428,7506727.67383),900913),4326));

astext

---

POINT(-5.69091796875412 55.7394816986905)

(1 row)

trunk=# SELECT ASTEXT(ST_Transform(ST_SetSRID(ST_Point(-633510.090428,7506727.67383),900913),4326));

astext

---

POINT(-5.69091796875412 55.7394816986905)

(1 row)

This is using a PROJ.4 I've compiled myself. Can anyone else recreate this?

[nicklas]

I get the same error with:

select ASTEXT(ST_Transform(ST_SetSRID(ST_Point(-305748.113141,7834489.65112),900913),4267));

I'm not sure if Paul meant that was when he gets the error, but it is for me. so it is not 4326 but 4327 for me /Nicklas

[mcayland]

Here is the output from another session:

trunk=# SELECT ASTEXT(ST_Transform(ST_SetSRID(ST_Point(-633510.090428,7506727.67383),900913),4326));

astext

---

POINT(-5.69091796875412 55.7394816986905)

(1 row)

trunk=# SELECT ASTEXT(ST_Transform(ST_SetSRID(ST_Point(-305748.113141,7834489.65112),900913),4326));

astext

---

POINT(-2.74658203125265 57.3620898615614)

(1 row)

trunk=# SELECT ASTEXT(ST_Transform(ST_SetSRID(ST_Point(-305748.113141,7834489.65112),900913),4269));

astext

---

POINT(-2.74658203125265 57.3620898624176)

(1 row)

trunk=# SELECT ASTEXT(ST_Transform(ST_SetSRID(ST_Point(-305748.113141,7834489.65112),900913),4267)); ERROR: transform: couldn't project point (-305748 7.83449e+06 0): failed to load NAD27-83 correction file (-38) HINT: PostGIS was unable to transform the point because either no grid shift files were found, or the point does not lie within the range for which the grid shift is defined. Refer to the ST_Transform() section of the PostGIS manual for details on how to configure PostGIS to alter this behaviour.

So is the bug here that 4267 should work or perhaps the point is genuinely out of range?

[mloskot]

I can not reproduce Paul's problem and I get similar output to Mark's:

mloskot@vb-ubuntu910-x64:~$ psql ticket318
psql (8.4.1)
Type "help" for help.

ticket318=# SELECT postgis_full_version();
                                         postgis_full_version                                          
-------------------------------------------------------------------------------------------------------
 POSTGIS="1.5.0SVN" GEOS="3.3.0-CAPI-1.6.1" PROJ="Rel. 4.6.1, 21 August 2008" LIBXML="2.7.5" USE_STATS
(1 row)

ticket318=# SELECT ASTEXT(ST_Transform(ST_SetSRID(ST_Point(-633510.090428,7506727.67383),900913),4326));
                  astext                   
-------------------------------------------
 POINT(-5.69091796875412 55.7394816986905)
(1 row)

ticket318=# SELECT ASTEXT(ST_Transform(ST_SetSRID(ST_Point(-305748.113141,7834489.65112),900913),4326));
                  astext                   
-------------------------------------------
 POINT(-2.74658203125265 57.3620898615614)
(1 row)

ticket318=# SELECT ASTEXT(ST_Transform(ST_SetSRID(ST_Point(-305748.113141,7834489.65112),900913),4269));
                  astext                   
-------------------------------------------
 POINT(-2.74658203125265 57.3620898624176)
(1 row)

ticket318=# SELECT ASTEXT(ST_Transform(ST_SetSRID(ST_Point(-305748.113141,7834489.65112),900913),4267));
ERROR:  transform: couldn't project point (-305748 7.83449e+06 0): failed to load NAD27-83 correction file (-38)
HINT:  PostGIS was unable to transform the point because either no grid shift files were found, or the point does not lie within the range for which the grid shift is defined. Refer to the ST_Transform() section of the PostGIS manual for details on how to configure PostGIS to alter this behaviour.
ticket318=# 

I use GEOS and PostGIS built from SVN trunk today morning on Linux Ubuntu 9.10 (64-bit)

[pramsey]

OK, seeing as this is a rarer condition than I thought, I'm rolling forward to the next point release. Fortunately I can see it so I can try and fix it, but it's not a blocker.

[arthurnederlof]

Hi, I have this problem as wel, running postgis 1.4 on the OSM data. ST_Transform(ST_SetSRID(ST_MakePoint("+lon+","+lat+"),4326),900913) is not a problem ST_TransForm(way,4326) goes wrong for all european points (error as above) Using srid=4267 seems to work ok, not sure if the resulting coor. are correct.

Is there a workaround/patch thinkable to solve this in 1.4?

regards, Arthur

[mcayland]

Hi Arthur,

If you look at the thread above, the issue is that so far only Paul can reproduce the error. You can help us by giving as much information about your system, OS, PostgreSQL version, GEOS version and PROJ version as possible, as well as complete test case.

HTH,

Mark.

[arthurnederlof]

Replying to mcayland:

Hi Arthur,

If you look at the thread above, the issue is that so far only Paul can reproduce the error. You can help us by giving as much information about your system, OS, PostgreSQL version, GEOS version and PROJ version as possible, as well as complete test case.

HTH,

Mark.

Mark, No problem, the probem is very consistent. below the enviromnet, testoutput and a small demoprogram. The used testcoordinates are from Cologne, germany, but all tested coordinates in Europe give the same problem. If you need more, please ask. Arthur

system: GNU bash, version 3.2.48(1)-release (x86_64-pc-linux-gnu) postgresql:"PostgreSQL 8.3.9 on x86_64-pc-linux-gnu, compiled by GCC gcc-4.3.real (Ubuntu 4.3.3-5ubuntu4) 4.3.3" ​postgis:"POSTGIS="1.4.1" GEOS="3.0.0-CAPI-1.4.1" PROJ="Rel. 4.6.0, 21 Dec 2007" USE_STATS"

Output testprogram: 15:18:13,445 INFO test::15 - start testrun

15:18:13,449 DEBUG DbfMaintenance::483 - Access to postGreSql on jdbc:postgresql://10.13.54.8:5432/osm, user= arthur 15:18:13,581 DEBUG test::54 - postgreql version = PostgreSQL 8.3.9 on x86_64-pc-linux-gnu, compiled by GCC gcc-4.3.real (Ubuntu 4.3.3-5ubuntu4) 4.3.3 15:18:13,582 DEBUG test::55 - postgis version = POSTGIS="1.4.1" GEOS="3.0.0-CAPI-1.4.1" PROJ="Rel. 4.6.0, 21 Dec 2007" USE_STATS 15:18:13,582 DEBUG test::60 - sql:select osm_id, name from planet_osm_polygon where (ST_Contains(way,ST_Transform(ST_SetSRID(ST_MakePoint(6.95,50.9333333),4326),900913))) limit 1 15:18:13,591 DEBUG test::63 - name = Köln, Stadt 15:18:13,592 DEBUG test::68 - sql:select osm_id, name, ST_AsGeoJSON( ST_TransForm(ST_setSRID(St_Box2D(way)::box2d,900913),4326),6 ) as box from planet_osm_polygon where (ST_Contains(way,ST_Transform(ST_SetSRID(ST_MakePoint(6.95,50.9333333),4326),900913))) limit 1 15:18:13,620 DEBUG test::75 - problem with sql : ERROR: transform: couldn't project point (753919 6.59133e+06 0): failed to load NAD27-83 correction file (-38)

Hint: PostGIS was unable to transform the point because either no grid shift files were found, or the point does not lie within the range for which the grid shift is defined. Refer to the ST_Transform() section of the PostGIS manual for details on how to configure PostGIS to alter this behaviour.

15:18:13,620 DEBUG DbfMaintenance::507 - disconnected to database. 15:18:13,621 INFO test::17 - Run postgistest OK

The testrountine itself:

public static String demontratePostgisProblem(double lat, double lon){

Connection conn=DbfMaintenance.getPgDbConn(); String sql="";

try {

Statement stmt= conn.createStatement();

stmt = conn.createStatement(); ResultSet rset; rset = stmt.executeQuery ("select version() as postresqlversion, postgis_full_version() as postgisversion");

while (rset.next()) {

logger.debug("postgreql version = "+rset.getString("postresqlversion")); logger.debug("postgis version = "+rset.getString("postgisversion"));

}

this does work

sql = "select osm_id, name from planet_osm_polygon"+

" where (ST_Contains(way,ST_Transform(ST_SetSRID(ST_MakePoint("+lon+","+lat+"),4326),900913))) limit 1";

logger.debug("sql:"+sql); rset = stmt.executeQuery (sql); while (rset.next()) {

logger.debug("name = "+rset.getString("name"));

} this one gives the error for 4326, 4267 works ok sql = "select osm_id, name, ST_AsGeoJSON( ST_TransForm(ST_setSRID(St_Box2D(way)::box2d,900913),4326),6 ) as box from planet_osm_polygon"+

" where (ST_Contains(way,ST_Transform(ST_SetSRID(ST_MakePoint("+lon+","+lat+"),4326),900913))) limit 1";

logger.debug("sql:"+sql); rset = stmt.executeQuery (sql); while (rset.next()) {

logger.debug("name = "+rset.getString("name")); logger.debug("box = "+rset.getString("box"));

}

} catch (SQLException e) {

logger.debug("problem with sql : "+e.getMessage());

} DbfMaintenance.closeDbConn(conn); return"OK";

}

[arthurnederlof]

Hallo, a more trivial example (osm form open street map, worldcountries has geometries in 4326 srid format: this works: select world.gid, osm.admin_level from planet_osm_polygon as osm,worldcountries as world where (osm.admin_level='4') AND (ST_Contains(ST_Transform(world.the_geom,900913),ST_Centroid(osm.way))) limit 10

this one does not work: select world.gid, osm.admin_level from planet_osm_polygon as osm,worldcountries as world where (osm.admin_level='4') AND (ST_Contains(world.the_geom,ST_Centroid(ST_Transform(osm.way,4326)))) limit 10

ERROR: transform: couldn't project point (5.1755e+06 6.1799e+06 0): failed to load NAD27-83 correction file (-38) HINT: PostGIS was unable to transform the point because either no grid shift files were found, or the point does not lie within the range for which the grid shift is defined. Refer to the ST_Transform() section of the PostGIS manual for details on how to configure PostGIS to alter this behaviour.

[mcayland]

Unfortunately, without a database dump of the relevant tables then it is impossible for me to re-create this here. Can you extract a small set of geometries that will re-create the bug into a separate table and then provide a pg_dump file that I can use for testing?

[arthurnederlof]

Replying to mcayland:

Unfortunately, without a database dump of the relevant tables then it is impossible for me to re-create this here. Can you extract a small set of geometries that will re-create the bug into a separate table and then provide a pg_dump file that I can use for testing?

I could give you access to my database if you want to. No special or private tables there. Arthur

[robe]

Guys,

I can produce this error on my PostGIS 1.5SVN windows install

SELECT ST_ASTEXT(ST_Transform(ST_SetSRID(ST_Point(-305748.113141,7834489.65112),900913),4267));


Yields:
couldn't project point (-305748 7.83449e+006 0): failed to load NAD27-83 correction file (-38)
HINT:  PostGIS was unable to transform the point because either no grid shift files were found, or the point does not lie within the range for which the grid shift is defined. Refer to the ST_Transform() section of the PostGIS manual for details on how to configure PostGIS to alter this behaviour.

This could possibly be the same problem we ran across in #387 that the proj4text is simply different, or possibly I'm using an out of date set of datum shift files.

I tried my 1.3.6 install, and that doesn't even have this entry. It is my understanding that this is a defunct srid put in a while ago before there was an official google mercator EPSG SRID, so its highly likely everyones is slightly different unless they reloaded their spatial_ref_sys.

Mine reads:

SELECT proj4text from spatial_ref_sys where srid = 900913;


+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

and I'm running proj 4.6.1

[robe]

Slight correction -- I see the one that fails for me is the NAD27 and not the other one that was complained about. But it would be still interesting to find out if we are all working against the same proj4text.

For record:
SELECT ST_AsText(ST_Transform(ST_SetSRID(ST_Point(-305748.113141,7834489.65112),900913),4326));

Yields:
POINT(-2.74658203125265 57.3620898615614)

[filbertkm]

I had the same problem. I upgraded from PostGIS 1.3.6 to the 1.5 svn version. I used spatial_ref_sys and the postgis upgrade script. I also used a 900913.sql script of my own (can't remember where I got it).

Since PostGIS 1.3.6, the Google Spherical Mercator projection has been assigned an EPSG code:

​http://spatialreference.org/ref/epsg/3785/

My database still had all the tables assigned as 900913. To resolve the problem, I dropped the 900913 entry from my spatial_ref_sys table:

Then I took the PostGIS spatial_ref_sys INSERT statement for 3785 from the spatial reference site. I updated my 900913.sql statement with this one, but changing the srid and auth_srid to 900913. I then readded 900913 to PostGIS and everything works great now.

[arthurnederlof]

Replying to arthurnederlof:

Replying to mcayland:

Unfortunately, without a database dump of the relevant tables then it is impossible for me to re-create this here. Can you extract a small set of geometries that will re-create the bug into a separate table and then provide a pg_dump file that I can use for testing?

I could give you access to my database if you want to. No special or private tables there. Arthur

Ok, I exctracted thes sql without using a database, reproducing the error. Changing the srid as descibed by filbertkm does not give any result on my side.

select ST_AsGeoJSON( ST_TransForm(ST_setSRID(St_Box2D(ST_Transform(ST_SetSRID(ST_MakePoint(6.95,50.9333333),4326),900913))::box2d,900913),4326),6 )

hop this helps understanding what goes on, Arthur

[mcayland]

Unfortunately this still works for me:

trunk=# select ST_AsGeoJSON( ST_TransForm(ST_setSRID(St_Box2D(ST_Transform(ST_SetSRID(ST_MakePoint(6.95,50.9333333),4326),900913))::box2d,900913),4326),6 );

st_asgeojson

---

{"type":"Polygon","coordinates":[6.95,50.933333],[6.95,50.933335],[6.95,50.933335],[6.95,50.933333],[6.95,50.933333]}

(1 row)

trunk=# select * from spatial_ref_sys where srid = 900913;

srid | auth_name | auth_srid | srtext | proj4text

---

900913 | spatialreferencing.org | 900913 | PROJCS["Popular Visualisation CRS / Mercator (deprecated)",GEOGCS["Popular Visualisation CRS",DATUM["Popular_Visualisation_Datum",SPHEROID["Popular Visualisation Sphere",6378137,0,AUTHORITY["EPSG","7059"]],TOWGS84[0,0,0,0,0,0,0],AUTHORITY["EPSG","6055"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.01745329251994328,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4055"]],UNIT["metre",1,AUTHORITY["EPSG","9001"]],PROJECTIONMercator_1SP,PARAMETER["central_meridian",0],PARAMETER["scale_factor",1],PARAMETER["false_easting",0],PARAMETER["false_northing",0],AUTHORITY["EPSG","3785"],AXIS["X",EAST],AXIS["Y",NORTH]] | +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

(1 row)

[mcayland]

Hi Arthur,

Where did you get your PROJ.4 installation from? Did you compile it yourself or install from packages?

ATB,

Mark.

[arthurnederlof]

Replying to mcayland:

Hi Arthur,

Where did you get your PROJ.4 installation from? Did you compile it yourself or install from packages?

ATB,

Mark.

The installation is from the Ubuntu repository, I did not build it myself. project 4.6.0-2 Arthur.

[pbalomiri]

Hi, I got this bug too. I am working on os/x The lib versions are:

template_postgis=# select * from postgis_version();
            postgis_version            
---------------------------------------
 1.5 USE_GEOS=1 USE_PROJ=1 USE_STATS=1

template_postgis=# select * from postgis_proj_version();
     postgis_proj_version      
-------------------------------
 Rel. 4.7.1, 23 September 2009

The Error consistently appears when using SRID 4326:

wien_python=# select st_transform(st_transform('srid=4326;POINT(1 1)'::geometry, 900913), 4326);
ERROR:  transform: couldn't project point (111319 111325 0): failed to load NAD27-83 correction file (-38)
HINT:  PostGIS was unable to transform the point because either no grid shift files were found, or the point does not lie within the range for which the grid shift is defined. Refer to the ST_Transform() section of the PostGIS manual for details on how to configure PostGIS to alter this behaviour.

For SRID 4267, however it works

wien_python=# select ST_Transform('SRID=900913;POINT(653103 6.63036e+06 0)'::geometry, 4267);

st_transform

---

01010000A0AB10000037CD83F1BA771740E31A12928E8649400000000000000000

(1 row)

[pbalomiri]

Sorry for the formatting err:

For SRID 4267, however it works

 wien_python=# select ST_Transform('SRID=900913;POINT(653103 6.63036e+06 0)'::geometry, 4267);
                             st_transform                            
 --------------------------------------------------------------------
  01010000A0AB10000037CD83F1BA771740E31A12928E8649400000000000000000
 (1 row)
 {{{
> 
}}}

[realityexists]

I'm also getting this error on PostGIS 1.5.2 on Windows, installed from the installer using Application Stack builder.

SELECT Transform(ST_GeomFromEWKT('SRID=4267;POINT(-79.867611 9.3553472)'), 4326)

ERROR: transform: couldn't project point (-79.8676 9.35535 0): failed to load NAD27-83 correction file (-38)

However, this succeeds:

SELECT Transform(Transform(ST_GeomFromEWKT('SRID=4267;POINT(-79.867611 9.3553472)'), 3395), 4326)

PostGIS full version: "POSTGIS="1.5.2" GEOS="3.2.2-CAPI-1.6.2" PROJ="Rel. 4.6.1, 21 August 2008" LIBXML="2.7.6" USE_STATS" PostgreSQL version 9.0.4 on Windows, 32-bit

[kbayram]

I came across this bug and Google brought me here. I've managed to fix the problem in my case and it turns out it was a extremely trivial problem.

If you do a "select * from spatial_ref_sys where srid = 900913;" in psql you may notice that the text in the srtext and proj4text has line feeds. It turns out that proj4 does not like these and anything after the first line in the proj4text field is ignored, namely the important nadgrids=@null which causes the error.

Removing the line feeds solved the problem for me.