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:
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 |
> 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 |
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 |
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 |
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 |
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 |
> 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 |
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. 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 |
Free forum by Nabble | Edit this page |