module NSWTopo::Contour

Constants

CREATE
DEFAULTS

Public Instance Methods

check_geos!() click to toggle source
# File lib/nswtopo/layer/contour.rb, line 38
def check_geos!
  json = OS.ogr2ogr "-dialect", "SQLite", "-sql", "SELECT geos_version() AS version", "-f", "GeoJSON", "/vsistdout/", "/vsistdin/" 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
get_features() click to toggle source
# File lib/nswtopo/layer/contour.rb, line 48
    def get_features
      @simplify ||= [0.5 * @interval / Math::tan(Math::PI * 85 / 180), 0.001 * 0.05 * @map.scale].min
      @index ||= 10 * @interval
      @params = {
        "Index" => { "stroke-width" => 2 * @params["stroke-width"] },
        "labels" => { "fill" => @fill || @params["stroke"] }
      }.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, @map.projection

        if @no_depression.nil?
          candidates = contours.select do |feature|
            feature.coordinates.last == feature.coordinates.first &&
            feature.coordinates.anticlockwise?
          end.each 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.coordinates.first.within?(other.coordinates) ||
              other.coordinates.first.within?(feature.coordinates)
            end
          end
        end

        contours.reject! do |feature|
          feature.coordinates.last == feature.coordinates.first &&
          feature.bounds.all? { |min, max| max - min < @knolls * @map.scale / 1000.0 }
        end

        contours.each do |feature|
          id, elevation, depression = feature.values_at "ID", "elevation", "depression"
          feature.properties.replace("id" => id, "elevation" => elevation, "modulo" => elevation % @index, "depression" => depression || 0)
        end

        contours.reject! do |feature|
          feature["elevation"].zero?
        end

        OS.ogr2ogr "-a_srs", @map.projection, "-nln", "contour", "-simplify", @simplify, *db_flags, db_path, "GeoJSON:/vsistdin/" do |stdin|
          stdin.write contours.to_json
        end

        if @thin
          slope_tif_path = temp_dir / "slope.tif"
          slope_vrt_path = temp_dir / "slope.vrt"
          min_length = @min_length * @map.scale / 1000.0

          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
          case multiplier
          when  4 then [ [1,3], 2 ]
          when  5 then [ [1,4], [2,3] ]
          when  6 then [ [1,4], [2,5], 3 ]
          when  7 then [ [2,5], [1,3,6], 4 ]
          when  8 then [ [1,3,5,7], [2,6], 4 ]
          when  9 then [ [1,4,7], [2,5,8], [3,6] ]
          when 10 then [ [2,5,8], [1,4,6,9], [3,7] ]
          else raise "contour thinning not available for specified index interval"
          end.inject(multiplier) do |count, (*drop)|
            angle = Math::atan(1000.0 * @index * @density / @map.scale / 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, @map.projection).each do |feature|
          elevation, modulo, depression = feature.values_at "elevation", "modulo", "depression"
          category = modulo.zero? ? %w[Index] : %w[Standard]
          category << "Depression" if depression == 1
          feature.clear
          feature["elevation"] = elevation
          feature["category"] = category
          feature["label"] = elevation.to_i.to_s if modulo.zero?
        end
      end
    end
margin() click to toggle source
# File lib/nswtopo/layer/contour.rb, line 34
def margin
  { mm: [3 * @smooth, 1].min }
end
to_s() click to toggle source
# File lib/nswtopo/layer/contour.rb, line 220
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