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 * @author Jan Peter Stotz
011 * @author Jason Huntley
012 */
013public class OsmMercator {
014
015    /**
016     * default tile size
017     */
018    public static final int DEFAUL_TILE_SIZE = 256;
019    /** maximum latitude (north) for mercator display */
020    public static final double MAX_LAT = 85.05112877980659;
021    /** minimum latitude (south) for mercator display */
022    public static final double MIN_LAT = -85.05112877980659;
023    /** equatorial earth radius for EPSG:3857 (Mercator) */
024    private static double EARTH_RADIUS = 6378137;
025
026    /**
027     * instance with tile size of 256 for easy conversions
028     */
029    public static final OsmMercator MERCATOR_256 = new OsmMercator();
030
031    /** tile size of the displayed tiles */
032    private int tileSize = DEFAUL_TILE_SIZE;
033
034    /**
035     * Creates instance with default tile size of 256
036     */
037    public OsmMercator() {
038    }
039
040    /**
041     * Creates instance with provided tile size.
042     * @param tileSize tile size in pixels
043     */
044    public OsmMercator(int tileSize) {
045        this.tileSize = tileSize;
046    }
047
048    public double radius(int aZoomlevel) {
049        return (tileSize * (1 << aZoomlevel)) / (2.0 * Math.PI);
050    }
051
052    /**
053     * Returns the absolut number of pixels in y or x, defined as: 2^Zoomlevel *
054     * tileSize where tileSize is the width of a tile in pixels
055     *
056     * @param aZoomlevel zoom level to request pixel data
057     * @return number of pixels
058     */
059    public int getMaxPixels(int aZoomlevel) {
060        return tileSize * (1 << aZoomlevel);
061    }
062
063    public int falseEasting(int aZoomlevel) {
064        return getMaxPixels(aZoomlevel) / 2;
065    }
066
067    public int falseNorthing(int aZoomlevel) {
068        return -1 * getMaxPixels(aZoomlevel) / 2;
069    }
070
071    /**
072     * Transform pixelspace to coordinates and get the distance.
073     *
074     * @param x1 the first x coordinate
075     * @param y1 the first y coordinate
076     * @param x2 the second x coordinate
077     * @param y2 the second y coordinate
078     *
079     * @param zoomLevel the zoom level
080     * @return the distance
081     */
082    public double getDistance(int x1, int y1, int x2, int y2, int zoomLevel) {
083        double la1 = yToLat(y1, zoomLevel);
084        double lo1 = xToLon(x1, zoomLevel);
085        double la2 = yToLat(y2, zoomLevel);
086        double lo2 = xToLon(x2, zoomLevel);
087
088        return getDistance(la1, lo1, la2, lo2);
089    }
090
091    /**
092     * Gets the distance using Spherical law of cosines.
093     *
094     * @param la1 the Latitude in degrees
095     * @param lo1 the Longitude in degrees
096     * @param la2 the Latitude from 2nd coordinate in degrees
097     * @param lo2 the Longitude from 2nd coordinate in degrees
098     * @return the distance
099     */
100    public double getDistance(double la1, double lo1, double la2, double lo2) {
101        double aStartLat = Math.toRadians(la1);
102        double aStartLong = Math.toRadians(lo1);
103        double aEndLat = Math.toRadians(la2);
104        double aEndLong = Math.toRadians(lo2);
105
106        double distance = Math.acos(Math.sin(aStartLat) * Math.sin(aEndLat)
107                + Math.cos(aStartLat) * Math.cos(aEndLat)
108                * Math.cos(aEndLong - aStartLong));
109
110        return EARTH_RADIUS * distance;
111    }
112
113    /**
114     * Transform longitude to pixelspace
115     *
116     * <p>
117     * Mathematical optimization<br>
118     * <code>
119     * x = radius(aZoomlevel) * toRadians(aLongitude) + falseEasting(aZoomLevel)<br>
120     * x = getMaxPixels(aZoomlevel) / (2 * PI) * (aLongitude * PI) / 180 + getMaxPixels(aZoomlevel) / 2<br>
121     * x = getMaxPixels(aZoomlevel) * aLongitude / 360 + 180 * getMaxPixels(aZoomlevel) / 360<br>
122     * x = getMaxPixels(aZoomlevel) * (aLongitude + 180) / 360<br>
123     * </code>
124     * </p>
125     *
126     * @param aLongitude
127     *            [-180..180]
128     * @param aZoomlevel zoom level
129     * @return [0..2^Zoomlevel*TILE_SIZE[
130     */
131    public double lonToX(double aLongitude, int aZoomlevel) {
132        int mp = getMaxPixels(aZoomlevel);
133        double x = (mp * (aLongitude + 180L)) / 360L;
134        return Math.min(x, mp);
135    }
136
137    /**
138     * Transforms latitude to pixelspace
139     * <p>
140     * Mathematical optimization<br>
141     * <code>
142     * log(u) := log((1.0 + sin(toRadians(aLat))) / (1.0 - sin(toRadians(aLat))<br>
143     *
144     * y = -1 * (radius(aZoomlevel) / 2 * log(u)))) - falseNorthing(aZoomlevel))<br>
145     * y = -1 * (getMaxPixel(aZoomlevel) / 2 * PI / 2 * log(u)) - -1 * getMaxPixel(aZoomLevel) / 2<br>
146     * y = getMaxPixel(aZoomlevel) / (-4 * PI) * log(u)) + getMaxPixel(aZoomLevel) / 2<br>
147     * y = getMaxPixel(aZoomlevel) * ((log(u) / (-4 * PI)) + 1/2)<br>
148     * </code>
149     * </p>
150     * @param aLat
151     *            [-90...90]
152     * @param aZoomlevel zoom level
153     * @return [0..2^Zoomlevel*TILE_SIZE[
154     */
155    public double latToY(double aLat, int aZoomlevel) {
156        if (aLat < MIN_LAT)
157            aLat = MIN_LAT;
158        else if (aLat > MAX_LAT)
159            aLat = MAX_LAT;
160        double sinLat = Math.sin(Math.toRadians(aLat));
161        double log = Math.log((1.0 + sinLat) / (1.0 - sinLat));
162        int mp = getMaxPixels(aZoomlevel);
163        double y = mp * (0.5 - (log / (4.0 * Math.PI)));
164        return Math.min(y, mp - 1);
165    }
166
167    /**
168     * Transforms pixel coordinate X to longitude
169     *
170     * <p>
171     * Mathematical optimization<br>
172     * <code>
173     * lon = toDegree((aX - falseEasting(aZoomlevel)) / radius(aZoomlevel))<br>
174     * lon = 180 / PI * ((aX - getMaxPixels(aZoomlevel) / 2) / getMaxPixels(aZoomlevel) / (2 * PI)<br>
175     * lon = 180 * ((aX - getMaxPixels(aZoomlevel) / 2) / getMaxPixels(aZoomlevel))<br>
176     * lon = 360 / getMaxPixels(aZoomlevel) * (aX - getMaxPixels(aZoomlevel) / 2)<br>
177     * lon = 360 * aX / getMaxPixels(aZoomlevel) - 180<br>
178     * </code>
179     * </p>
180     * @param aX
181     *            [0..2^Zoomlevel*TILE_WIDTH[
182     * @param aZoomlevel zoom level
183     * @return ]-180..180[
184     */
185    public double xToLon(int aX, int aZoomlevel) {
186        return ((360d * aX) / getMaxPixels(aZoomlevel)) - 180.0;
187    }
188
189    /**
190     * Transforms pixel coordinate Y to latitude
191     *
192     * @param aY
193     *            [0..2^Zoomlevel*TILE_WIDTH[
194     * @param aZoomlevel zoom level
195     * @return [MIN_LAT..MAX_LAT] is about [-85..85]
196     */
197    public double yToLat(int aY, int aZoomlevel) {
198        aY += falseNorthing(aZoomlevel);
199        double latitude = (Math.PI / 2) - (2 * Math.atan(Math.exp(-1.0 * aY / radius(aZoomlevel))));
200        return -1 * Math.toDegrees(latitude);
201    }
202
203}