Custom projection wildly inaccurate

classic Classic list List threaded Threaded
8 messages Options
Reply | Threaded
Open this post in threaded view
|

Custom projection wildly inaccurate

Thy, Kristian
Hi list,
 
(Apologies for sending HTML email - I'm stuck on OWA ...)
 
I'm trying to reproject data in an obsolete Danish orthogonal projection called System34 in MapServer. To this end, I've found a .prj file for ArcMap distributed by the Danish National Survey & Cadastre, which defines the following parameters:
 
PROJCS["Approximation of non-projective System 34 Sjælland",
   GEOGCS["GCS_WGS_1984",
      DATUM["D_WGS_1984", SPHEROID["WGS_1984",6378137.0,298.257223563]],
      PRIMEM["Greenwich",0.0],
      UNIT["Degree",0.0174532925199433]
   ],
   PROJECTION["Orthographic"],
   PARAMETER["False_Easting",-100000.0],
   PARAMETER["False_Northing",100000.0],
   PARAMETER["Longitude_Of_Center",12.105661],
   PARAMETER["Latitude_Of_Center",55.29235],
   UNIT["Meter",1.0]
]
 
To get MapServer to reproject my data, I added the following code to ms4w\proj\nad\epsg:
 
<34005> +proj=ortho +lat_0=55.29235 +lon_0=12.105661 +x_0=-100000.0 +y_0=100000.0
 
Problem is that this transposes the data many kilometres north and slightly west (and perhaps also distorts it slightly, I'm not sure). The effect can be seen in these three maps:
 
http://quovadis.dk/img/proj/
 
All three maps are drawn in EPSG:32632. The green line data is natively EPSG:32632. The blue points are fetched from a WMS service in EPSG:4326, EPSG:32632 and my custom transformation respectively. The difference between lat/lon and UTM is noticeable but acceptable, but my own projection is way off. I have confirmed the distortion with different data sources in System 34 - this WMS and an orthophoto from another source - so I'm certain the error is on my end. Can anybody tell me which rookie mistake I've commited? ;-)
 
 
 
best regards,
Kristian Thy


This email and any attached files are confidential and copyright protected. If you are not the addressee, any dissemination of this communication is strictly prohibited. Unless otherwise expressly agreed in writing, nothing stated in this communication shall be legally binding.


_______________________________________________
Proj mailing list
[hidden email]
http://lists.maptools.org/mailman/listinfo/proj
Reply | Threaded
Open this post in threaded view
|

Re: Custom projection wildly inaccurate

OvV_HN
> From: "Thy, Kristian"
> Date: Sat, 18 Nov 2006 01:17:16 +0100
> Subject: [Proj] Custom projection wildly inaccurate

 > I'm trying to reproject data in an obsolete Danish orthogonal projection
> called System34 in MapServer. To this end, I've found a .prj file for ArcMap
> distributed by the Danish National Survey & Cadastre, which defines the
> following parameters:
....
 
> Problem is that this transposes the data many kilometres north and slightly
> west (and perhaps also distorts it slightly, I'm not sure). The effect can be
> seen in these three maps:
...

Another problem is that System 34 is not coupled to any projection. I have
heard of a transverse mercator approximation, but an orthographic?
The Danish use a 11th (13th?) degree polynomial to convert System 34
coordinates to UTM. I haven't seen this polynomial roaming on the internet,
so I cannot comment on the accuracy of approximations.
There are 3 zones in System 34 and the coordinate system is west-oriented.

Kort & Matrikelstyrelsen - Mere om System 34
<http://www.kms.dk/Referencenetogopmaaling/Referencesystemer/System_34/mere_
om_system34.htm>
On-line transformations:
<http://valdemar.kms.dk/trf/>

_______________________________________________
Proj mailing list
[hidden email]
http://lists.maptools.org/mailman/listinfo/proj
Reply | Threaded
Open this post in threaded view
|

RE: Custom projection wildly inaccurate

Thy, Kristian
From: Oscar van Vlijmen
> Another problem is that System 34 is not coupled to any
> projection. I have heard of a transverse mercator approximation,
> but an orthographic?

Now that you mention it, this seems extremely strange. But I just went
with what was in the official .prj files. You wouldn't happen to have
any information about the transverse mercator approximation?

> The Danish use a 11th (13th?) degree polynomial to convert System 34
> coordinates to UTM. I haven't seen this polynomial roaming on
> the internet, so I cannot comment on the accuracy of approximations.

I had them in my notes from a geodesy course in uni, but they have
unfortunately been mislaid. At any rate, I am not sure of the copyright
implications of publishing these polynomials from the tech note I got
them in.

However, using the coordinate transformer thingy on valdemar.kms.dk,
I get the following results for the top left corner of a TIFF orthophoto
I have:

System 34 Sjælland (Zealand): (-86714.42, 162554.98)
Transformed to lon/lat (WGS84): (12d20'42.97906"E, 55d52'40.50251"N)

(The negative X is due to the system being west-oriented, but for ease
of use a "normal" east-oriented X-axis is used with all X coords then
being negative. These coordinates are confirmed using KMS's KMSTrans
program, a desktop application for conversion of coordinates.)

Using cs2cs with the orthographic projection inferred from the .prj file:

C:\>cs2cs -E +proj=ortho +lat_0=55.29235 +lon_0=12.105661 +x_0=-100000.0 +y_0=100000.0 +to +proj=latlong +datum=WGS84 cs.txt
-86714.420 162554.980   12d19'5.829"E   56d1'56.147"N 14670.382

Clearly, that projection doesn't fit the bill. But when I pull the same
TIFF into ArcMap and examine the spatial reference for it, this is what
I see:

http://quovadis.dk/img/proj/srs.jpg

The orthophoto at bottom right is the original in System 34, and also the
one for which the SRS is shown in the dialog. At top right you see the
same orthophoto reprojected using GDAL with what I thought were the
parameters shown in the .prj file.

Any help in unravelling this mystery is much appreciated. I suspect that
the problem lies in my interpretation of the ArcMap .prj file.

best regards,
Kristian Thy


This email and any attached files are confidential and copyright protected. If you are not the addressee, any dissemination of this communication is strictly prohibited. Unless otherwise expressly agreed in writing, nothing stated in this communication shall be legally binding.

_______________________________________________
Proj mailing list
[hidden email]
http://lists.maptools.org/mailman/listinfo/proj
Reply | Threaded
Open this post in threaded view
|

RE: Custom projection wildly inaccurate

Thy, Kristian
From: Thy, Kristian
> You wouldn't happen to have any information about the transverse mercator approximation?

I think I just stumbled upon something useful:

http://www.gfy.ku.dk/~cct/4geo_k24.htm - at the bottom:

"Det Danske koordinatsystem System 34 er ikke baseret på nogen
 kortprojektion i matematisk forstand. Basis er en tilnærmet
 transversal konform cylinderprojektion (samt 1 for Bornholm).
 Jorden er regnet kugleformet, baseret på middelkrumningsradius
 for Hayford-ellipsoiden (a = 6378388 m, 1/f = 297) i bredden
 56o 08' for Jylland og 55o 20' for Sjælland. Det trigonometriske
 punkt Agri Baunehøj har fået koordinaterne (Y, X) = ( 200 km,
 200 km), så alle danske koordinater er positive. X - aksen peger
 mod Nord og Y - aksen peger mod vest. Afbildningen er iøvrigt
 fastlagt så retningen fra Agri til Lysnet er 24o31'14".17. "

Transl.:

"The Danish coordinate system System 34 is not based on a
 map projection in the mathematical sense. It's basis is an
 approximate transverse conformal cylindrical projection
 (plus another for Bornholm). Earth is defined as spherical,
 based on the average curvature radius for the Hayford
 ellipsoid (a = 6378388 m, 1/f = 297) at latitude 56d08' for
 Jutland and 55d20' for Zealand. The trigonometric point Agri
 Baunehøj is given coordinates (Y, X) = (200km, 200km) such
 that all Danish coordinates are positive. The X axis points
 North and the Y axis points west. The projection is defined
 such that the bearing from Agri to Lysnet is 24d31'14.17". "

For practical purposes, what is officially the X and Y axis in
this projection is switched to conform with other coordinate
systems, and the X axis is then counted as negative from its
origin East of Denmark, meaning that Agri has coordinates
(-200000, 200000) in my data.

The last sentence implies that the Y axis in this projection is
not parallel to geographical North (great). I can't find a value
for this deviance, but I'm willing to forgo that for now, if
someone can help me make up parameters for this projection.

I've tried

C:\>cs2cs -E +proj=tmerc +lat_0=55.29235 +lon_0=12.105661 +x_0=-100000.0 +y_0=100000.0 +to +proj=latlong +datum=WGS84 cs.txt

which gives me

-86714.420 162554.980   12d19'4.069"E   55d51'14.525"N 0.000

but I'm still a bit away from the correct result of
12d20'42.97906" 55d52'40.50251" - I can see I'm missing the
scale factor from the parameters (+k), but I can't find a
stated value for this anywhere.

Kristian


This email and any attached files are confidential and copyright protected. If you are not the addressee, any dissemination of this communication is strictly prohibited. Unless otherwise expressly agreed in writing, nothing stated in this communication shall be legally binding.

_______________________________________________
Proj mailing list
[hidden email]
http://lists.maptools.org/mailman/listinfo/proj
Reply | Threaded
Open this post in threaded view
|

Re: Custom projection wildly inaccurate

OvV_HN
In reply to this post by Thy, Kristian
> From: "Thy, Kristian"
> Date: Mon, 20 Nov 2006 10:13:00 +0100
> Subject: RE: [Proj] Custom projection wildly inaccurate
 
> From: Oscar van Vlijmen
>> Another problem is that System 34 is not coupled to any
>> projection. I have heard of a transverse mercator approximation,
>> but an orthographic?
> Now that you mention it, this seems extremely strange. But I just went
> with what was in the official .prj files. You wouldn't happen to have
> any information about the transverse mercator approximation?

There seem to be several approximations. Most straight (read: what you can
do per proj/cs2cs command line alone) projective approximations are bad.
There is really only one way to do it correctly: use the official polynomial
solution.

An interesting approximation, using transverse mercator with some pre- &
postprocessing, can be read in:
Kp2000 - en nødvendighed eller ikke?
Dansk Vejtidsskrift, (6/7) 2001, p. 18-19
Leif Kahl Kristensen, Institut for fysik og astronomi
<http://asp.vejtid.dk/Artikler/2001/06-07%5C2930.pdf>

He gives an approximation for two System34 zones: Sjælland and Jylland/Fynn.
Using his approximation I got for the coordinates you mentioned:
lon =12d 20m 42.97906s; lat = 55d 52m 40.50251s; wgs84 datum.
-> x = 86714.465; y = 162554.937 m
The values you mentioned were:
-86714.42, 162554.98
Kristensen states that his procedure results in an average difference in X
of 5 cm, in Y of 9 cm, with a maximum difference of 18 cm.
For his Jylland/Fynn approximation the differences are larger.

>> The Danish use a 11th (13th?) degree polynomial to convert System 34
>> coordinates to UTM. I haven't seen this polynomial roaming on
>> the internet, so I cannot comment on the accuracy of approximations.
> I had them in my notes from a geodesy course in uni, but they have
> unfortunately been mislaid. At any rate, I am not sure of the copyright
> implications of publishing these polynomials from the tech note I got
> them in.
Algorithms and formulae are usually not copyrightable, but they are
patentable. So, still publishable, more or less, but if patented, commercial
exploitation and sometimes other forms of exploitation are prohibited, at
least in the countries where the patent is valid.
Note: I'm no lawyer, so this is no legal advise!


> I think I just stumbled upon something useful:
> http://www.gfy.ku.dk/~cct/4geo_k24.htm - at the bottom:
I've read this too, but I didn't give it a try.





_______________________________________________
Proj mailing list
[hidden email]
http://lists.maptools.org/mailman/listinfo/proj
Reply | Threaded
Open this post in threaded view
|

RE: Custom projection wildly inaccurate

Thy, Kristian
From: Oscar van Vlijmen
> There seem to be several approximations. Most straight (read:
> what you can do per proj/cs2cs command line alone) projective
> approximations are bad.
> There is really only one way to do it correctly: use the
> official polynomial solution.

Yes, but that's not easily translatable for use in Proj-based
applications like GDAL and MapServer (as far I know, at least).
We will most likely have to convert the data manually, but in
the meantime it would be nice to have a bad-and-fast projection
for demonstration and testing purposes.

> An interesting approximation, using transverse mercator with
> some pre- & postprocessing, can be read in:
> Kp2000 - en nødvendighed eller ikke?
> Dansk Vejtidsskrift, (6/7) 2001, p. 18-19
> Leif Kahl Kristensen, Institut for fysik og astronomi
> <http://asp.vejtid.dk/Artikler/2001/06-07%5C2930.pdf>
>
> He gives an approximation for two System34 zones: Sjælland
> and Jylland/Fynn.
> Using his approximation I got for the coordinates you mentioned:
> lon =12d 20m 42.97906s; lat = 55d 52m 40.50251s; wgs84 datum.
> -> x = 86714.465; y = 162554.937 m
> The values you mentioned were:
> -86714.42, 162554.98
> Kristensen states that his procedure results in an average
> difference in X of 5 cm, in Y of 9 cm, with a maximum difference
> of 18 cm.
> For his Jylland/Fynn approximation the differences are larger.

Interesting. For our application, that error would be acceptable.
My next question would then be if I can then somehow use the
formulae presented to warp an orthophoto using gdalwarp, but I
guess that's a question for another list :-)

> Algorithms and formulae are usually not copyrightable, but they are
> patentable. So, still publishable, more or less, but if
> patented, commercial exploitation and sometimes other forms
> of exploitation are prohibited, at least in the countries
> where the patent is valid.
> Note: I'm no lawyer, so this is no legal advise!

Makes sense. But them again, IANAL (either) :-)

best regards,
Kristian


This email and any attached files are confidential and copyright protected. If you are not the addressee, any dissemination of this communication is strictly prohibited. Unless otherwise expressly agreed in writing, nothing stated in this communication shall be legally binding.

_______________________________________________
Proj mailing list
[hidden email]
http://lists.maptools.org/mailman/listinfo/proj
Reply | Threaded
Open this post in threaded view
|

Re: Custom projection wildly inaccurate

Greg Troxel

  > There seem to be several approximations. Most straight (read:
  > what you can do per proj/cs2cs command line alone) projective
  > approximations are bad.
  > There is really only one way to do it correctly: use the
  > official polynomial solution.

  Yes, but that's not easily translatable for use in Proj-based
  applications like GDAL and MapServer (as far I know, at least).
  We will most likely have to convert the data manually, but in
  the meantime it would be nice to have a bad-and-fast projection
  for demonstration and testing purposes.

One thing that would work, but perhaps set a bad precedent, is to
define a new projection in the proj source, that knows the
polynomials.  This would be a particular transform, rather than a
generic projection with parameters.  Alternatively, there could be a
polynomial transform with vast numbers of parameters.


--
    Greg Troxel <[hidden email]>

_______________________________________________
Proj mailing list
[hidden email]
http://lists.maptools.org/mailman/listinfo/proj

attachment0 (191 bytes) Download Attachment
Reply | Threaded
Open this post in threaded view
|

Re: Custom projection wildly inaccurate

Martin Vermeer
On Mon, 2006-11-20 at 08:36 -0500, Greg Troxel wrote:

>   > There seem to be several approximations. Most straight (read:
>   > what you can do per proj/cs2cs command line alone) projective
>   > approximations are bad.
>   > There is really only one way to do it correctly: use the
>   > official polynomial solution.
>
>   Yes, but that's not easily translatable for use in Proj-based
>   applications like GDAL and MapServer (as far I know, at least).
>   We will most likely have to convert the data manually, but in
>   the meantime it would be nice to have a bad-and-fast projection
>   for demonstration and testing purposes.
>
> One thing that would work, but perhaps set a bad precedent, is to
> define a new projection in the proj source, that knows the
> polynomials.  This would be a particular transform, rather than a
> generic projection with parameters.  Alternatively, there could be a
> polynomial transform with vast numbers of parameters.
Alternatively, generate a shift grid containing the shifts produced by
the polynomial. This is the approach I am considering for the Finnish
kkj, which is officially represented by a multi-triangle (Delaunay)
affine transformation.

- Martin


_______________________________________________
Proj mailing list
[hidden email]
http://lists.maptools.org/mailman/listinfo/proj

signature.asc (196 bytes) Download Attachment