1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
|
// Copyright 2018 Google Inc. All rights reserved.
//
// Licensed under the Apache License, Version 2.0 (the "License");
// you may not use this file except in compliance with the License.
// You may obtain a copy of the License at
//
// http://www.apache.org/licenses/LICENSE-2.0
//
// Unless required by applicable law or agreed to in writing, software
// distributed under the License is distributed on an "AS IS" BASIS,
// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
// See the License for the specific language governing permissions and
// limitations under the License.
package s2
import (
"math"
"github.com/golang/geo/r2"
"github.com/golang/geo/s1"
)
// Projection defines an interface for different ways of mapping between s2 and r2 Points.
// It can also define the coordinate wrapping behavior along each axis.
type Projection interface {
// Project converts a point on the sphere to a projected 2D point.
Project(p Point) r2.Point
// Unproject converts a projected 2D point to a point on the sphere.
//
// If wrapping is defined for a given axis (see below), then this method
// should accept any real number for the corresponding coordinate.
Unproject(p r2.Point) Point
// FromLatLng is a convenience function equivalent to Project(LatLngToPoint(ll)),
// but the implementation is more efficient.
FromLatLng(ll LatLng) r2.Point
// ToLatLng is a convenience function equivalent to LatLngFromPoint(Unproject(p)),
// but the implementation is more efficient.
ToLatLng(p r2.Point) LatLng
// Interpolate returns the point obtained by interpolating the given
// fraction of the distance along the line from A to B.
// Fractions < 0 or > 1 result in extrapolation instead.
Interpolate(f float64, a, b r2.Point) r2.Point
// WrapDistance reports the coordinate wrapping distance along each axis.
// If this value is non-zero for a given axis, the coordinates are assumed
// to "wrap" with the given period. For example, if WrapDistance.Y == 360
// then (x, y) and (x, y + 360) should map to the same Point.
//
// This information is used to ensure that edges takes the shortest path
// between two given points. For example, if coordinates represent
// (latitude, longitude) pairs in degrees and WrapDistance().Y == 360,
// then the edge (5:179, 5:-179) would be interpreted as spanning 2 degrees
// of longitude rather than 358 degrees.
//
// If a given axis does not wrap, its WrapDistance should be set to zero.
WrapDistance() r2.Point
}
// PlateCarreeProjection defines the "plate carree" (square plate) projection,
// which converts points on the sphere to (longitude, latitude) pairs.
// Coordinates can be scaled so that they represent radians, degrees, etc, but
// the projection is always centered around (latitude=0, longitude=0).
//
// Note that (x, y) coordinates are backwards compared to the usual (latitude,
// longitude) ordering, in order to match the usual convention for graphs in
// which "x" is horizontal and "y" is vertical.
type PlateCarreeProjection struct {
xWrap float64
toRadians float64 // Multiplier to convert coordinates to radians.
fromRadians float64 // Multiplier to convert coordinates from radians.
}
// NewPlateCarreeProjection constructs a plate carree projection where the
// x-coordinates (lng) span [-xScale, xScale] and the y coordinates (lat)
// span [-xScale/2, xScale/2]. For example if xScale==180 then the x
// range is [-180, 180] and the y range is [-90, 90].
//
// By default coordinates are expressed in radians, i.e. the x range is
// [-Pi, Pi] and the y range is [-Pi/2, Pi/2].
func NewPlateCarreeProjection(xScale float64) Projection {
return &PlateCarreeProjection{
xWrap: 2 * xScale,
toRadians: math.Pi / xScale,
fromRadians: xScale / math.Pi,
}
}
// Project converts a point on the sphere to a projected 2D point.
func (p *PlateCarreeProjection) Project(pt Point) r2.Point {
return p.FromLatLng(LatLngFromPoint(pt))
}
// Unproject converts a projected 2D point to a point on the sphere.
func (p *PlateCarreeProjection) Unproject(pt r2.Point) Point {
return PointFromLatLng(p.ToLatLng(pt))
}
// FromLatLng returns the LatLng projected into an R2 Point.
func (p *PlateCarreeProjection) FromLatLng(ll LatLng) r2.Point {
return r2.Point{
X: p.fromRadians * ll.Lng.Radians(),
Y: p.fromRadians * ll.Lat.Radians(),
}
}
// ToLatLng returns the LatLng projected from the given R2 Point.
func (p *PlateCarreeProjection) ToLatLng(pt r2.Point) LatLng {
return LatLng{
Lat: s1.Angle(p.toRadians * pt.Y),
Lng: s1.Angle(p.toRadians * math.Remainder(pt.X, p.xWrap)),
}
}
// Interpolate returns the point obtained by interpolating the given
// fraction of the distance along the line from A to B.
func (p *PlateCarreeProjection) Interpolate(f float64, a, b r2.Point) r2.Point {
return a.Mul(1 - f).Add(b.Mul(f))
}
// WrapDistance reports the coordinate wrapping distance along each axis.
func (p *PlateCarreeProjection) WrapDistance() r2.Point {
return r2.Point{p.xWrap, 0}
}
// MercatorProjection defines the spherical Mercator projection. Google Maps
// uses this projection together with WGS84 coordinates, in which case it is
// known as the "Web Mercator" projection (see Wikipedia). This class makes
// no assumptions regarding the coordinate system of its input points, but
// simply applies the spherical Mercator projection to them.
//
// The Mercator projection is finite in width (x) but infinite in height (y).
// "x" corresponds to longitude, and spans a finite range such as [-180, 180]
// (with coordinate wrapping), while "y" is a function of latitude and spans
// an infinite range. (As "y" coordinates get larger, points get closer to
// the north pole but never quite reach it.) The north and south poles have
// infinite "y" values. (Note that this will cause problems if you tessellate
// a Mercator edge where one endpoint is a pole. If you need to do this, clip
// the edge first so that the "y" coordinate is no more than about 5 * maxX.)
type MercatorProjection struct {
xWrap float64
toRadians float64 // Multiplier to convert coordinates to radians.
fromRadians float64 // Multiplier to convert coordinates from radians.
}
// NewMercatorProjection constructs a Mercator projection with the given maximum
// longitude axis value corresponding to a range of [-maxLng, maxLng].
// The horizontal and vertical axes are scaled equally.
func NewMercatorProjection(maxLng float64) Projection {
return &MercatorProjection{
xWrap: 2 * maxLng,
toRadians: math.Pi / maxLng,
fromRadians: maxLng / math.Pi,
}
}
// Project converts a point on the sphere to a projected 2D point.
func (p *MercatorProjection) Project(pt Point) r2.Point {
return p.FromLatLng(LatLngFromPoint(pt))
}
// Unproject converts a projected 2D point to a point on the sphere.
func (p *MercatorProjection) Unproject(pt r2.Point) Point {
return PointFromLatLng(p.ToLatLng(pt))
}
// FromLatLng returns the LatLng projected into an R2 Point.
func (p *MercatorProjection) FromLatLng(ll LatLng) r2.Point {
// This formula is more accurate near zero than the log(tan()) version.
// Note that latitudes of +/- 90 degrees yield "y" values of +/- infinity.
sinPhi := math.Sin(float64(ll.Lat))
y := 0.5 * math.Log((1+sinPhi)/(1-sinPhi))
return r2.Point{p.fromRadians * float64(ll.Lng), p.fromRadians * y}
}
// ToLatLng returns the LatLng projected from the given R2 Point.
func (p *MercatorProjection) ToLatLng(pt r2.Point) LatLng {
// This formula is more accurate near zero than the atan(exp()) version.
x := p.toRadians * math.Remainder(pt.X, p.xWrap)
k := math.Exp(2 * p.toRadians * pt.Y)
var y float64
if math.IsInf(k, 0) {
y = math.Pi / 2
} else {
y = math.Asin((k - 1) / (k + 1))
}
return LatLng{s1.Angle(y), s1.Angle(x)}
}
// Interpolate returns the point obtained by interpolating the given
// fraction of the distance along the line from A to B.
func (p *MercatorProjection) Interpolate(f float64, a, b r2.Point) r2.Point {
return a.Mul(1 - f).Add(b.Mul(f))
}
// WrapDistance reports the coordinate wrapping distance along each axis.
func (p *MercatorProjection) WrapDistance() r2.Point {
return r2.Point{p.xWrap, 0}
}
|
