headers
Havard Tveite | 20 Dec 16:10

transform() does not do work for NGO1948 zone III -> WGS84 longlat (epsg:27393->epsg:4326)

I have a table with coordinates in NGO1948, Zone 3 (epgs:27393).
I would like to do transformations on the fly to WGS84 longlat
(epsg:4326).

"Reprojection" into WGS84 longlat (epsg:4326) works OK using
Mapserver (fetching the data from postgres using the default
epsg:27393 and projecting to epsg:4326 in Mapserver). No
noticable errors in the output map, so I guess proj.4 is OK.

When using transform() in postgis, however, the longitude of origin
(+lon_0=-10.72291666666667) seems to be ignored, and there seems to
be no datum shift performed
(+towgs84=278.3,93,474.5,7.889,0.05,-6.61,6.21).
Longitude of origin (+lat_0=58) seems to be used correctly.

The following queries "demonstrate" the problem:

select AsText(the_geom) from mndplngo;
 > POINT(2836.394208896 185526.991186176)

select AsText(transform(the_geom,4326)) from mndplngo;
 > POINT(0.0464201635443672 59.6661916850467)

(the longitude of origin is 10.72291666666667, so a longitude
value of about 10.77 should be expected)

select * from geometry_columns;
 >  f_table_catalog | f_table_schema | f_table_name | f_geometry_column | coord_dimension | srid  | type
 > -----------------+----------------+--------------+-------------------+-----------------+-------+-------
 >                  | public         | mndplngo     | the_geom          |               2 | 27393 | POINT
(Continue reading)

Permalink | Reply |
headers
Mark Cave-Ayland | 20 Dec 17:46

RE: transform() does not do work for NGO1948 zone III -> WGS84 longlat (epsg:27393->epsg:4326)


> -----Original Message-----
> From: postgis-users-bounces <at> postgis.refractions.net [mailto:postgis-users-
> bounces <at> postgis.refractions.net] On Behalf Of Havard Tveite
> Sent: 20 December 2005 15:12
> To: postgis-users <at> postgis.refractions.net
> Subject: [postgis-users] transform() does not do work for NGO1948 zone III
> -> WGS84 longlat (epsg:27393->epsg:4326)

(cut)

> When using transform() in postgis, however, the longitude of origin
> (+lon_0=-10.72291666666667) seems to be ignored, and there seems to
> be no datum shift performed
> (+towgs84=278.3,93,474.5,7.889,0.05,-6.61,6.21).
> Longitude of origin (+lat_0=58) seems to be used correctly.

Hi Havard,

Actually, looking through the source lwgeom_transform.c I see the following
comment around line 500:

//this is *exactly* the same as PROJ.4's pj_transform(), but it doesn't
//do the datum shift.

So does this mean that the PostGIS transform() function won't do datum
shifts? Can anyone with a bit more PROJ4 knowledge say if this is correct,
and if so, why PostGIS can't support it?

Kind regards,
(Continue reading)

Permalink | Reply |
headers
strk | 20 Dec 22:19
Favicon

Re: transform() does not do work for NGO1948 zone III -> WGS84 longlat (epsg:27393->epsg:4326)

On Tue, Dec 20, 2005 at 04:46:16PM -0000, Mark Cave-Ayland wrote:
...
> Actually, looking through the source lwgeom_transform.c I see the following
> comment around line 500:
> 
> //this is *exactly* the same as PROJ.4's pj_transform(), but it doesn't
> //do the datum shift.
> 
> So does this mean that the PostGIS transform() function won't do datum
> shifts? Can anyone with a bit more PROJ4 knowledge say if this is correct,
> and if so, why PostGIS can't support it?

I'm not responsible for that code, anyway with a closer look, that
function is only invoked as a second attempt if pj_transform()
returns -38. Line 930:

                if (pj_errno == -38)  //2nd chance

Finding meaning of -38 would be the next thing to do ...

--strk;
Permalink | Reply |
headers
Frank Warmerdam | 20 Dec 22:56
Favicon

