module NSWTopo
  module Contour
    using Helpers
    include VectorRender, DEM, Log

    CREATE = %w[interval index auxiliary smooth simplify thin density min-length no-depression knolls fill epsg]
    DEFAULTS = YAML.load <<~YAML
      interval: 5
      smooth: 0.2
      density: 4.0
      min-length: 2.0
      knolls: 0.2
      section: 100
      stroke: hsl(40,100%,25%)
      stroke-width: 0.08
      Auxiliary:
        stroke-dasharray: 0.5 0.5
        stroke-dashoffset: 0.5
      Depression:
        symbolise:
          interval: 2.0
          line:
            stroke-width: 0.12
            y2: -0.3
      labels:
        font-size: 1.4
        letter-spacing: 0.05
        orientation: downhill
        collate: true
        min-radius: 5
        max-turn: 20
        sample: 10
        minimum-area: 70
        separation:
          self: 40
          other: 15
          along: 100
    YAML

    def margin
      { mm: [3 * @smooth, 1].max }
    end

    def check_geos!
      json = OS.ogr2ogr "-dialect", "SQLite", "-sql", "SELECT geos_version() AS version", "-f", "GeoJSON", "-lco", "RFC7946=NO", "/vsistdout/", "GeoJSON:/vsistdin?buffer_limit=-1" do |stdin|
        stdin.write GeoJSON::Collection.new.to_json
      end
      raise unless version = JSON.parse(json).dig("features", 0, "properties", "version")
      raise unless (version.split(?-).first.split(?.).map(&:to_i) <=> [3, 3]) >= 0
    rescue OS::Error, JSON::ParserError, RuntimeError
      raise "contour thinning requires GDAL with SpatiaLite and GEOS support"
    end

    def get_features
      @simplify ||= [@map.to_mm(0.5 * @interval) / Math::tan(Math::PI * 85 / 180), 0.05].min
      @index ||= 10 * @interval
      @params = {
        "Index" => { "stroke-width" => 2 * @params["stroke-width"] },
        "labels" => { "fill" => @fill || "black" }
      }.deep_merge(@params)

      check_geos! if @thin
      raise "%im index interval not a multiple of %im contour interval" % [@index, @interval] unless @index % @interval == 0

      Dir.mktmppath do |temp_dir|
        dem_path, blur_path = temp_dir / "dem.tif", temp_dir / "dem.blurred.tif"

        if @smooth.zero?
          get_dem temp_dir, blur_path
        else
          get_dem temp_dir, dem_path
          blur_dem dem_path, blur_path
        end

        db_flags = @thin ? %w[-f SQLite -dsco SPATIALITE=YES] : ["-f", "ESRI Shapefile"]
        db_path = temp_dir / "contour"

        log_update "%s: generating contour lines" % @name
        json = OS.gdal_contour "-q", "-a", "elevation", "-i", @interval, "-f", "GeoJSON", "-lco", "RFC7946=NO", blur_path, "/vsistdout/"
        contours = GeoJSON::Collection.load(json, projection: @map.projection).map! do |feature|
          id, elevation = feature.values_at "ID", "elevation"
          properties = { "id" => id, "elevation" => elevation, "modulo" => elevation % @index, "depression" => feature.closed? && feature.anticlockwise? ? 1 : 0}
          feature.with_properties(properties)
        end
        raise "no elevation data found in map area" unless contours.any?

        if @no_depression.nil?
          candidates = contours.select do |feature|
            feature["depression"] == 1
          end
          index = RTree.load(candidates, &:bounds)

          contours.reject! do |feature|
            next unless feature["depression"] == 1
            index.search(feature.bounds).none? do |other|
              next if other == feature
              feature.to_polygon.contains?(other.first) || other.to_polygon.contains?(feature.first)
            end
          end
        end

        contours.reject! do |feature|
          feature.closed? &&
          feature.bounds.all? { |min, max| max - min < @knolls }
        end.reject! do |feature|
          feature["elevation"].zero?
        end

        contours.each_slice(100).inject(nil) do |update, features|
          OS.ogr2ogr "-a_srs", @map.projection, "-nln", "contour", *update, "-simplify", @simplify, *db_flags, db_path, "GeoJSON:/vsistdin?buffer_limit=-1" do |stdin|
            stdin.write contours.with_features(features).to_json
          end
          %w[-update -append]
        end

        if @thin
          slope_tif_path = temp_dir / "slope.tif"
          slope_vrt_path = temp_dir / "slope.vrt"

          log_update "%s: generating slope masks" % @name
          OS.gdaldem "slope", blur_path, slope_tif_path, "-compute_edges"
          json = OS.gdalinfo "-json", slope_tif_path
          width, height = JSON.parse(json)["size"]
          srcwin = [ -2, -2, width + 4, height + 4 ]
          OS.gdal_translate "-srcwin", *srcwin, "-a_nodata", "none", "-of", "VRT", slope_tif_path, slope_vrt_path

          multiplier = @index / @interval
          Enumerator.new do |yielder|
            keep = 0...multiplier
            until keep.one?
              keep, drop = keep.count.even? ? keep.each_slice(2).entries.transpose : [[0], keep.drop(1)]
              yielder << drop
            end
          end.inject(multiplier) do |count, drop|
            angle = Math::atan(@index * @density / count) * 180.0 / Math::PI
            mask_path = temp_dir / "mask.#{count}.sqlite"

            OS.gdal_contour "-nln", "ring", "-a", "angle", "-fl", angle, *db_flags, slope_vrt_path, mask_path

            OS.ogr2ogr "-update", "-nln", "mask", "-nlt", "MULTIPOLYGON", mask_path, mask_path, "-dialect", "SQLite", "-sql", <<~SQL
              SELECT
                ST_Buffer(ST_Buffer(ST_Polygonize(geometry), #{0.5 * @min_length}, 6), #{-0.5 * @min_length}, 6) AS geometry
              FROM ring
            SQL

            drop.each do |index|
              OS.ogr2ogr "-nln", "mask", "-update", "-append", "-explodecollections", "-q", db_path, mask_path, "-dialect", "SQLite", "-sql", <<~SQL
                SELECT geometry, #{index * @interval} AS modulo
                FROM mask
              SQL
            end

            count - drop.count
          end

          log_update "%s: thinning contour lines" % @name
          OS.ogr2ogr "-nln", "divided", "-update", "-explodecollections", db_path, db_path, "-dialect", "SQLite", "-sql", <<~SQL
            WITH intersecting(contour, mask) AS (
              SELECT contour.rowid, mask.rowid
              FROM contour
              INNER JOIN mask
              ON
                mask.modulo = contour.modulo AND
                contour.rowid IN (
                  SELECT rowid FROM SpatialIndex
                  WHERE
                    f_table_name = 'contour' AND
                    search_frame = mask.geometry
                ) AND
                ST_Relate(contour.geometry, mask.geometry, 'T********')
            )

            SELECT contour.geometry, contour.id, contour.elevation, contour.modulo, contour.depression, 1 AS unmasked, 1 AS unaltered
            FROM contour
            LEFT JOIN intersecting ON intersecting.contour = contour.rowid
            WHERE intersecting.contour IS NULL

            UNION SELECT ExtractMultiLinestring(ST_Difference(contour.geometry, ST_Collect(mask.geometry))) AS geometry, contour.id, contour.elevation, contour.modulo, contour.depression, 1 AS unmasked, 0 AS unaltered
            FROM contour
            INNER JOIN intersecting ON intersecting.contour = contour.rowid
            INNER JOIN mask ON intersecting.mask = mask.rowid
            GROUP BY contour.rowid
            HAVING min(ST_Relate(contour.geometry, mask.geometry, '**T******'))

            UNION SELECT ExtractMultiLinestring(ST_Intersection(contour.geometry, ST_Collect(mask.geometry))) AS geometry, contour.id, contour.elevation, contour.modulo, contour.depression, 0 AS unmasked, 0 AS unaltered
            FROM contour
            INNER JOIN intersecting ON intersecting.contour = contour.rowid
            INNER JOIN mask ON intersecting.mask = mask.rowid
            GROUP BY contour.rowid
          SQL

          OS.ogr2ogr "-nln", "thinned", "-update", "-explodecollections", db_path, db_path, "-dialect", "SQLite", "-sql", <<~SQL
            SELECT ST_LineMerge(ST_Collect(geometry)) AS geometry, id, elevation, modulo, depression, unaltered
            FROM divided
            WHERE unmasked OR ST_Length(geometry) < #{@min_length}
            GROUP BY id, elevation, modulo, unaltered
          SQL

          OS.ogr2ogr "-nln", "contour", "-update", "-overwrite", db_path, db_path, "-dialect", "SQLite", "-sql", <<~SQL
            SELECT geometry, id, elevation, modulo, depression
            FROM thinned
            WHERE unaltered OR ST_Length(geometry) > #{@min_length}
          SQL
        end

        json = OS.ogr2ogr "-f", "GeoJSON", "-lco", "RFC7946=NO", "/vsistdout/", db_path, "contour"
        GeoJSON::Collection.load(json, projection: @map.projection).map! do |feature|
          elevation, modulo, depression = feature.values_at "elevation", "modulo", "depression"
          category = case
          when @auxiliary && elevation % (2 * @interval) != 0 then %w[Auxiliary]
          when modulo.zero? then %w[Index]
          else %w[Standard]
          end
          category << "Depression" if depression == 1

          properties = Hash[]
          properties["elevation"] = elevation
          properties["category"] = category
          properties["label"] = elevation.to_i.to_s if modulo.zero?
          feature.with_properties(properties)
        end
      end
    end

    def to_s
      elevations = features.map do |feature|
        [feature["elevation"], feature["category"].include?("Index")]
      end.uniq.sort_by(&:first)
      range = elevations.map(&:first).minmax
      interval, index = %i[itself last].map do |selector|
        elevations.select(&selector).map(&:first).each_cons(2).map { |e0, e1| e1 - e0 }.min
      end
      [["%im intervals", interval], ["%im indices", index], ["%im-%im elevation", (range if range.all?)]].select(&:last).map do |label, value|
        label % value
      end.join(", ").prepend("%s: " % @name)
    end
  end
end
