Understand shapefiles in different projections (notably OSGB)
[potlatch2.git] / com / gradoservice / proj4as / Proj4as.as
1 package com.gradoservice.proj4as
2 {
3         import com.gradoservice.proj4as.proj.AbstractProjProjection;
4         
5         public class Proj4as
6         {
7                 
8                 static public const defaultDatum:String = 'WGS84';
9                 static public const WGS84:ProjProjection = new ProjProjection('WGS84');
10
11                 
12                 public function Proj4as()
13                 {
14                 }
15                 
16                 public function transform(source:ProjProjection, dest:ProjProjection, point:ProjPoint):ProjPoint
17                 {
18                 
19                 if (source==null || dest==null || point==null)
20                 {
21                         trace("Parametrs not created!");
22                         return null;  
23                 }
24                 
25                 if (!source.readyToUse || !dest.readyToUse) 
26                 {
27             trace("Proj4js initialization for "+source.srsCode+" not yet complete");
28                 return point;
29                 }
30         
31         // Workaround for Spherical Mercator
32         if ((source.srsProjNumber =="900913" && dest.datumCode != "WGS84") ||
33             (dest.srsProjNumber == "900913" && source.datumCode != "WGS84")) {
34             var wgs84:ProjProjection = WGS84;
35             this.transform(source, wgs84, point);
36             source = wgs84;
37         }
38
39         // Transform source points to long/lat, if they aren't already.
40         if ( source.projName=="longlat") {
41             point.x *= ProjConstants.D2R;  // convert degrees to radians
42             point.y *= ProjConstants.D2R;
43         } else {
44             if (source.to_meter) {
45                 point.x *= source.to_meter;
46                 point.y *= source.to_meter;
47             }
48             source.inverse(point); // Convert Cartesian to longlat
49         }
50
51         // Adjust for the prime meridian if necessary
52         if (source.from_greenwich) { 
53             point.x += source.from_greenwich; 
54         }
55
56         // Convert datums if needed, and if possible.
57         point = this.datum_transform( source.datum, dest.datum, point );
58
59         // Adjust for the prime meridian if necessary
60         if (dest.from_greenwich) {
61             point.x -= dest.from_greenwich;
62         }
63
64         if( dest.projName=="longlat" ) {             
65             // convert radians to decimal degrees
66             point.x *= ProjConstants.R2D;
67             point.y *= ProjConstants.R2D;
68         } else  {               // else project
69             dest.forward(point);
70             if (dest.to_meter) {
71                 point.x /= dest.to_meter;
72                 point.y /= dest.to_meter;
73             }
74         }
75         return point;                   
76                 }
77                 
78                 
79     public function datum_transform( source:Datum, dest:Datum, point:ProjPoint ):ProjPoint 
80     {
81
82       // Short cut if the datums are identical.
83       if( source.compare_datums( dest ) ) {
84           return point; // in this case, zero is sucess,
85                     // whereas cs_compare_datums returns 1 to indicate TRUE
86                     // confusing, should fix this
87       }
88
89       // Explicitly skip datum transform by setting 'datum=none' as parameter for either source or dest
90       if( source.datum_type == ProjConstants.PJD_NODATUM
91           || dest.datum_type == ProjConstants.PJD_NODATUM) {
92           return point;
93       }
94
95       // If this datum requires grid shifts, then apply it to geodetic coordinates.
96       if( source.datum_type == ProjConstants.PJD_GRIDSHIFT )
97       {
98        trace("ERROR: Grid shift transformations are not implemented yet.");
99         /*
100           pj_apply_gridshift( pj_param(source.params,"snadgrids").s, 0,
101                               point_count, point_offset, x, y, z );
102           CHECK_RETURN;
103
104           src_a = SRS_WGS84_SEMIMAJOR;
105           src_es = 0.006694379990;
106         */
107       }
108
109       if( dest.datum_type == ProjConstants.PJD_GRIDSHIFT )
110       {
111         trace("ERROR: Grid shift transformations are not implemented yet.");
112         /*
113           dst_a = ;
114           dst_es = 0.006694379990;
115         */
116       }
117
118       // Do we need to go through geocentric coordinates?
119       if( source.es != dest.es || source.a != dest.a
120           || source.datum_type == ProjConstants.PJD_3PARAM
121           || source.datum_type == ProjConstants.PJD_7PARAM
122           || dest.datum_type == ProjConstants.PJD_3PARAM
123           || dest.datum_type == ProjConstants.PJD_7PARAM)
124       {
125
126         // Convert to geocentric coordinates.
127         source.geodetic_to_geocentric( point );
128         // CHECK_RETURN;
129
130         // Convert between datums
131         if( source.datum_type == ProjConstants.PJD_3PARAM || source.datum_type == ProjConstants.PJD_7PARAM ) {
132           source.geocentric_to_wgs84(point);
133           // CHECK_RETURN;
134         }
135
136         if( dest.datum_type == ProjConstants.PJD_3PARAM || dest.datum_type == ProjConstants.PJD_7PARAM ) {
137           dest.geocentric_from_wgs84(point);
138           // CHECK_RETURN;
139         }
140
141         // Convert back to geodetic coordinates
142         dest.geocentric_to_geodetic( point );
143           // CHECK_RETURN;
144       }
145
146       // Apply grid shift to destination if required
147       if( dest.datum_type == ProjConstants.PJD_GRIDSHIFT )
148       {
149         trace("ERROR: Grid shift transformations are not implemented yet.");
150         // pj_apply_gridshift( pj_param(dest.params,"snadgrids").s, 1, point);
151         // CHECK_RETURN;
152       }
153       return point;
154     }           
155
156         }
157 }