Re: transform() does not do work for NGO1948 zone III -> WGS84 longlat (epsg:27393->epsg:4326)

On 12/20/05, strk <at> refractions.net <strk <at> refractions.net> wrote:
> On Tue, Dec 20, 2005 at 04:46:16PM -0000, Mark Cave-Ayland wrote:
> ...
> > Actually, looking through the source lwgeom_transform.c I see the following
> > comment around line 500:
> >
> > //this is *exactly* the same as PROJ.4's pj_transform(), but it doesn't
> > //do the datum shift.
> >
> > So does this mean that the PostGIS transform() function won't do datum
> > shifts? Can anyone with a bit more PROJ4 knowledge say if this is correct,
> > and if so, why PostGIS can't support it?
>
> I'm not responsible for that code, anyway with a closer look, that
> function is only invoked as a second attempt if pj_transform()
> returns -38. Line 930:
>
>                 if (pj_errno == -38)  //2nd chance
>
> Finding meaning of -38 would be the next thing to do ...

Strk,

Error -38 is:
 	"failed to load NAD27-83 correction file",  	/* -38 */

So basically, if the datum shift files are not found, then the
code fallback to trying a conversion without the datum shift.

Best regards,
(Continue reading)

Permalink | Reply |
headers
strk | 21 Dec 09:19
Favicon

Re: transform() does not do work for NGO1948 zone III -> WGS84 longlat (epsg:27393->epsg:4326)

On Tue, Dec 20, 2005 at 04:56:47PM -0500, Frank Warmerdam wrote:
...
> Error -38 is:
>  	"failed to load NAD27-83 correction file",  	/* -38 */
> 
> So basically, if the datum shift files are not found, then the
> code fallback to trying a conversion without the datum shift.

Thank you Frank,
maybe we should raise a WARNING in that case ? I mean, is
"shift files not found" condition a sign of incomplete
installation of PROJ4 ?

--strk;

 /"\    ASCII Ribbon Campaign
 \ /    Respect for low technology.
  X     Keep e-mail messages readable by any computer system.
 / \    Keep it ASCII.
Permalink | Reply |
headers
Mark Cave-Ayland | 21 Dec 11:10

RE: transform() does not do work for NGO1948 zone III-> WGS84 longlat (epsg:27393->epsg:4326)


> -----Original Message-----
> From: postgis-users-bounces <at> postgis.refractions.net [mailto:postgis-users-
> bounces <at> postgis.refractions.net] On Behalf Of strk <at> refractions.net
> Sent: 21 December 2005 08:19
> To: PostGIS Users Discussion
> Subject: Re: [postgis-users] transform() does not do work for NGO1948 zone
> III-> WGS84 longlat (epsg:27393->epsg:4326)
> 
> On Tue, Dec 20, 2005 at 04:56:47PM -0500, Frank Warmerdam wrote:
> ...
> > Error -38 is:
> >  	"failed to load NAD27-83 correction file",  	/* -38 */
> >
> > So basically, if the datum shift files are not found, then the
> > code fallback to trying a conversion without the datum shift.
> 
> Thank you Frank,
> maybe we should raise a WARNING in that case ? I mean, is
> "shift files not found" condition a sign of incomplete
> installation of PROJ4 ?

I would say a warning should definitely be given if the transformation being
done is not the one that was requested by the user.

Kind regards,

Mark.

------------------------
(Continue reading)

Permalink | Reply |
headers
strk | 21 Dec 12:57
Favicon

Re: transform() does not do work for NGO1948 zone III-> WGS84 longlat (epsg:27393->epsg:4326)

