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