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}