find_gsk3_bounding_boxes.rb

# sudo aptitude install libproj-dev 
# gem install proj4rb

# This Ruby script :
#   finds all the citygml files
#   reads them
#   find the 2D convex hull
#   writes it to a kml file

require 'fileutils'
require 'proj4'
require 'builder'

require_relative 'convex_hull'

Resolution = 0.1
RepositoryDir = "../TestRepository/"

class Point
  attr_reader :x,:y,:z
  @@dest_prj  = Proj4::Projection.new("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs")
  @@src_prj = Proj4::Projection.new("+proj=tmerc +lat_0=0 +lon_0=9 +k=1 +x_0=3500000 +y_0=0 +ellps=bessel +datum=potsdam +units=m +no_defs")
  def initialize(l)
    @x,@y,@z = l.map{|v| v.to_f}
  end


  def to_s
    [x,y,z].join(',')
  end

  def to_wgs84
    wgs84_point = @@src_prj.transform(@@dest_prj, proj4_point)
    ['%.7f' % (wgs84_point.x*180/Math::PI), '%.7f' % (wgs84_point.y*180/Math::PI),0].join(',')
  end

  def <=>(p2)
    [x,y]<=>[p2.x,p2.y]
  end

  private

  def proj4_point
    Proj4::Point.new(x,y)
  end
end



xml = Builder::XmlMarkup.new(:target=>STDOUT, :indent=>2)
xml.instruct! :xml, :version=>"1.0", :encoding=>"UTF-8"

xml.kml do
  xml.Document do
    FileUtils.cd(RepositoryDir){
      xmls_and_gmls = Dir["**/*"].select{|f|
        File.extname(f).downcase=~/(g|x)ml/
      }.sort_by{|f| File.basename(f)}

      citygmls = xmls_and_gmls.select{|xml|
        File.open(xml,'r'){|f|
          f.read(1000)
        }=~/citymodel/i
      }

      citygmls.each{|citygml|
        basename = File.basename(citygml)
        project = File.dirname(citygml).split('/').last.gsub('.simstadt','')

        content = File.read(citygml).force_encoding("ISO-8859-1").encode("utf-8", replace: nil)
        coordinates = content.gsub(/(low|upp)erCorner.*?Corner/,'').scan(/(?<![\d\.])(3\d\d\d\d\d\d[\.\d]*) (5\d\d\d\d\d\d[\.\d]*) ([\d\.]+)/)
        next if coordinates.empty?
        xml.Placemark do
          xml.name basename
          xml.description project+">"+basename
          xml.Polygon do
            xml.outerBoundaryIs do
              xml.LinearRing do
                xml.tessellate 1
                xyz_points = coordinates.map{|l| Point.new(l)}
                xy_points = xyz_points.uniq{|p| [(p.x/Resolution).round,(p.y/Resolution).round]}
                hull = ConvexHull.calculate(xy_points).map{|p| p.to_wgs84}
                xml.coordinates hull.join(' ')
              end
            end
          end
        end
      }
    }
  end
end