Преобразовать точку широты / долготы в пиксели (x,y) на проекции меркатора
Я пытаюсь преобразовать точку широты / долготы в точку 2d, чтобы я мог отобразить ее на изображении мира, которое является проекцией меркатора.
Я видел различные способы сделать это и несколько вопросов о переполнении стека - я опробовал разные фрагменты кода, и хотя я получаю правильную долготу в пикселях, широта всегда отключена - кажется, что она становится более разумной.
Мне нужна формула, чтобы принять во внимание размер изображения, ширину и т. Д.
Я пробовал этот кусок кода:
double minLat = -85.05112878;
double minLong = -180;
double maxLat = 85.05112878;
double maxLong = 180;
// Map image size (in points)
double mapHeight = 768.0;
double mapWidth = 991.0;
// Determine the map scale (points per degree)
double xScale = mapWidth/ (maxLong - minLong);
double yScale = mapHeight / (maxLat - minLat);
// position of map image for point
double x = (lon - minLong) * xScale;
double y = - (lat + minLat) * yScale;
System.out.println("final coords: " + x + " " + y);
В примере, который я пытаюсь использовать, широта примерно равна 30px. Любая помощь или совет?
Обновить
На основании этого вопроса: Lat / lon to xy
Я пытался использовать предоставленный код, но у меня все еще есть некоторые проблемы с преобразованием широты, долгота в порядке.
int mapWidth = 991;
int mapHeight = 768;
double mapLonLeft = -180;
double mapLonRight = 180;
double mapLonDelta = mapLonRight - mapLonLeft;
double mapLatBottom = -85.05112878;
double mapLatBottomDegree = mapLatBottom * Math.PI / 180;
double worldMapWidth = ((mapWidth / mapLonDelta) * 360) / (2 * Math.PI);
double mapOffsetY = (worldMapWidth / 2 * Math.log((1 + Math.sin(mapLatBottomDegree)) / (1 - Math.sin(mapLatBottomDegree))));
double x = (lon - mapLonLeft) * (mapWidth / mapLonDelta);
double y = 0.1;
if (lat < 0) {
lat = lat * Math.PI / 180;
y = mapHeight - ((worldMapWidth / 2 * Math.log((1 + Math.sin(lat)) / (1 - Math.sin(lat)))) - mapOffsetY);
} else if (lat > 0) {
lat = lat * Math.PI / 180;
lat = lat * -1;
y = mapHeight - ((worldMapWidth / 2 * Math.log((1 + Math.sin(lat)) / (1 - Math.sin(lat)))) - mapOffsetY);
System.out.println("y before minus: " + y);
y = mapHeight - y;
} else {
y = mapHeight / 2;
}
System.out.println(x);
System.out.println(y);
При использовании исходного кода, если значение широты положительное, оно вернуло отрицательную точку, поэтому я немного изменил его и проверил с крайними широтами - которые должны быть точкой 0 и точкой 766, это работает нормально. Однако, когда я пытаюсь использовать другое значение широты, например, 58,07 (к северу от Великобритании), оно отображается как север Испании.
7 ответов
Проекция карты Меркатора является особым предельным случаем конформной проекции Ламберта с коническим отображением с экватором в качестве единой стандартной параллели. Все остальные параллели широты являются прямыми линиями, а меридианы также являются прямыми линиями под прямыми углами к экватору, которые расположены на одинаковом расстоянии. Это основа для поперечной и косой форм проекции. Он мало используется для картографирования, но почти универсален для навигационных карт. Помимо того, что оно конформно, оно обладает особым свойством, заключающимся в том, что прямые линии, нарисованные на нем, являются линиями постоянного значения. Таким образом, навигаторы могут получить свой курс из угла, который прямая линия курса составляет с меридианами. [1]
Формулы для получения проецируемых координат на восток и север от сферической широты φ и долготы λ:
E = FE + R (λ – λₒ)
N = FN + R ln[tan(π/4 + φ/2)]
где λ O - долгота естественного происхождения, а FE и FN - ложное восточное и ложное северное направление. В сферическом Меркаторе эти значения фактически не используются, поэтому вы можете упростить формулу
Пример псевдокода, поэтому его можно адаптировать к любому языку программирования.
latitude = 41.145556; // (φ)
longitude = -73.995; // (λ)
mapWidth = 200;
mapHeight = 100;
// get x value
x = (longitude+180)*(mapWidth/360)
// convert from degrees to radians
latRad = latitude*PI/180;
// get y value
mercN = ln(tan((PI/4)+(latRad/2)));
y = (mapHeight/2)-(mapWidth*mercN/(2*PI));
Источники:
- Комитет по геоматике OGP, Руководство № 7, часть 2: Преобразование и преобразование координат
- Вывод проекции Меркатора
- Национальный Атлас: Картографические Проекции
- Проекция карты Меркатор
РЕДАКТИРОВАТЬ Создал рабочий пример в PHP (потому что я сосу на Java)
https://github.com/mfeldheim/mapStuff.git
EDIT2
Отличная анимация проекции Меркатора
Вы не можете просто переместить из долготы / широты в x / y таким образом, потому что мир не плоский. Вы смотрели на этот пост? Преобразование долготы / широты в координаты X/Y
ОБНОВЛЕНИЕ - 18.01.13
Я решил нанести удар, и вот как я это делаю:
public class MapService {
// CHANGE THIS: the output path of the image to be created
private static final String IMAGE_FILE_PATH = "/some/user/path/map.png";
// CHANGE THIS: image width in pixel
private static final int IMAGE_WIDTH_IN_PX = 300;
// CHANGE THIS: image height in pixel
private static final int IMAGE_HEIGHT_IN_PX = 500;
// CHANGE THIS: minimum padding in pixel
private static final int MINIMUM_IMAGE_PADDING_IN_PX = 50;
// formula for quarter PI
private final static double QUARTERPI = Math.PI / 4.0;
// some service that provides the county boundaries data in longitude and latitude
private CountyService countyService;
public void run() throws Exception {
// configuring the buffered image and graphics to draw the map
BufferedImage bufferedImage = new BufferedImage(IMAGE_WIDTH_IN_PX,
IMAGE_HEIGHT_IN_PX,
BufferedImage.TYPE_INT_RGB);
Graphics2D g = bufferedImage.createGraphics();
Map<RenderingHints.Key, Object> map = new HashMap<RenderingHints.Key, Object>();
map.put(RenderingHints.KEY_INTERPOLATION, RenderingHints.VALUE_INTERPOLATION_BICUBIC);
map.put(RenderingHints.KEY_RENDERING, RenderingHints.VALUE_RENDER_QUALITY);
map.put(RenderingHints.KEY_ANTIALIASING, RenderingHints.VALUE_ANTIALIAS_ON);
RenderingHints renderHints = new RenderingHints(map);
g.setRenderingHints(renderHints);
// min and max coordinates, used in the computation below
Point2D.Double minXY = new Point2D.Double(-1, -1);
Point2D.Double maxXY = new Point2D.Double(-1, -1);
// a list of counties where each county contains a list of coordinates that form the county boundary
Collection<Collection<Point2D.Double>> countyBoundaries = new ArrayList<Collection<Point2D.Double>>();
// for every county, convert the longitude/latitude to X/Y using Mercator projection formula
for (County county : countyService.getAllCounties()) {
Collection<Point2D.Double> lonLat = new ArrayList<Point2D.Double>();
for (CountyBoundary countyBoundary : county.getCountyBoundaries()) {
// convert to radian
double longitude = countyBoundary.getLongitude() * Math.PI / 180;
double latitude = countyBoundary.getLatitude() * Math.PI / 180;
Point2D.Double xy = new Point2D.Double();
xy.x = longitude;
xy.y = Math.log(Math.tan(QUARTERPI + 0.5 * latitude));
// The reason we need to determine the min X and Y values is because in order to draw the map,
// we need to offset the position so that there will be no negative X and Y values
minXY.x = (minXY.x == -1) ? xy.x : Math.min(minXY.x, xy.x);
minXY.y = (minXY.y == -1) ? xy.y : Math.min(minXY.y, xy.y);
lonLat.add(xy);
}
countyBoundaries.add(lonLat);
}
// readjust coordinate to ensure there are no negative values
for (Collection<Point2D.Double> points : countyBoundaries) {
for (Point2D.Double point : points) {
point.x = point.x - minXY.x;
point.y = point.y - minXY.y;
// now, we need to keep track the max X and Y values
maxXY.x = (maxXY.x == -1) ? point.x : Math.max(maxXY.x, point.x);
maxXY.y = (maxXY.y == -1) ? point.y : Math.max(maxXY.y, point.y);
}
}
int paddingBothSides = MINIMUM_IMAGE_PADDING_IN_PX * 2;
// the actual drawing space for the map on the image
int mapWidth = IMAGE_WIDTH_IN_PX - paddingBothSides;
int mapHeight = IMAGE_HEIGHT_IN_PX - paddingBothSides;
// determine the width and height ratio because we need to magnify the map to fit into the given image dimension
double mapWidthRatio = mapWidth / maxXY.x;
double mapHeightRatio = mapHeight / maxXY.y;
// using different ratios for width and height will cause the map to be stretched. So, we have to determine
// the global ratio that will perfectly fit into the given image dimension
double globalRatio = Math.min(mapWidthRatio, mapHeightRatio);
// now we need to readjust the padding to ensure the map is always drawn on the center of the given image dimension
double heightPadding = (IMAGE_HEIGHT_IN_PX - (globalRatio * maxXY.y)) / 2;
double widthPadding = (IMAGE_WIDTH_IN_PX - (globalRatio * maxXY.x)) / 2;
// for each country, draw the boundary using polygon
for (Collection<Point2D.Double> points : countyBoundaries) {
Polygon polygon = new Polygon();
for (Point2D.Double point : points) {
int adjustedX = (int) (widthPadding + (point.getX() * globalRatio));
// need to invert the Y since 0,0 starts at top left
int adjustedY = (int) (IMAGE_HEIGHT_IN_PX - heightPadding - (point.getY() * globalRatio));
polygon.addPoint(adjustedX, adjustedY);
}
g.drawPolygon(polygon);
}
// create the image file
ImageIO.write(bufferedImage, "PNG", new File(IMAGE_FILE_PATH));
}
}
РЕЗУЛЬТАТ: Ширина изображения = 600 пикселей, высота изображения = 600 пикселей, отступы изображения = 50 пикселей
РЕЗУЛЬТАТ: Ширина изображения = 300 пикселей, Высота изображения = 500 пикселей, Заполнение изображения = 50 пикселей
Java-версия исходного кода JavaScript API Карт Google v3 выглядит следующим образом, без проблем работает
public final class GoogleMapsProjection2
{
private final int TILE_SIZE = 256;
private PointF _pixelOrigin;
private double _pixelsPerLonDegree;
private double _pixelsPerLonRadian;
public GoogleMapsProjection2()
{
this._pixelOrigin = new PointF(TILE_SIZE / 2.0,TILE_SIZE / 2.0);
this._pixelsPerLonDegree = TILE_SIZE / 360.0;
this._pixelsPerLonRadian = TILE_SIZE / (2 * Math.PI);
}
double bound(double val, double valMin, double valMax)
{
double res;
res = Math.max(val, valMin);
res = Math.min(res, valMax);
return res;
}
double degreesToRadians(double deg)
{
return deg * (Math.PI / 180);
}
double radiansToDegrees(double rad)
{
return rad / (Math.PI / 180);
}
PointF fromLatLngToPoint(double lat, double lng, int zoom)
{
PointF point = new PointF(0, 0);
point.x = _pixelOrigin.x + lng * _pixelsPerLonDegree;
// Truncating to 0.9999 effectively limits latitude to 89.189. This is
// about a third of a tile past the edge of the world tile.
double siny = bound(Math.sin(degreesToRadians(lat)), -0.9999,0.9999);
point.y = _pixelOrigin.y + 0.5 * Math.log((1 + siny) / (1 - siny)) *- _pixelsPerLonRadian;
int numTiles = 1 << zoom;
point.x = point.x * numTiles;
point.y = point.y * numTiles;
return point;
}
PointF fromPointToLatLng(PointF point, int zoom)
{
int numTiles = 1 << zoom;
point.x = point.x / numTiles;
point.y = point.y / numTiles;
double lng = (point.x - _pixelOrigin.x) / _pixelsPerLonDegree;
double latRadians = (point.y - _pixelOrigin.y) / - _pixelsPerLonRadian;
double lat = radiansToDegrees(2 * Math.atan(Math.exp(latRadians)) - Math.PI / 2);
return new PointF(lat, lng);
}
public static void main(String []args)
{
GoogleMapsProjection2 gmap2 = new GoogleMapsProjection2();
PointF point1 = gmap2.fromLatLngToPoint(41.850033, -87.6500523, 15);
System.out.println(point1.x+" "+point1.y);
PointF point2 = gmap2.fromPointToLatLng(point1,15);
System.out.println(point2.x+" "+point2.y);
}
}
public final class PointF
{
public double x;
public double y;
public PointF(double x, double y)
{
this.x = x;
this.y = y;
}
}
Только JAVA?
Код Python здесь! См. Преобразование точки широты / долготы в пиксели (x,y) на проекции меркатора.
import math
from numpy import log as ln
# Define the size of map
mapWidth = 200
mapHeight = 100
def convert(latitude, longitude):
# get x value
x = (longitude + 180) * (mapWidth / 360)
# convert from degrees to radians
latRad = (latitude * math.pi) / 180
# get y value
mercN = ln(math.tan((math.pi / 4) + (latRad / 2)))
y = (mapHeight / 2) - (mapWidth * mercN / (2 * math.pi))
return x, y
print(convert(41.145556, 121.2322))
Ответ:
(167.35122222222225, 24.877939817552335)
Я хотел бы отметить, что код в границах процедуры должен читать
double bound(double val, double valMin, double valMax)
{
double res;
res = Math.max(val, valMin);
res = Math.min(res, valMax);
return res;
}
Я здесь новичок, просто чтобы написать, так как я слежу за сообществом в течение нескольких лет. Я счастлив, что могу внести свой вклад.
На поиски этого у меня ушел практически день, и ваш вопрос побудил меня продолжить поиски.
Я пришел к следующей функции, которая работает! Источники для этой статьи: https://towardsdatascience.com/geotiff-coordinate-querying-with-javascript-5e6caaaf88cf
var bbox = [minLong, minLat, maxLong, maxLat];
var pixelWidth = mapWidth;
var pixelHeight = mapHeight;
var bboxWidth = bbox[2] - bbox[0];
var bboxHeight = bbox[3] - bbox[1];
var convertToXY = function(latitude, longitude) {
var widthPct = ( longitude - bbox[0] ) / bboxWidth;
var heightPct = ( latitude - bbox[1] ) / bboxHeight;
var x = Math.floor( pixelWidth * widthPct );
var y = Math.floor( pixelHeight * ( 1 - heightPct ) );
return { x, y };
}
public static String getTileNumber(final double lat, final double lon, final int zoom) {
int xtile = (int)Math.floor( (lon + 180) / 360 * (1<<zoom) ) ;
int ytile = (int)Math.floor( (1 - Math.log(Math.tan(Math.toRadians(lat)) + 1 / Math.cos(Math.toRadians(lat))) / Math.PI) / 2 * (1<<zoom) ) ;
if (xtile < 0)
xtile=0;
if (xtile >= (1<<zoom))
xtile=((1<<zoom)-1);
if (ytile < 0)
ytile=0;
if (ytile >= (1<<zoom))
ytile=((1<<zoom)-1);
return("" + zoom + "/" + xtile + "/" + ytile);
}
}