📗Develop a DEM1A parser and PMTiles

[hfu]

Current progress

Current code

require 'nokogiri'
require 'gdal'
require 'zip'

def convert(input, tif_path)
  puts "#{tif_path}"
  doc = Nokogiri::XML(input) {|config| config.huge}
  lower_corner = doc.at_xpath('//gml:lowerCorner').text.split.map(&:to_f)
  upper_corner = doc.at_xpath('//gml:upperCorner').text.split.map(&:to_f)
  high = doc.at_xpath('//gml:GridEnvelope/gml:high').text.split.map(&:to_i)
  start_point = doc.at_xpath('//gml:GridFunction/gml:startPoint').text.split.map(&:to_i)
  width = high[0] + 1
  height = high[1] + 1

  #tuple_list = doc.at_xpath('//gml:tupleList').text.strip.split("\n").map do |line|
  #  type, value = line.split(',')
  #  value.to_f
  #end

  tuple_list = []
  (start_point[1] * width + start_point[0]).times {|i|
    tuple_list << -9999.0
  }
  input.each_line do |line|
    r = line.split(',')
    tuple_list << r[1].to_f if r.size == 2
  end
  (width * height - tuple_list.size).times {|i|
    tuple_list << -9999.0
  }

  #data = Array.new(width * height, -9999)
  #tuple_list.each_with_index do |value, index|
  #  data[index] = value
  #end

  puts "Start Point: #{start_point}"
  puts "Lower Corner: #{lower_corner}"
  puts "Upper Corner: #{upper_corner}"
  puts "Width: #{width}, Height: #{height}"
  puts "Tuple List size: #{tuple_list.size}, pixels: #{width * height}"
  puts

  driver = Gdal::Gdal.get_driver_by_name('GTiff')
  File.delete(tif_path) if File.exist?(tif_path)
  dataset = driver.create(tif_path, width, height, 1, Gdal::Gdalconst::GDT_FLOAT32)

  min_y, min_x = lower_corner
  max_y, max_x = upper_corner
  pixel_width = (max_x - min_x) / width
  pixel_height = (max_y - min_y) / height
  geo_transform = [min_x, pixel_width, 0, max_y, 0, -pixel_height]
  dataset.set_geo_transform(geo_transform)

  srs = Gdal::Osr::SpatialReference.new
  srs.import_from_epsg(6668)
  dataset.set_projection(srs.export_to_wkt)

  band = dataset.get_raster_band(1)

  flattened_data = tuple_list.pack('f*')
  band.write_raster(0, 0, width, height, flattened_data)

  band.set_no_data_value(-9999.0)

  band.flush_cache
  dataset = nil
end

zip_path = "FG-GML-5741-53-DEM1A.zip"

Zip::File.open(zip_path) do |zip_file|
  zip_file.each do |entry|
    src_name = entry.name
    next unless src_name.downcase.end_with?('.xml')
    dst_name = src_name.sub('.xml', '.tif')
    input = entry.get_input_stream.read
    convert(input, dst_name)
  end
end

[hfu]

New version under development

https://github.com/UNopenGIS/gmldem2tif/blob/main/gmldem2tif.rb

[hfu]

File downloading

• [x] 5740
• [x] 5741
• [x] 5841
• [x] 5941
• [x] 5942
• [x] 6041
• [x] 5540
• [x] 5541
• [x] 5640
• [x] 5641

[hfu]

DEM1A Terrain Tiles in PMTiles are getting ready. To be published soon!

[hfu]

A version in 2024-05-27 morning

• a.tif 15GB
• mbtiles 687MB
• pmtiles 666MB

Updating a.tif by:

gdal_merge.py -co COMPRESS=LZW -n -9999.0 -o a.tif dst/*.tif

[hfu]

gdal_merge.py -co COMPRESS=LZW -co BIGTIFF=YES -n -9999.0 -o a.tif dst/*.tif 

[hfu]

a.tif created

[hfu]

I am happy to relase 1m GSD terrain tiles sourced from Geospatial Information Authority of Japan (GSI) in Dhaka, because we received the use approval from GSI recently. https://observablehq.com/d/c0a3f48cf7111cd2

[hfu]

I am closing this.

I thank @smellman for https://github.com/UNopenGIS/gmldem2tif/pull/1. It would be a new chapter of this story!