001 package org.openstreetmap.gui.jmapviewer;
002
003 // License: GPL. Copyright 2007 by Tim Haussmann
004
005 /**
006 * This class implements the Mercator Projection as it is used by Openstreetmap
007 * (and google). It provides methods to translate coordinates from 'map space'
008 * into latitude and longitude (on the WGS84 ellipsoid) and vice versa. Map
009 * space is measured in pixels. The origin of the map space is the top left
010 * corner. The map space origin (0,0) has latitude ~85 and longitude -180
011 *
012 * @author Tim Haussmann
013 *
014 */
015
016 public class OsmMercator {
017
018 private static int TILE_SIZE = 256;
019 public static final double MAX_LAT = 85.05112877980659;
020 public static final double MIN_LAT = -85.05112877980659;
021 private static double EARTH_RADIUS = 6378137; // equatorial earth radius for EPSG:3857 (Mercator)
022
023 public static double radius(int aZoomlevel) {
024 return (TILE_SIZE * (1 << aZoomlevel)) / (2.0 * Math.PI);
025 }
026
027 /**
028 * Returns the absolut number of pixels in y or x, defined as: 2^Zoomlevel *
029 * TILE_WIDTH where TILE_WIDTH is the width of a tile in pixels
030 *
031 * @param aZoomlevel zoom level to request pixel data
032 * @return number of pixels
033 */
034 public static int getMaxPixels(int aZoomlevel) {
035 return TILE_SIZE * (1 << aZoomlevel);
036 }
037
038 public static int falseEasting(int aZoomlevel) {
039 return getMaxPixels(aZoomlevel) / 2;
040 }
041
042 public static int falseNorthing(int aZoomlevel) {
043 return (-1 * getMaxPixels(aZoomlevel) / 2);
044 }
045
046 /**
047 * Transform pixelspace to coordinates and get the distance.
048 *
049 * @param x1 the first x coordinate
050 * @param y1 the first y coordinate
051 * @param x2 the second x coordinate
052 * @param y2 the second y coordinate
053 *
054 * @param zoomLevel the zoom level
055 * @return the distance
056 * @author Jason Huntley
057 */
058 public static double getDistance(int x1, int y1, int x2, int y2, int zoomLevel) {
059 double la1 = YToLat(y1, zoomLevel);
060 double lo1 = XToLon(x1, zoomLevel);
061 double la2 = YToLat(y2, zoomLevel);
062 double lo2 = XToLon(x2, zoomLevel);
063
064 return getDistance(la1, lo1, la2, lo2);
065 }
066
067 /**
068 * Gets the distance using Spherical law of cosines.
069 *
070 * @param la1 the Latitude in degrees
071 * @param lo1 the Longitude in degrees
072 * @param la2 the Latitude from 2nd coordinate in degrees
073 * @param lo2 the Longitude from 2nd coordinate in degrees
074 * @return the distance
075 * @author Jason Huntley
076 */
077 public static double getDistance(double la1, double lo1, double la2, double lo2) {
078 double aStartLat = Math.toRadians(la1);
079 double aStartLong = Math.toRadians(lo1);
080 double aEndLat =Math.toRadians(la2);
081 double aEndLong = Math.toRadians(lo2);
082
083 double distance = Math.acos(Math.sin(aStartLat) * Math.sin(aEndLat)
084 + Math.cos(aStartLat) * Math.cos(aEndLat)
085 * Math.cos(aEndLong - aStartLong));
086
087 return (EARTH_RADIUS * distance);
088 }
089
090 /**
091 * Transform longitude to pixelspace
092 *
093 * <p>
094 * Mathematical optimization<br>
095 * <code>
096 * x = radius(aZoomlevel) * toRadians(aLongitude) + falseEasting(aZoomLevel)<br>
097 * x = getMaxPixels(aZoomlevel) / (2 * PI) * (aLongitude * PI) / 180 + getMaxPixels(aZoomlevel) / 2<br>
098 * x = getMaxPixels(aZoomlevel) * aLongitude / 360 + 180 * getMaxPixels(aZoomlevel) / 360<br>
099 * x = getMaxPixels(aZoomlevel) * (aLongitude + 180) / 360<br>
100 * </code>
101 * </p>
102 *
103 * @param aLongitude
104 * [-180..180]
105 * @return [0..2^Zoomlevel*TILE_SIZE[
106 * @author Jan Peter Stotz
107 */
108 public static int LonToX(double aLongitude, int aZoomlevel) {
109 int mp = getMaxPixels(aZoomlevel);
110 int x = (int) ((mp * (aLongitude + 180l)) / 360l);
111 x = Math.min(x, mp - 1);
112 return x;
113 }
114
115 /**
116 * Transforms latitude to pixelspace
117 * <p>
118 * Mathematical optimization<br>
119 * <code>
120 * log(u) := log((1.0 + sin(toRadians(aLat))) / (1.0 - sin(toRadians(aLat))<br>
121 *
122 * y = -1 * (radius(aZoomlevel) / 2 * log(u)))) - falseNorthing(aZoomlevel))<br>
123 * y = -1 * (getMaxPixel(aZoomlevel) / 2 * PI / 2 * log(u)) - -1 * getMaxPixel(aZoomLevel) / 2<br>
124 * y = getMaxPixel(aZoomlevel) / (-4 * PI) * log(u)) + getMaxPixel(aZoomLevel) / 2<br>
125 * y = getMaxPixel(aZoomlevel) * ((log(u) / (-4 * PI)) + 1/2)<br>
126 * </code>
127 * </p>
128 * @param aLat
129 * [-90...90]
130 * @return [0..2^Zoomlevel*TILE_SIZE[
131 * @author Jan Peter Stotz
132 */
133 public static int LatToY(double aLat, int aZoomlevel) {
134 if (aLat < MIN_LAT)
135 aLat = MIN_LAT;
136 else if (aLat > MAX_LAT)
137 aLat = MAX_LAT;
138 double sinLat = Math.sin(Math.toRadians(aLat));
139 double log = Math.log((1.0 + sinLat) / (1.0 - sinLat));
140 int mp = getMaxPixels(aZoomlevel);
141 int y = (int) (mp * (0.5 - (log / (4.0 * Math.PI))));
142 y = Math.min(y, mp - 1);
143 return y;
144 }
145
146 /**
147 * Transforms pixel coordinate X to longitude
148 *
149 * <p>
150 * Mathematical optimization<br>
151 * <code>
152 * lon = toDegree((aX - falseEasting(aZoomlevel)) / radius(aZoomlevel))<br>
153 * lon = 180 / PI * ((aX - getMaxPixels(aZoomlevel) / 2) / getMaxPixels(aZoomlevel) / (2 * PI)<br>
154 * lon = 180 * ((aX - getMaxPixels(aZoomlevel) / 2) / getMaxPixels(aZoomlevel))<br>
155 * lon = 360 / getMaxPixels(aZoomlevel) * (aX - getMaxPixels(aZoomlevel) / 2)<br>
156 * lon = 360 * aX / getMaxPixels(aZoomlevel) - 180<br>
157 * </code>
158 * </p>
159 * @param aX
160 * [0..2^Zoomlevel*TILE_WIDTH[
161 * @return ]-180..180[
162 * @author Jan Peter Stotz
163 */
164 public static double XToLon(int aX, int aZoomlevel) {
165 return ((360d * aX) / getMaxPixels(aZoomlevel)) - 180.0;
166 }
167
168 /**
169 * Transforms pixel coordinate Y to latitude
170 *
171 * @param aY
172 * [0..2^Zoomlevel*TILE_WIDTH[
173 * @return [MIN_LAT..MAX_LAT] is about [-85..85]
174 */
175 public static double YToLat(int aY, int aZoomlevel) {
176 aY += falseNorthing(aZoomlevel);
177 double latitude = (Math.PI / 2) - (2 * Math.atan(Math.exp(-1.0 * aY / radius(aZoomlevel))));
178 return -1 * Math.toDegrees(latitude);
179 }
180
181 }