PostGIS users

Havard Tveite | 20 Dec 16:11

Picon

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

All the relevant information seems to be present in the spatial_ref_sys table:

select * from spatial_ref_sys where srid=27393;
 > 27393 | EPSG      |     27393 | PROJCS["NGO 1948 (Oslo) / NGO zone III",GEOGCS["NGO 1948
(Oslo)",DATUM["NGO_1948_Oslo",SPHEROID["Bessel
Modified",6377492.018,299.1528128,AUTHORITY["EPSG","7005"]],TOWGS84[278.3,93,474.5,7.889,0.05,-6.61,6.21],AUTHORITY["EPSG","6817"]],PRIMEM["Oslo",10.72291666666667,AUTHORITY["EPSG","8913"]],UNIT["degree",0.01745329251994328,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4817"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",58],PARAMETER["central_meridian",0],PARAMETER["scale_factor",1],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AUTHORITY["EPSG","27393"]]
| +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

Is there a problem with transform(), or is it the spatial_ref_sys
entry that is wrong, or am I doing something wrong here?

I have postgresql CVS from some days ago and postgis CVS from today
(solaris 2.7, gcc 3.3.2):
select postgis_full_version();
POSTGIS="1.1.0" GEOS="2.1.4" PROJ="Rel. 4.4.9, 29 Oct 2004" USE_STATS

--

-- 
Håvard Tveite
Department of Mathematical Sciences and Technology, UMB
Drøbakveien 14, POBox 5003, N-1432 Ås, NORWAY
Phone: +47 64965483 Fax: +47 64965401 http://www.umb.no/imt

Mark Cave-Ayland | 20 Dec 17:46

Picon

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,

Mark.

------------------------
WebBased Ltd
17 Research Way
Plymouth
PL6 8BT

T: +44 (0)1752 797131
F: +44 (0)1752 791023

http://www.webbased.co.uk   
http://www.infomapper.com
http://www.swtc.co.uk  

This email and any attachments are confidential to the intended recipient
and may also be privileged. If you are not the intended recipient please
delete it from your system and notify the sender. You should not copy it or
use it for any purpose nor disclose or distribute its contents to any other
person.

strk | 20 Dec 22:19

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;

Frank Warmerdam | 20 Dec 22:56

Picon

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,
--
---------------------------------------+--------------------------------------
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

Patrick | 21 Dec 10:23

Picon

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

strk | 21 Dec 09:19

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.

Frank Warmerdam | 21 Dec 16:31

Picon

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

Mark Cave-Ayland | 21 Dec 11:10

Picon

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.

------------------------
WebBased Ltd
17 Research Way
Plymouth
PL6 8BT

T: +44 (0)1752 797131
F: +44 (0)1752 791023

http://www.webbased.co.uk   
http://www.infomapper.com
http://www.swtc.co.uk  

This email and any attachments are confidential to the intended recipient
and may also be privileged. If you are not the intended recipient please
delete it from your system and notify the sender. You should not copy it or
use it for any purpose nor disclose or distribute its contents to any other
person.

strk | 21 Dec 12:57

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");

But, I cannot reproduce the problem with this query:

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

Even if I temporarly remove /usr/local/share/proj  ...
Output is:

	POINT(0.0464201635443672 59.6661916850467)

postgis_proj_version(): Rel. 4.4.8, 3 May 2004

Tips on how to reproduce the error (or is it possible PROJ4 would
make the two attempts by itself ?) ?

--strk;

Havard Tveite | 22 Dec 13:07

Picon

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
>>>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");
> 
> But, I cannot reproduce the problem with this query:
> 
> 	select AsText(
> 		transform('SRID=27393;POINT(2836.394208896 185526.991186176)',
> 			4326)
> 	);
> 
> Even if I temporarly remove /usr/local/share/proj  ...
> Output is:
> 
> 	POINT(0.0464201635443672 59.6661916850467)
> 
> postgis_proj_version(): Rel. 4.4.8, 3 May 2004
> 
> Tips on how to reproduce the error (or is it possible PROJ4 would
> make the two attempts by itself ?) ?
> 
> --strk;
> _______________________________________________
> postgis-users mailing list
> postgis-users <at> postgis.refractions.net
> http://postgis.refractions.net/mailman/listinfo/postgis-users
> 

--

-- 
Håvard Tveite
Department of Mathematical Sciences and Technology, UMB
Drøbakveien 14, POBox 5003, N-1432 Ås, NORWAY
Phone: +47 64965483 Fax: +47 64965401 http://www.umb.no/imt

Frank Warmerdam | 22 Dec 15:46

Picon

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

Havard Tveite | 23 Dec 15:19

Picon

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]]
| +proj=tmerc +lat_0=58 +lon_0=10.72291666666667 +k=1.000000 +x_0=0 +y_0=0 +a=6377492.018
+b=6356173.508712696 +pm=10.72291666666667 +units=m