On Wed, Dec 21, 2005 at 10:10:01AM -0000, Mark Cave-Ayland wrote:
> 
> > -----Original Message-----
> > From: postgis-users-bounces <at> postgis.refractions.net [mailto:postgis-users-
> > bounces <at> postgis.refractions.net] On Behalf Of strk <at> refractions.net
> > Sent: 21 December 2005 08:19
> > To: PostGIS Users Discussion
> > Subject: Re: [postgis-users] transform() does not do work for NGO1948 zone
> > III-> WGS84 longlat (epsg:27393->epsg:4326)
> > 
> > On Tue, Dec 20, 2005 at 04:56:47PM -0500, Frank Warmerdam wrote:
> > ...
> > > Error -38 is:
> > >  	"failed to load NAD27-83 correction file",  	/* -38 */
> > >
> > > So basically, if the datum shift files are not found, then the
> > > code fallback to trying a conversion without the datum shift.
> > 
> > Thank you Frank,
> > maybe we should raise a WARNING in that case ? I mean, is
> > "shift files not found" condition a sign of incomplete
> > installation of PROJ4 ?
> 
> I would say a warning should definitely be given if the transformation being
> done is not the one that was requested by the user.

I've added the warning:

	elog(WARNING, "transform: failed to load NAD27-83 correction file");
(Continue reading)

Permalink | Reply |
headers
Havard Tveite | 22 Dec 13:07

Re: transform() does not do work for NGO1948 zone III-> WGS84 longlat (epsg:27393->epsg:4326)

It is nice to provide warnings, but for epsg:27393 -> epsg:4326
the following problems persist:

1) Longitude of origin (+lon_0=-10.72291666666667) is ignored
    (as demonstrated by strk's query).
2) No datum shift seems to be performed.

Since the transformation runs smoothly using UMN Mapserver (
with dynamically linked proj - so it should be 4.4.9), there
is probably a problem somewhere.

I have an old PostGIS running (0.8 over 7.4.1 - proj statically
linked "old" version), and for that setup, the transformation
seems to be OK (including datum shift).

select AsText(
transform('SRID=27393;POINT(2836.394208896 185526.991186176)',4326)
);
 >  POINT(10.7732463298862 59.665688397594)

Håvard

strk <at> refractions.net wrote:
> On Wed, Dec 21, 2005 at 10:10:01AM -0000, Mark Cave-Ayland wrote:
> 
>>>-----Original Message-----
>>>From: postgis-users-bounces <at> postgis.refractions.net [mailto:postgis-users-
>>>bounces <at> postgis.refractions.net] On Behalf Of strk <at> refractions.net
>>>Sent: 21 December 2005 08:19
>>>To: PostGIS Users Discussion
(Continue reading)

Permalink | Reply |
headers
strk | 22 Dec 13:22
Favicon

Re: transform() does not do work for NGO1948 zone III-> WGS84 longlat (epsg:27393->epsg:4326)

mmm... I checked using transform_geometry() to exclude
the new proj4 caching code... It still gives the same
result:

strk=# select astext(transform('SRID=27393;POINT(2836.394208896 185526.991186176)',4326));
                   astext
--------------------------------------------
 POINT(0.0464201635443672 59.6661916850467)
(1 row)

strk=# select astext(transform_geometry('SRID=27393;POINT(2836.394208896 185526.991186176)',
get_proj4_from_srid(27393), get_proj4_from_srid(4326), 4326));
                   astext
--------------------------------------------
 POINT(0.0464201635443672 59.6661916850467)
(1 row)

This is with  Rel. 4.4.8, 3 May 2004, PostGIS 1.1.0.

PROJ4 text for 27393 is as follows:

 +proj=tmerc +lat_0=58 +lon_0=-10.72291666666667 +k=1.000000 +x_0=0 +y_0=0 +a=6377492.018
+b=6356173.508712696 +towgs84=278.3,93,474.5,7.889,0.05,-6.61,6.21 +pm=oslo +units=m +no_defs

PROJ4 text for 4326 is as follows:

 +proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs

This is both in the db and in the spatial_ref_sys.sql script in 1.1.0.
(Continue reading)

Permalink | Reply |
headers
Frank Warmerdam | 22 Dec 15:46
Favicon

Re: transform() does not do work for NGO1948 zone III-> WGS84 longlat (epsg:27393->epsg:4326)

