1
+ ( function ( ) {
2
+ var gju = this . gju = { } ;
3
+
4
+ // Export the geojson object for **CommonJS**
5
+ if ( typeof module !== 'undefined' && module . exports ) {
6
+ module . exports = gju ;
7
+ }
8
+
9
+ // adapted from http://www.kevlindev.com/gui/math/intersection/Intersection.js
10
+ gju . lineStringsIntersect = function ( l1 , l2 ) {
11
+ var intersects = [ ] ;
12
+ for ( var i = 0 ; i <= l1 . coordinates . length - 2 ; ++ i ) {
13
+ for ( var j = 0 ; j <= l2 . coordinates . length - 2 ; ++ j ) {
14
+ var a1 = {
15
+ x : l1 . coordinates [ i ] [ 1 ] ,
16
+ y : l1 . coordinates [ i ] [ 0 ]
17
+ } ,
18
+ a2 = {
19
+ x : l1 . coordinates [ i + 1 ] [ 1 ] ,
20
+ y : l1 . coordinates [ i + 1 ] [ 0 ]
21
+ } ,
22
+ b1 = {
23
+ x : l2 . coordinates [ j ] [ 1 ] ,
24
+ y : l2 . coordinates [ j ] [ 0 ]
25
+ } ,
26
+ b2 = {
27
+ x : l2 . coordinates [ j + 1 ] [ 1 ] ,
28
+ y : l2 . coordinates [ j + 1 ] [ 0 ]
29
+ } ,
30
+ ua_t = ( b2 . x - b1 . x ) * ( a1 . y - b1 . y ) - ( b2 . y - b1 . y ) * ( a1 . x - b1 . x ) ,
31
+ ub_t = ( a2 . x - a1 . x ) * ( a1 . y - b1 . y ) - ( a2 . y - a1 . y ) * ( a1 . x - b1 . x ) ,
32
+ u_b = ( b2 . y - b1 . y ) * ( a2 . x - a1 . x ) - ( b2 . x - b1 . x ) * ( a2 . y - a1 . y ) ;
33
+ if ( u_b != 0 ) {
34
+ var ua = ua_t / u_b ,
35
+ ub = ub_t / u_b ;
36
+ if ( 0 <= ua && ua <= 1 && 0 <= ub && ub <= 1 ) {
37
+ intersects . push ( {
38
+ 'type' : 'Point' ,
39
+ 'coordinates' : [ a1 . x + ua * ( a2 . x - a1 . x ) , a1 . y + ua * ( a2 . y - a1 . y ) ]
40
+ } ) ;
41
+ }
42
+ }
43
+ }
44
+ }
45
+ if ( intersects . length == 0 ) intersects = false ;
46
+ return intersects ;
47
+ }
48
+
49
+ // adapted from http://jsfromhell.com/math/is-point-in-poly
50
+ gju . pointInPolygon = function ( point , polygon ) {
51
+ var x = point . coordinates [ 1 ] ,
52
+ y = point . coordinates [ 0 ] ,
53
+ poly = polygon . coordinates [ 0 ] ; //TODO: support polygons with holes
54
+ for ( var c = false , i = - 1 , l = poly . length , j = l - 1 ; ++ i < l ; j = i ) {
55
+ var px = poly [ i ] [ 1 ] ,
56
+ py = poly [ i ] [ 0 ] ,
57
+ jx = poly [ j ] [ 1 ] ,
58
+ jy = poly [ j ] [ 0 ] ;
59
+ if ( ( ( py <= y && y < jy ) || ( jy <= y && y < py ) ) && ( x < ( jx - px ) * ( y - py ) / ( jy - py ) + px ) ) {
60
+ c = [ point ] ;
61
+ }
62
+ }
63
+ return c ;
64
+ }
65
+
66
+ gju . numberToRadius = function ( number ) {
67
+ return number * Math . PI / 180 ;
68
+ }
69
+
70
+ gju . numberToDegree = function ( number ) {
71
+ return number * 180 / Math . PI ;
72
+ }
73
+
74
+ // written with help from @tautologe
75
+ gju . drawCircle = function ( radiusInMeters , centerPoint , steps ) {
76
+ var center = [ centerPoint . coordinates [ 1 ] , centerPoint . coordinates [ 0 ] ] ,
77
+ dist = ( radiusInMeters / 1000 ) / 6371 ,
78
+ // convert meters to radiant
79
+ radCenter = [ gju . numberToRadius ( center [ 0 ] ) , gju . numberToRadius ( center [ 1 ] ) ] ,
80
+ steps = steps || 15 ,
81
+ // 15 sided circle
82
+ poly = [ [ center [ 0 ] , center [ 1 ] ] ] ;
83
+ for ( var i = 0 ; i < steps ; i ++ ) {
84
+ var brng = 2 * Math . PI * i / steps ;
85
+ var lat = Math . asin ( Math . sin ( radCenter [ 0 ] ) * Math . cos ( dist )
86
+ + Math . cos ( radCenter [ 0 ] ) * Math . sin ( dist ) * Math . cos ( brng ) ) ;
87
+ var lng = radCenter [ 1 ] + Math . atan2 ( Math . sin ( brng ) * Math . sin ( dist ) * Math . cos ( radCenter [ 0 ] ) ,
88
+ Math . cos ( dist ) - Math . sin ( radCenter [ 0 ] ) * Math . sin ( lat ) ) ;
89
+ poly [ i ] = [ ] ;
90
+ poly [ i ] [ 1 ] = gju . numberToDegree ( lat ) ;
91
+ poly [ i ] [ 0 ] = gju . numberToDegree ( lng ) ;
92
+ }
93
+ return {
94
+ "type" : "Polygon" ,
95
+ "coordinates" : [ poly ]
96
+ } ;
97
+ }
98
+
99
+ // assumes rectangle starts at lower left point
100
+ gju . rectangleCentroid = function ( rectangle ) {
101
+ var bbox = rectangle . coordinates [ 0 ] ;
102
+ var xmin = bbox [ 0 ] [ 0 ] ,
103
+ ymin = bbox [ 0 ] [ 1 ] ,
104
+ xmax = bbox [ 2 ] [ 0 ] ,
105
+ ymax = bbox [ 2 ] [ 1 ] ;
106
+ var xwidth = xmax - xmin ;
107
+ var ywidth = ymax - ymin ;
108
+ return {
109
+ 'type' : 'Point' ,
110
+ 'coordinates' : [ xmin + xwidth / 2 , ymin + ywidth / 2 ]
111
+ } ;
112
+ }
113
+
114
+ // from http://www.movable-type.co.uk/scripts/latlong.html
115
+ gju . pointDistance = function ( pt1 , pt2 ) {
116
+ var lon1 = pt1 . coordinates [ 0 ] ,
117
+ lat1 = pt1 . coordinates [ 1 ] ,
118
+ lon2 = pt2 . coordinates [ 0 ] ,
119
+ lat2 = pt2 . coordinates [ 1 ] ,
120
+ dLat = gju . numberToRadius ( lat2 - lat1 ) ,
121
+ dLon = gju . numberToRadius ( lon2 - lon1 ) ,
122
+ a = Math . pow ( Math . sin ( dLat / 2 ) , 2 ) + Math . cos ( gju . numberToRadius ( lat1 ) )
123
+ * Math . cos ( gju . numberToRadius ( lat2 ) ) * Math . pow ( Math . sin ( dLon / 2 ) , 2 ) ,
124
+ c = 2 * Math . atan2 ( Math . sqrt ( a ) , Math . sqrt ( 1 - a ) ) ;
125
+ return ( 6371 * c ) * 1000 ; // returns meters
126
+ } ,
127
+
128
+ // checks if geometry lies entirely within a circle
129
+ // works with Point, LineString, Polygon
130
+ gju . geometryWithinRadius = function ( geometry , center , radius ) {
131
+ if ( geometry . type == 'Point' ) {
132
+ return gju . pointDistance ( geometry , center ) <= radius ;
133
+ } else if ( geometry . type == 'LineString' || geometry . type == 'Polygon' ) {
134
+ var point = { } ;
135
+ var coordinates ;
136
+ if ( geometry . type == 'Polygon' ) {
137
+ // it's enough to check the exterior ring of the Polygon
138
+ coordinates = geometry . coordinates [ 0 ] ;
139
+ } else {
140
+ coordinates = geometry . coordinates ;
141
+ }
142
+ for ( var i in coordinates ) {
143
+ point . coordinates = coordinates [ i ] ;
144
+ if ( gju . pointDistance ( point , center ) > radius ) {
145
+ return false ;
146
+ }
147
+ }
148
+ }
149
+ return true ;
150
+ }
151
+
152
+ // adapted from http://paulbourke.net/geometry/polyarea/javascript.txt
153
+ gju . area = function ( polygon ) {
154
+ var area = 0 ;
155
+ // TODO: polygon holes at coordinates[1]
156
+ var points = polygon . coordinates [ 0 ] ;
157
+ var j = points . length - 1 ;
158
+ var p1 , p2 ;
159
+
160
+ for ( var i = 0 ; i < points . length ; j = i ++ ) {
161
+ var p1 = {
162
+ x : points [ i ] [ 1 ] ,
163
+ y : points [ i ] [ 0 ]
164
+ } ;
165
+ var p2 = {
166
+ x : points [ j ] [ 1 ] ,
167
+ y : points [ j ] [ 0 ]
168
+ } ;
169
+ area += p1 . x * p2 . y ;
170
+ area -= p1 . y * p2 . x ;
171
+ }
172
+
173
+ area /= 2 ;
174
+ return area ;
175
+ } ,
176
+
177
+ // adapted from http://paulbourke.net/geometry/polyarea/javascript.txt
178
+ gju . centroid = function ( polygon ) {
179
+ var f , x = 0 ,
180
+ y = 0 ;
181
+ // TODO: polygon holes at coordinates[1]
182
+ var points = polygon . coordinates [ 0 ] ;
183
+ var j = points . length - 1 ;
184
+ var p1 , p2 ;
185
+
186
+ for ( var i = 0 ; i < points . length ; j = i ++ ) {
187
+ var p1 = {
188
+ x : points [ i ] [ 1 ] ,
189
+ y : points [ i ] [ 0 ]
190
+ } ;
191
+ var p2 = {
192
+ x : points [ j ] [ 1 ] ,
193
+ y : points [ j ] [ 0 ]
194
+ } ;
195
+ f = p1 . x * p2 . y - p2 . x * p1 . y ;
196
+ x += ( p1 . x + p2 . x ) * f ;
197
+ y += ( p1 . y + p2 . y ) * f ;
198
+ }
199
+
200
+ f = gju . area ( polygon ) * 6 ;
201
+ return {
202
+ 'type' : 'Point' ,
203
+ 'coordinates' : [ y / f , x / f ]
204
+ } ;
205
+ } ,
206
+
207
+ gju . simplify = function ( source , kink ) { /* source[] array of geojson points */
208
+ /* kink in metres, kinks above this depth kept */
209
+ /* kink depth is the height of the triangle abc where a-b and b-c are two consecutive line segments */
210
+ kink = kink || 20 ;
211
+ source = source . map ( function ( o ) {
212
+ return {
213
+ lng : o . coordinates [ 0 ] ,
214
+ lat : o . coordinates [ 1 ]
215
+ }
216
+ } ) ;
217
+
218
+ var n_source , n_stack , n_dest , start , end , i , sig ;
219
+ var dev_sqr , max_dev_sqr , band_sqr ;
220
+ var x12 , y12 , d12 , x13 , y13 , d13 , x23 , y23 , d23 ;
221
+ var F = ( Math . PI / 180.0 ) * 0.5 ;
222
+ var index = new Array ( ) ; /* aray of indexes of source points to include in the reduced line */
223
+ var sig_start = new Array ( ) ; /* indices of start & end of working section */
224
+ var sig_end = new Array ( ) ;
225
+
226
+ /* check for simple cases */
227
+
228
+ if ( source . length < 3 ) return ( source ) ; /* one or two points */
229
+
230
+ /* more complex case. initialize stack */
231
+
232
+ n_source = source . length ;
233
+ band_sqr = kink * 360.0 / ( 2.0 * Math . PI * 6378137.0 ) ; /* Now in degrees */
234
+ band_sqr *= band_sqr ;
235
+ n_dest = 0 ;
236
+ sig_start [ 0 ] = 0 ;
237
+ sig_end [ 0 ] = n_source - 1 ;
238
+ n_stack = 1 ;
239
+
240
+ /* while the stack is not empty ... */
241
+ while ( n_stack > 0 ) {
242
+
243
+ /* ... pop the top-most entries off the stacks */
244
+
245
+ start = sig_start [ n_stack - 1 ] ;
246
+ end = sig_end [ n_stack - 1 ] ;
247
+ n_stack -- ;
248
+
249
+ if ( ( end - start ) > 1 ) { /* any intermediate points ? */
250
+
251
+ /* ... yes, so find most deviant intermediate point to
252
+ either side of line joining start & end points */
253
+
254
+ x12 = ( source [ end ] . lng ( ) - source [ start ] . lng ( ) ) ;
255
+ y12 = ( source [ end ] . lat ( ) - source [ start ] . lat ( ) ) ;
256
+ if ( Math . abs ( x12 ) > 180.0 ) x12 = 360.0 - Math . abs ( x12 ) ;
257
+ x12 *= Math . cos ( F * ( source [ end ] . lat ( ) + source [ start ] . lat ( ) ) ) ; /* use avg lat to reduce lng */
258
+ d12 = ( x12 * x12 ) + ( y12 * y12 ) ;
259
+
260
+ for ( i = start + 1 , sig = start , max_dev_sqr = - 1.0 ; i < end ; i ++ ) {
261
+
262
+ x13 = source [ i ] . lng ( ) - source [ start ] . lng ( ) ;
263
+ y13 = source [ i ] . lat ( ) - source [ start ] . lat ( ) ;
264
+ if ( Math . abs ( x13 ) > 180.0 ) x13 = 360.0 - Math . abs ( x13 ) ;
265
+ x13 *= Math . cos ( F * ( source [ i ] . lat ( ) + source [ start ] . lat ( ) ) ) ;
266
+ d13 = ( x13 * x13 ) + ( y13 * y13 ) ;
267
+
268
+ x23 = source [ i ] . lng ( ) - source [ end ] . lng ( ) ;
269
+ y23 = source [ i ] . lat ( ) - source [ end ] . lat ( ) ;
270
+ if ( Math . abs ( x23 ) > 180.0 ) x23 = 360.0 - Math . abs ( x23 ) ;
271
+ x23 *= Math . cos ( F * ( source [ i ] . lat ( ) + source [ end ] . lat ( ) ) ) ;
272
+ d23 = ( x23 * x23 ) + ( y23 * y23 ) ;
273
+
274
+ if ( d13 >= ( d12 + d23 ) ) dev_sqr = d23 ;
275
+ else if ( d23 >= ( d12 + d13 ) ) dev_sqr = d13 ;
276
+ else dev_sqr = ( x13 * y12 - y13 * x12 ) * ( x13 * y12 - y13 * x12 ) / d12 ; // solve triangle
277
+ if ( dev_sqr > max_dev_sqr ) {
278
+ sig = i ;
279
+ max_dev_sqr = dev_sqr ;
280
+ }
281
+ }
282
+
283
+ if ( max_dev_sqr < band_sqr ) { /* is there a sig. intermediate point ? */
284
+ /* ... no, so transfer current start point */
285
+ index [ n_dest ] = start ;
286
+ n_dest ++ ;
287
+ } else { /* ... yes, so push two sub-sections on stack for further processing */
288
+ n_stack ++ ;
289
+ sig_start [ n_stack - 1 ] = sig ;
290
+ sig_end [ n_stack - 1 ] = end ;
291
+ n_stack ++ ;
292
+ sig_start [ n_stack - 1 ] = start ;
293
+ sig_end [ n_stack - 1 ] = sig ;
294
+ }
295
+ } else { /* ... no intermediate points, so transfer current start point */
296
+ index [ n_dest ] = start ;
297
+ n_dest ++ ;
298
+ }
299
+ }
300
+
301
+ /* transfer last point */
302
+ index [ n_dest ] = n_source - 1 ;
303
+ n_dest ++ ;
304
+
305
+ /* make return array */
306
+ var r = new Array ( ) ;
307
+ for ( var i = 0 ; i < n_dest ; i ++ )
308
+ r . push ( source [ index [ i ] ] ) ;
309
+
310
+ return r . map ( function ( o ) {
311
+ return {
312
+ type : "Point" ,
313
+ coordinates : [ o . lng , o . lat ]
314
+ }
315
+ } ) ;
316
+ }
317
+
318
+ // http://www.movable-type.co.uk/scripts/latlong.html#destPoint
319
+ gju . destinationPoint = function ( pt , brng , dist ) {
320
+ dist = dist / 6371 ; // convert dist to angular distance in radians
321
+ brng = gju . numberToRadius ( brng ) ;
322
+
323
+ var lat1 = gju . numberToRadius ( pt . coordinates [ 0 ] ) ;
324
+ var lon1 = gju . numberToRadius ( pt . coordinates [ 1 ] ) ;
325
+
326
+ var lat2 = Math . asin ( Math . sin ( lat1 ) * Math . cos ( dist ) +
327
+ Math . cos ( lat1 ) * Math . sin ( dist ) * Math . cos ( brng ) ) ;
328
+ var lon2 = lon1 + Math . atan2 ( Math . sin ( brng ) * Math . sin ( dist ) * Math . cos ( lat1 ) ,
329
+ Math . cos ( dist ) - Math . sin ( lat1 ) * Math . sin ( lat2 ) ) ;
330
+ lon2 = ( lon2 + 3 * Math . PI ) % ( 2 * Math . PI ) - Math . PI ; // normalise to -180..+180º
331
+
332
+ return {
333
+ 'type' : 'Point' ,
334
+ 'coordinates' : [ gju . numberToDegree ( lat2 ) , gju . numberToDegree ( lon2 ) ]
335
+ } ;
336
+ } ;
337
+
338
+ } ) ( ) ;
0 commit comments