2

I am new to GIS but I am working on a project mapping a building using floor plans and JOSM.

Using the PicLayer plugin I added the .png floor plan as a layer.

I wrote a .cal file to set the scale, rotation, and center position of the image.

In the JOSM default Mercator projection, the image moved only slightly when the center was changed (using lat and lon decimal degrees). The image was still centered somewhere around zero. In WGS84 the image moved to the correct location but was sheared and stretched.

I found this reference on converting lat/lon to a Mercator projection: http://wiki.openstreetmap.org/wiki/Mercator

The Python code under Elliptical Mercator correctly converts my longitude, but the latitude is off by about 0.5% (a few metres).

import math

def merc_x(lon):
  r_major=6378137.000
  return r_major*math.radians(lon)

def merc_y(lat):
  if lat>89.5:lat=89.5
  if lat<-89.5:lat=-89.5
  r_major=6378137.000
  r_minor=6356752.3142
  temp=r_minor/r_major
  eccent=math.sqrt(1-temp**2)
  phi=math.radians(lat)
  sinphi=math.sin(phi)
  con=eccent*sinphi
  com=eccent/2
  con=((1.0-con)/(1.0+con))**com
  ts=math.tan((math.pi/2-phi)/2)/con
  y=0-r_major*math.log(ts)
  return y

My nodes and ways were lined up correctly with the image with this latitude and longitude in WGS84, so I think they are correct and the conversion is not. I am in the northern hemisphere and using an area near 49 N, -122 E.

  • So I was looking at the elliptical Mercator code because the Python for spherical Mercator didn't show X and gave me the wrong Y value. In the ActionScript section (wiki.openstreetmap.org/wiki/…), it says that the X coordinate is the longitude multiplied by a chosen scale. I guessed that I could take the output X value from the elliptical Mercator conversion and divide it by the actual longitude to get the scale, and it worked. I then took this value and multiplied it by the result of lat2y to get the Mercator Y value. Will post as answer. – Josh Vazquez Jul 28 '15 at 6:12
1

Take a look at Proj4 https://github.com/OSGeo/proj.4/wiki python's implementation pyproj https://github.com/jswhit/pyproj

You could specify src and dest projections (like EPSG:4326 and EPSG:3857) and convert coordinates between them. And pyproj will made all raw calculations.

1

Longitude conversion from elliptical Mercator (same formula in both systems?:

import math
def merc_x(lon):
  r_major=6378137.000
  return r_major*math.radians(lon)

Mercator X now known. Divide this by the longitude to get the scale value.

Latitude conversion from spherical Mercator:

import math
def lat2y(a):
  return 180.0/math.pi*math.log(math.tan(math.pi/4.0+a*(math.pi/180.0)/2.0))

Multiply the result by the scale value to get Mercator Y value.

In one function:

import math
def merc(lat, lon):
    r_major = 6378137.000
    x = r_major * math.radians(lon)
    scale = x/lon
    y = 180.0/math.pi * math.log(math.tan(math.pi/4.0 + lat * (math.pi/180.0)/2.0)) * scale
    return (x, y)
1

Note: this solution targets array inputs. Inspired by the accepted @JoshVazquez's solution.

If lat and lon are both arrays (or dataframe columns)...

The current accepted answer would work very well when lat and lon are scalar parameters.

I've however recently bumped into a scenario where lat and lon are numpy arrays (or pandas dataframe columns). To make this work for arrays / dataframe columns, I've essentially modified @JoshVazquez's solution above - essentially replacing all the math functions with numpy functions.

import numpy as np
def merc_from_arrays(lats, lons):
    r_major = 6378137.000
    x = r_major * np.radians(lons)
    scale = x/lons
    y = 180.0/np.pi * np.log(np.tan(np.pi/4.0 + lats * (np.pi/180.0)/2.0)) * scale
    return (x, y)

lats = np.asarray([54.984105, 56])
lons = np.asarray([-3.193693, -2.2])

xs, ys = merc_from_arrays(lats, lons)

print('xs: {}'.format(xs))  # => [-355520.27851004 -244902.8797452 ]
print('ys: {}'.format(ys))  # => [ 7358781.82857011  7558415.65608178]
  • 1
    There is a bug in this code: in the mere_from_arrays function, you reference variables lon and lat, but those do not exist.Presumably you meant lats and lons, the arguments to the function. – ibrewster Feb 7 at 0:20
  • thank you @ibrewster - typo corrected now. – Atlas7 Feb 9 at 10:35

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.