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