On 12/22/05, Havard Tveite <havard.tveite <at> umb.no> wrote:
> It is nice to provide warnings, but for epsg:27393 -> epsg:4326
> the following problems persist:
>
> 1) Longitude of origin (+lon_0=-10.72291666666667) is ignored
>     (as demonstrated by strk's query).

Havard,

Are you sure that the the +lon_0 is being ignored?  It seems more
likely that it is being undone by the +pm=oslo.  There are a variety of
problems with prime meridian handling in the EPSG->proj.4 translation
process.  I would suggest modifying the PROJ.4 definition to not
include the +pm=oslo and see if that helps.

> 2) No datum shift seems to be performed.
>
> Since the transformation runs smoothly using UMN Mapserver (
> with dynamically linked proj - so it should be 4.4.9), there
> is probably a problem somewhere.

You are likely right, I'm not sure what the story is.

Best regards,
--
---------------------------------------+--------------------------------------
I set the clouds in motion - turn up   | Frank Warmerdam, warmerdam <at> pobox.com
light and sound - activate the windows | http://pobox.com/~warmerdam
and watch the world go round - Rush    | Geospatial Programmer for Rent
(Continue reading)

Permalink | Reply |
headers
Michael Fuhr | 22 Dec 17:54
Favicon

Re: transform() does not do work for NGO1948 zone III-> WGS84 longlat (epsg:27393->epsg:4326)

