Understand shapefiles in different projections (notably OSGB)
[potlatch2.git] / com / gradoservice / proj4as / proj / ProjSoMerc.as
1 /*******************************************************************************\r
2 NAME                       SWISS OBLIQUE MERCATOR\r
3 \r
4 PURPOSE:        Swiss projection.\r
5 WARNING:  X and Y are inverted (weird) in the swiss coordinate system. Not\r
6    here, since we want X to be horizontal and Y vertical.\r
7 \r
8 ALGORITHM REFERENCES\r
9 1. "Formules et constantes pour le Calcul pour la\r
10  projection cylindrique conforme Г  axe oblique et pour la transformation entre\r
11  des systГЁmes de rГ©fГ©rence".\r
12  http://www.swisstopo.admin.ch/internet/swisstopo/fr/home/topics/survey/sys/refsys/switzerland.parsysrelated1.31216.downloadList.77004.DownloadFile.tmp/swissprojectionfr.pdf\r
13 \r
14 *******************************************************************************/\r
15 \r
16 \r
17 package com.gradoservice.proj4as.proj\r
18 {\r
19         import com.gradoservice.proj4as.ProjPoint;\r
20         import com.gradoservice.proj4as.ProjConstants;\r
21         import com.gradoservice.proj4as.Datum;\r
22                 \r
23         public class ProjSoMerc extends AbstractProjProjection\r
24         {\r
25                 private var lambda0:Number;\r
26                 private var R:Number;\r
27                 private var b0:Number;\r
28                 private var K:Number;\r
29                 \r
30                 \r
31                 public function ProjSoMerc(data:ProjParams)\r
32                 {\r
33                         super(data);\r
34                 }\r
35                 \r
36                 \r
37                   override public function init():void\r
38                   {\r
39                     var phy0:Number = this.lat0;\r
40                     this.lambda0 = this.long0;\r
41                     var sinPhy0:Number = Math.sin(phy0);\r
42                     var semiMajorAxis:Number = this.a;\r
43                     var invF:Number = this.rf;\r
44                     var flattening:Number = 1 / invF;\r
45                     var e2:Number = 2 * flattening - Math.pow(flattening, 2);\r
46                     var e:Number = this.e = Math.sqrt(e2);\r
47                     this.R = semiMajorAxis * Math.sqrt(1 - e2) / (1 - e2 * Math.pow(sinPhy0, 2.0));\r
48                     this.alpha = Math.sqrt(1 + e2 / (1 - e2) * Math.pow(Math.cos(phy0), 4.0));\r
49                     this.b0 = Math.asin(sinPhy0 / this.alpha);\r
50                     this.K = Math.log(Math.tan(Math.PI / 4.0 + this.b0 / 2.0))\r
51                             - this.alpha\r
52                             * Math.log(Math.tan(Math.PI / 4.0 + phy0 / 2.0))\r
53                             + this.alpha\r
54                             * e / 2\r
55                             * Math.log((1 + e * sinPhy0)\r
56                             / (1 - e * sinPhy0));\r
57                   }\r
58                 \r
59                 \r
60                   override public function forward(p:ProjPoint):ProjPoint\r
61                   {\r
62                     var Sa1:Number = Math.log(Math.tan(Math.PI / 4.0 - p.y / 2.0));\r
63                     var Sa2:Number = this.e / 2.0\r
64                             * Math.log((1 + this.e * Math.sin(p.y))\r
65                             / (1 - this.e * Math.sin(p.y)));\r
66                     var S:Number = -this.alpha * (Sa1 + Sa2) + this.K;\r
67                 \r
68                         // spheric latitude\r
69                     var b:Number = 2.0 * (Math.atan(Math.exp(S)) - Math.PI / 4.0);\r
70                 \r
71                         // spheric longitude\r
72                     var I:Number = this.alpha * (p.x - this.lambda0);\r
73                 \r
74                         // psoeudo equatorial rotation\r
75                     var rotI:Number = Math.atan(Math.sin(I)\r
76                             / (Math.sin(this.b0) * Math.tan(b) +\r
77                                Math.cos(this.b0) * Math.cos(I)));\r
78                 \r
79                     var rotB:Number = Math.asin(Math.cos(this.b0) * Math.sin(b) -\r
80                                          Math.sin(this.b0) * Math.cos(b) * Math.cos(I));\r
81                 \r
82                     p.y = this.R / 2.0\r
83                             * Math.log((1 + Math.sin(rotB)) / (1 - Math.sin(rotB)))\r
84                             + this.y0;\r
85                     p.x = this.R * rotI + this.x0;\r
86                     return p;\r
87                   }\r
88                 \r
89                  override public function inverse(p:ProjPoint):ProjPoint\r
90                  {\r
91                     var Y:Number = p.x - this.x0;\r
92                     var X:Number = p.y - this.y0;\r
93                 \r
94                     var rotI:Number = Y / this.R;\r
95                     var rotB:Number = 2 * (Math.atan(Math.exp(X / this.R)) - Math.PI / 4.0);\r
96                 \r
97                     var b:Number = Math.asin(Math.cos(this.b0) * Math.sin(rotB)\r
98                             + Math.sin(this.b0) * Math.cos(rotB) * Math.cos(rotI));\r
99                     var I:Number = Math.atan(Math.sin(rotI)\r
100                             / (Math.cos(this.b0) * Math.cos(rotI) - Math.sin(this.b0)\r
101                             * Math.tan(rotB)));\r
102                 \r
103                     var lambda:Number = this.lambda0 + I / this.alpha;\r
104                 \r
105                     var S:Number = 0.0;\r
106                     var phy:Number = b;\r
107                     var prevPhy:Number = -1000.0;\r
108                     var iteration:Number = 0;\r
109                     while (Math.abs(phy - prevPhy) > 0.0000001)\r
110                     {\r
111                       if (++iteration > 20)\r
112                       {\r
113                         trace("omercFwdInfinity");\r
114                         return null;\r
115                       }\r
116                       //S = Math.log(Math.tan(Math.PI / 4.0 + phy / 2.0));\r
117                       S = 1.0\r
118                               / this.alpha\r
119                               * (Math.log(Math.tan(Math.PI / 4.0 + b / 2.0)) - this.K)\r
120                               + this.e\r
121                               * Math.log(Math.tan(Math.PI / 4.0\r
122                               + Math.asin(this.e * Math.sin(phy))\r
123                               / 2.0));\r
124                       prevPhy = phy;\r
125                       phy = 2.0 * Math.atan(Math.exp(S)) - Math.PI / 2.0;\r
126                     }\r
127                 \r
128                     p.x = lambda;\r
129                     p.y = phy;\r
130                     return p;\r
131                   }\r
132                 \r
133                 \r
134         }\r
135 }