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    }