On Thu, Dec 22, 2005 at 09:46:48AM -0500, Frank Warmerdam wrote:
> On 12/22/05, Havard Tveite <havard.tveite <at> umb.no> wrote:
> > It is nice to provide warnings, but for epsg:27393 -> epsg:4326
> > the following problems persist:
> >
> > 1) Longitude of origin (+lon_0=-10.72291666666667) is ignored
> >     (as demonstrated by strk's query).
> 
> Are you sure that the the +lon_0 is being ignored?  It seems more
> likely that it is being undone by the +pm=oslo.  There are a variety of
> problems with prime meridian handling in the EPSG->proj.4 translation
> process.  I would suggest modifying the PROJ.4 definition to not
> include the +pm=oslo and see if that helps.

I don't think anybody has mentioned the PROJ.4 program cs2cs, which
you can use to do conversions from the command line; this is a
convenient way to remove PostGIS from the conversion process to see
if it's the problem.

Here's an example that converts from 27393 to 4326:

echo 2836.394208896 185526.991186176 | 
cs2cs +proj=tmerc +lat_0=58 +lon_0=-10.72291666666667 \
      +k=1.000000 +x_0=0 +y_0=0 +a=6377492.018 +b=6356173.508712696 \
      +towgs84=278.3,93,474.5,7.889,0.05,-6.61,6.21 \
      +pm=oslo \
      +units=m +no_defs \
      +to \
      +proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs \
      -f "%.8f"
(Continue reading)

Permalink | Reply |
headers
Havard Tveite | 23 Dec 15:19

Re: transform() does not do work for NGO1948 zone III-> WGS84 longlat (epsg:27393->epsg:4326)

An important omission in my earlier emails (sorry...):
The proj.4 epsg +lon_0 parameters (longitude offset) for the NGO
projections have been strange for the last couple of proj.4
releases. This has been reported as a bug in proj.4 bugzilla.
I have fixed the proj.4 epgs parameter file by hand to overcome
the problems.

I guess that the PostGIS spatial_ref_sys table is not
generated from the local proj.4 epsg file...
The longitude offset probably stems from this.

The epsg parameter file changes seems to have been made some
time after proj-4.4.7 (I think this version had correct values
for both latitude offset and towgs84).
Where the epsg "errors" come from, I don't know. It could be
an error introduced when generating the proj.4 epsg file, or
it could be an error in the EPSG sources. The fact that "NGO 1948"
and "NGO 1948 (Oslo)" use different prime meridians (Greenwich
and Oslo respectively) might be one of the sources of confusion.

I can still not explain the datum shift problems (or what
ever problem it might be).

Here are my local values for PostGIS 0.8 and 1.1:

select * from spatial_ref_sys where srid = 27393;

741, postgis 0.8:
 > 27393 | EPSG      |     27393 |
PROJCS["NGO_1948_Oslo_Norway_Zone_3",GEOGCS["GCS_NGO_1948_Oslo",DATUM["D_NGO_1948",SPHEROID["Bessel_Modified",6377492.018,299.1528128]],PRIMEM["Oslo",10.72291666666667],UNIT["Degree",0.017453292519943295]],PROJECTION["Transverse_Mercator"],PARAMETER["False_Easting",0],PARAMETER["False_Northing",0],PARAMETER["Central_Meridian",0],PARAMETER["Scale_Factor",1],PARAMETER["Latitude_Of_Origin",58],UNIT["Meter",1]]
(Continue reading)

Permalink | Reply |
headers
Havard Tveite | 3 Jan 19:55

Re: transform() does not do work for NGO1948 zone III-> WGS84 longlat (epsg:27393->epsg:4326)

Just to close this...

I have tested some more, and it seems as if transform() in
PostGIS 1.0.x is working as expected. When I changed +lon_0
in proj4text of the NGO projections in the spatial_ref_sys
to the correct values, the result is consistent with proj.4
cs2cs.
I have also checked with an independent tranformation
service (Norwegian Mapping Authority), and it gives about
the same results.

As soon as the epsg-file bug of proj.4 has been fixed,
PostGIS users with NGO1948 data should be happy.

Sorry about the noice...

Havard Tveite wrote:
> An important omission in my earlier emails (sorry...):
> The proj.4 epsg +lon_0 parameters (longitude offset) for the NGO
> projections have been strange for the last couple of proj.4
> releases. This has been reported as a bug in proj.4 bugzilla.
> I have fixed the proj.4 epgs parameter file by hand to overcome
> the problems.
> 
> I guess that the PostGIS spatial_ref_sys table is not
> generated from the local proj.4 epsg file...
> The longitude offset probably stems from this.
> 
> The epsg parameter file changes seems to have been made some
> time after proj-4.4.7 (I think this version had correct values
(Continue reading)

Permalink | Reply |
headers
Frank Warmerdam | 21 Dec 16:31
Favicon

Re: transform() does not do work for NGO1948 zone III -> WGS84 longlat (epsg:27393->epsg:4326)

On 12/21/05, strk <at> refractions.net <strk <at> refractions.net> wrote:
> Thank you Frank,
> maybe we should raise a WARNING in that case ? I mean, is
> "shift files not found" condition a sign of incomplete
> installation of PROJ4 ?

Strk,

In PROJ.4 an error is raised because pj_transform() terminates.
In your case you might try translating this to a warning, since you
try to workaround the issue by doing the no-datum-shift conversion.

On 12/21/05, Patrick <pvanlaake <at> users.sourceforge.net> wrote:
> Why would the error be thrown in the first place, seeing that we are going
> from NGO1948/Bessel to WGS84/WGS84? Where does the grid shift enter?

I don't think error -38 is occuring in your case.  I believe Strk just
pointed out that this is the only situation in which the fallback logic
in the postgis transform() function is used.

Best regards,
--
---------------------------------------+--------------------------------------
I set the clouds in motion - turn up   | Frank Warmerdam, warmerdam <at> pobox.com
light and sound - activate the windows | http://pobox.com/~warmerdam
and watch the world go round - Rush    | Geospatial Programmer for Rent
Permalink | Reply |
headers
Patrick | 21 Dec 10:23

Re: transform() does not do work for NGO1948 zone III-> WGS84 longlat (epsg:27393->epsg:4326)


"Frank Warmerdam" <warmerdam <at> pobox.com> wrote in message 
news:931f8ea90512201356v75c57d08ub18b6827460ee16e <at> mail.gmail.com...
>
> Error -38 is:
>   "failed to load NAD27-83 correction file",  /* -38 */
>
> So basically, if the datum shift files are not found, then the
> code fallback to trying a conversion without the datum shift.

Frank,

Why would the error be thrown in the first place, seeing that we are going 
from NGO1948/Bessel to WGS84/WGS84? Where does the grid shift enter?

Patrick
Permalink | Reply |

Gmane