Skip to content

Instantly share code, notes, and snippets.

@clody69
Created October 4, 2013 19:05
Show Gist options
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save clody69/6831028 to your computer and use it in GitHub Desktop.
Save clody69/6831028 to your computer and use it in GitHub Desktop.
Geoconverter from tile (z,y,zoom), to geodetic coordinates (lat, long) to elliptical mercator UTM in meters (x,y). http://www.maptiler.org/google-maps-coordinates-tile-bounds-projection/ http://wiki.openstreetmap.org/wiki/Slippy_map_tilenames http://wiki.openstreetmap.org/wiki/Mercator
module GeoConverter
R_MAJOR = 6378137.0 #Equatorial Radius, WGS84
R_MINOR = 6356752.314245179 #defined as constant
def to_rad(d)
d * (Math::PI / 180.0)
end
def to_deg(r)
r / (Math::PI / 180.0)
end
class Tile
include GeoConverter
attr_reader :x, :y, :z
def initialize(x, y, z)
@x, @y, @z = x, y, z
@nw = LatLon.new( tile2lat(y, z ), tile2long(x,z ) )
@se = LatLon.new( tile2lat(y + 1, z ), tile2long(x + 1,z ) )
end
def north
@nw.lat
end
def west
@nw.lon
end
def south
@se.lat
end
def east
@se.lon
end
def to_boundingBox
{top: north, left: west, bottom: south, right: east}
end
def to_boundingBoxUTM
puts @nw.to_utm
tl = Hash[ [:left, :top].zip(@nw.to_utm.values())]
br = Hash[[:right, :bottom].zip(@se.to_utm.values())]
tl.merge(br)
end
def tile2lat(y, z)
n = Math::PI - ( 2.0 * Math::PI * y) / (2.0 ** z)
to_deg( Math.atan( Math.sinh( n ) ) )
end
def tile2long(x, z)
x / (2.0 ** z) * 360.0 -180
end
end
class LatLon
include GeoConverter
attr_reader :lat, :lon
# Create a new coordinate instance based on latitude and longitude
def initialize(lat, lon)
raise Exception, "Invalid longitude #{lon}" unless (-180.0...180.0).member? lon
@lat, @lon = lat, lon
end
# Textual representation of the coordinate
def to_s
north_south = if @lat >= 0.0 then 'N' else 'S' end
east_west = if @lon >= 0.0 then 'E' else 'W' end
'%0.6f%s %0.6f%s' % [@lat.abs, north_south, @lon.abs, east_west]
end
# Convert the coordinate in latutude/longitude into the UTM coordinate system
def to_utm()
x = R_MAJOR * to_rad(@lon)
@lat = 89.5 if (@lat > 89.5)
@lat = -89.5 if (@lat < -89.5)
temp = R_MINOR / R_MAJOR
es = 1.0 - (temp * temp)
eccent = Math.sqrt(es)
phi = to_rad(@lat)
sinphi = Math.sin(phi)
con = eccent * sinphi
com = 0.5 * eccent
con2 = ((1.0-con)/(1.0+con)) ** com
ts = Math.tan(0.5 * (Math::PI*0.5 - phi))/con2
y = 0 - R_MAJOR * Math.log(ts);
{x: x, y: y}
end
end
end
puts GeoConverter::LatLon.new(60.15244221438077,24.9169921875).to_utm
puts GeoConverter::Tile.new(9326, 4744, 14).to_boundingBox
puts GeoConverter::Tile.new(9326, 4744, 14).to_boundingBoxUTM
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment