国土地理院のPNG標高タイルからGeoTIFF生成

はじめてのCrystal。といってもほとんどRuby。

Crystal本体の他に必要なもの
GitHub - stumpycr/stumpy_png: Read/Write PNG images in pure Crystal
GDAL
stumpy_pngのインストールはリンク先参照。
GDALはaptで。

$ sudo apt install gdal-bin

crystalでgdal用の入力データを作成。(Esri ASCII ラスター形式)
gdal_translateで処理。

素材は地理院タイルから

require "http/client"
require "stumpy_png"
include Math

def elevation(canvas)
    tile = [] of Array(Float64)
    0.upto(canvas.height - 1) do |h|
        row = [] of Float64
        0.upto(canvas.width - 1) do |w|
            r, g, b = canvas[w, h].to_rgb8  # r,g,bはUInt8
            v = 65536*r.to_f + 256*g.to_f + b.to_f
            if v < 2**23
                v = v*0.01
            elsif v > 2**23
                v = (v-2**24)*0.01
            else # v == 2**23
                v = -9999.0
            end
            row.push v
        end
        tile.push row
    end

    return tile
end

def header(z, x, y)
    # タイル座標 -> WEBメルカトル
    r = 6378137.0
    z = z - 1
    xm = PI*r*(x.to_f/(2**z) - 1.0)
    ym = PI*r*(1.0 - (y+1).to_f/(2**z))   # (y+1)で左下隅
    dt = PI*r/(2**z)/256

    return <<-HERE
    ncols\t#{256}
    nrows\t#{256}
    xllcorner\t#{xm}
    yllcorner\t#{ym}
    cellsize\t#{dt}
    NODATA_value\t-9999
    HERE
end


z = 14      # ズームレベル
x = 14505   # タイルX座標  
y = 6469    # タイルY座標

response = HTTP::Client.get "https://cyberjapandata.gsi.go.jp/xyz/dem_png/#{z}/#{x}/#{y}.png"
File.open("dem-test.png", "w") do |f|
    f.print(response.body)
end

canvas = StumpyPNG.read("dem-test.png")
tile = elevation(canvas)

File.open("dem-grid.txt", "w") do |f|
    f.puts header(z, x, y)
    tile.each do |h|
        f.puts h.join(" ")
    end
end

$ gdal_translate -of GTiff -a_nodata -9999 -a_srs EPSG:3857 dem-grid.txt dem-test.tif

(表示のためPNGに変換)

変換にはSAGA GISを使った。

$ sudo apt install saga

saga_cmd io_grid_image 0 -COL_PALETTE 2 -GRID input.tif -FILE output.png