経緯度からWeb Mercator座標を計算するコード
経緯度からWeb Mercator座標を計算するコードを 経緯度とGoogle Mercatorとの変換式を確認してみる(OpenLayers & Proj4js)。 - Relevant, Timely, and Accurate から取り出して、タイル矩形を EPSG:3857 の GeoTIFF に吐き出せるようにするコードを完成しました。
経緯度からWeb Mercator座標を計算するコード
経緯度とGoogle Mercatorとの変換式を確認してみる(OpenLayers & Proj4js)。 - Relevant, Timely, and Accurate の JavaScript コードを Ruby コードに翻訳しただけです。
def lng2mx(lng) lng * 20037508.34 / 180; end def lat2my(lat) Math::log(Math::tan((90 + lat) * Math::PI / 360)) / (Math::PI / 180) * 20037508.34 / 180 end
mx は Mercator's x、my は Mercator's y のつもりです。いずれにせよ gdal_translate に読ませるくらいの用途しかなさそうで、ブラウザ上では使わない値なので、命名もてきとうです。
XYZスキームの1タイルを1画素にしたGeoTIFF画像を作るgdal_translateコマンドを作成するコード
こんな感じになると思います。1ピクセルずれてたりしたらすみません。多分大丈夫かとは思いますが...。
cmd = "gdal_translate -of GTiff -a_srs EPSG:3857 -a_ullr #{lng2mx(x2lng(table.min(:x), n))} #{lat2my(y2lat(table.min(:y), n))} #{lng2mx(x2lng(table.max(:x) + 1, n))} #{lat2my(y2lat(table.max(:y) + 1, n))} #{z}#{type}.png #{z}#{type}.tif"
コード全体についても、共有します。かなり無色透明なコードのはず。
#image_converter.rb require 'rubygems' require 'sequel' require 'RMagick' ### z = 15 type = 'std' ### DB = Sequel.sqlite("#{z}#{type}.sqlite") table = DB["check#{z}#{type}".to_sym] class Range def size self.max - self.min end end def x2lng(x, n) 360.0 * x / n - 180.0 end def lng2mx(lng) lng * 20037508.34 / 180; end def y2lat(y, n) rad = Math::atan(Math::sinh(Math::PI * (1 - 2.0 * y / n))) rad * 360.0 / (2 * Math::PI) end def lat2my(lat) Math::log(Math::tan((90 + lat) * Math::PI / 360)) / (Math::PI / 180) * 20037508.34 / 180 end n = 2 ** z xr = (table.min(:x)..table.max(:x)) yr = (table.min(:y)..table.max(:y)) img = Magick::Image.new(xr.size + 1, yr.size + 1) { self.background_color = 'gray' } lut = {200 => 'green', 404 => 'orange'} table.each {|r| color = lut[r[:code]] || 'red' img.store_pixels(r[:x] - xr.min, r[:y] - yr.min, 1, 1, [Magick::Pixel::from_color(color)]) } img.write("#{z}#{type}.png") cmd = "gdal_translate -of GTiff -a_srs EPSG:3857 -a_ullr #{lng2mx(x2lng(table.min(:x), n))} #{lat2my(y2lat(table.min(:y), n))} #{lng2mx(x2lng(table.max(:x) + 1, n))} #{lat2my(y2lat(table.max(:y) + 1, n))} #{z}#{type}.png #{z}#{type}.tif" p cmd system cmd