82 cvs, postgis cvs:
 > 27393 | EPSG      |     27393 | PROJCS["NGO 1948 (Oslo) / NGO zone III",GEOGCS["NGO 1948
(Oslo)",DATUM["NGO_1948_Oslo",SPHEROID["Bessel
Modified",6377492.018,299.1528128,AUTHORITY["EPSG","7005"]],TOWGS84[278.3,93,474.5,7.889,0.05,-6.61,6.21],AUTHORITY["EPSG","6817"]],PRIMEM["Oslo",10.72291666666667,AUTHORITY["EPSG","8913"]],UNIT["degree",0.01745329251994328,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4817"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",58],PARAMETER["central_meridian",0],PARAMETER["scale_factor",1],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AUTHORITY["EPSG","27393"]]
| +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

My share/proj/epsg (assuming that +lon_0 is supposed to give
the offset from +pm):
# NGO 1948 (Oslo) / NGO zone III
<27393> +proj=tmerc +lat_0=58 +lon_0=0 +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  <>

Michael Fuhr suggested to run cs2cs. With the correct lon_0
(0 for NGO 1984 zone 3), the results are:
echo 2836.394208896 185526.991186176 |
cs2cs +proj=tmerc +lat_0=58 +lon_0=0 \
       +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"
 > 10.76845023     59.66650218 0.30776738

Which still seems to be quite far (a few hundred meters) off.

Very strange that my UMN Mapserver (4.6.1) seems to perform a much
better transformation (within a couple of meteres)... (I have
checked visually how the NGO 1948 zone 3 dataset lines up with
datasets with other coordinate systems when transformed to WGS84
longlat).

Perhaps I will have to do some more homework before pursuing this
further...

Thank you for your efforts!

Haavard

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).
> 
> 
> 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
> _______________________________________________
> postgis-users mailing list
> postgis-users <at> postgis.refractions.net
> http://postgis.refractions.net/mailman/listinfo/postgis-users
> 

--

-- 
Håvard Tveite
Department of Mathematical Sciences and Technology, UMB
Drøbakveien 14, POBox 5003, N-1432 Ås, NORWAY
Phone: +47 64965483 Fax: +47 64965401 http://www.umb.no/imt

Havard Tveite | 3 Jan 19:55

Picon

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
> 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]] 
> | +proj=tmerc +lat_0=58 +lon_0=10.72291666666667 +k=1.000000 +x_0=0 
> +y_0=0 +a=6377492.018 +b=6356173.508712696 +pm=10.72291666666667 +units=m
> 
> 82 cvs, postgis cvs:
>  > 27393 | EPSG      |     27393 | PROJCS["NGO 1948 (Oslo) / NGO zone 
> III",GEOGCS["NGO 1948 (Oslo)",DATUM["NGO_1948_Oslo",SPHEROID["Bessel 
>
Modified",6377492.018,299.1528128,AUTHORITY["EPSG","7005"]],TOWGS84[278.3,93,474.5,7.889,0.05,-6.61,6.21],AUTHORITY["EPSG","6817"]],PRIMEM["Oslo",10.72291666666667,AUTHORITY["EPSG","8913"]],UNIT["degree",0.01745329251994328,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4817"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",58],PARAMETER["central_meridian",0],PARAMETER["scale_factor",1],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AUTHORITY["EPSG","27393"]] 
> | +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
> 
> My share/proj/epsg (assuming that +lon_0 is supposed to give
> the offset from +pm):
> # NGO 1948 (Oslo) / NGO zone III
> <27393> +proj=tmerc +lat_0=58 +lon_0=0 +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  <>
> 
> Michael Fuhr suggested to run cs2cs. With the correct lon_0
> (0 for NGO 1984 zone 3), the results are:
> echo 2836.394208896 185526.991186176 |
> cs2cs +proj=tmerc +lat_0=58 +lon_0=0 \
>       +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"
>  > 10.76845023     59.66650218 0.30776738
> 
> Which still seems to be quite far (a few hundred meters) off.
> 
> Very strange that my UMN Mapserver (4.6.1) seems to perform a much
> better transformation (within a couple of meteres)... (I have
> checked visually how the NGO 1948 zone 3 dataset lines up with
> datasets with other coordinate systems when transformed to WGS84
> longlat).
> 
> Perhaps I will have to do some more homework before pursuing this
> further...
> 
> Thank you for your efforts!
> 
> 
> Haavard
> 
> 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).
>>
>>
>>
>> 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
>> _______________________________________________
>> postgis-users mailing list
>> postgis-users <at> postgis.refractions.net
>> http://postgis.refractions.net/mailman/listinfo/postgis-users
>>
> 

--

-- 
Håvard Tveite
Department of Mathematical Sciences and Technology, UMB
Drøbakveien 14, POBox 5003, N-1432 Ås, NORWAY
Phone: +47 64965483 Fax: +47 64965401 http://www.umb.no/imt

Michael Fuhr | 22 Dec 17:54

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"

0.04642016      59.66619169 -6.08805941

These are the same numbers people have been seeing in PostGIS (I
think the third coordinate is Z).  Here's the same example,
omitting +pm=oslo:

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 \
      +units=m +no_defs \
      +to \
      +proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs \
      -f "%.8f"

-10.67553750    59.66595692 -17.39326135

Here's the example with +pm=oslo and without +lon_0:

echo 2836.394208896 185526.991186176 | 
cs2cs +proj=tmerc +lat_0=58 \
      +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"

10.76845023     59.66650218 0.30776738

These aren't exactly the numbers Håvard got in the old version
of PostGIS but they're getting close.

--

-- 
Michael Fuhr

strk | 22 Dec 13:22

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.

Can you check corrispondent values in your old *working* setup ?

--strk;

On Thu, Dec 22, 2005 at 01:07:17PM +0100, Havard Tveite 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).
> 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

Return

Return to gmane.comp.gis.postgis.

Advertisement

Project Web Page

Search Archive

Language

Change language

Options

Current view: All messages / Showing whole messages / Hiding cited text.
Change to Threads only, shortened messages, or show cited text.

Post a message
NNTP Newsgroup
Classic Gmane web interface
XML

XML

RSS Feed
List Information

About Gmane