Module: NSWTopo::Contour

Includes:
DEM, Log, Vector
Defined in:
lib/nswtopo/layer/contour.rb

Constant Summary collapse

CREATE = 
%w[interval index auxiliary smooth simplify thin density min-length no-depression knolls fill]
DEFAULTS = 
YAML.load "interval: 5\nsmooth: 0.2\ndensity: 4.0\nmin-length: 2.0\nknolls: 0.2\nsection: 100\nstroke: hsl(40,100%,25%)\nstroke-width: 0.08\nAuxiliary:\n  stroke-dasharray: 0.5 0.5\n  stroke-dashoffset: 0.5\nDepression:\n  symbolise:\n    interval: 2.0\n    line:\n      stroke-width: 0.12\n      y2: -0.3\nlabels:\n  font-size: 1.4\n  letter-spacing: 0.05\n  orientation: downhill\n  collate: true\n  min-radius: 5\n  max-turn: 20\n  sample: 10\n  minimum-area: 70\n  separation:\n    self: 40\n    other: 15\n    along: 100\n"

Constants included from Vector

Vector::FONT_SCALED_ATTRIBUTES, Vector::MARGIN, Vector::SVG_ATTRIBUTES

Constants included from Log

Log::FAILURE, Log::NEUTRAL, Log::SUCCESS, Log::UPDATE

Instance Method Summary collapse

Methods included from Vector

#categorise, #create, #drawing_features, #features, #filename, #labeling_features, #params_for, #render, #svg_path_data

Methods included from DEM

#blur_dem, #get_dem

Methods included from GDALGlob

#gdal_raster?, #gdal_rasters

Methods included from Log

#log_abort, #log_neutral, #log_success, #log_update, #log_warn

Instance Method Details

#check_geos! ⇒ Object 



42
43
44
45
46
47
48
49
50
 
# File 'lib/nswtopo/layer/contour.rb', line 42

def check_geos!
  json = OS.ogr2ogr "-dialect", "SQLite", "-sql", "SELECT geos_version() AS version", "-f", "GeoJSON", "-lco", "RFC7946=NO", "/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 ⇒ Object 



52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
 
# File 'lib/nswtopo/layer/contour.rb', line 52

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 || @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, projection: @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 }
    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

    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/" do |stdin|
        stdin.write GeoJSON::Collection.new(projection: @map.projection, 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", "          SELECT\n            ST_Buffer(ST_Buffer(ST_Polygonize(geometry), \#{0.5 * @min_length}, 6), \#{-0.5 * @min_length}, 6) AS geometry\n          FROM ring\n        SQL\n\n        drop.each do |index|\n          OS.ogr2ogr \"-nln\", \"mask\", \"-update\", \"-append\", \"-explodecollections\", \"-q\", db_path, mask_path, \"-dialect\", \"SQLite\", \"-sql\", <<~SQL\n            SELECT geometry, \#{index * @interval} AS modulo\n            FROM mask\n          SQL\n        end\n\n        count - drop.count\n      end\n\n      log_update \"%s: thinning contour lines\" % @name\n      OS.ogr2ogr \"-nln\", \"divided\", \"-update\", \"-explodecollections\", db_path, db_path, \"-dialect\", \"SQLite\", \"-sql\", <<~SQL\n        WITH intersecting(contour, mask) AS (\n          SELECT contour.rowid, mask.rowid\n          FROM contour\n          INNER JOIN mask\n          ON\n            mask.modulo = contour.modulo AND\n            contour.rowid IN (\n              SELECT rowid FROM SpatialIndex\n              WHERE\n                f_table_name = 'contour' AND\n                search_frame = mask.geometry\n            ) AND\n            ST_Relate(contour.geometry, mask.geometry, 'T********')\n        )\n\n        SELECT contour.geometry, contour.id, contour.elevation, contour.modulo, contour.depression, 1 AS unmasked, 1 AS unaltered\n        FROM contour\n        LEFT JOIN intersecting ON intersecting.contour = contour.rowid\n        WHERE intersecting.contour IS NULL\n\n        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\n        FROM contour\n        INNER JOIN intersecting ON intersecting.contour = contour.rowid\n        INNER JOIN mask ON intersecting.mask = mask.rowid\n        GROUP BY contour.rowid\n        HAVING min(ST_Relate(contour.geometry, mask.geometry, '**T******'))\n\n        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\n        FROM contour\n        INNER JOIN intersecting ON intersecting.contour = contour.rowid\n        INNER JOIN mask ON intersecting.mask = mask.rowid\n        GROUP BY contour.rowid\n      SQL\n\n      OS.ogr2ogr \"-nln\", \"thinned\", \"-update\", \"-explodecollections\", db_path, db_path, \"-dialect\", \"SQLite\", \"-sql\", <<~SQL\n        SELECT ST_LineMerge(ST_Collect(geometry)) AS geometry, id, elevation, modulo, depression, unaltered\n        FROM divided\n        WHERE unmasked OR ST_Length(geometry) < \#{@min_length}\n        GROUP BY id, elevation, modulo, unaltered\n      SQL\n\n      OS.ogr2ogr \"-nln\", \"contour\", \"-update\", \"-overwrite\", db_path, db_path, \"-dialect\", \"SQLite\", \"-sql\", <<~SQL\n        SELECT geometry, id, elevation, modulo, depression\n        FROM thinned\n        WHERE unaltered OR ST_Length(geometry) > \#{@min_length}\n      SQL\n    end\n\n    json = OS.ogr2ogr \"-f\", \"GeoJSON\", \"-lco\", \"RFC7946=NO\", \"/vsistdout/\", db_path, \"contour\"\n    GeoJSON::Collection.load(json, projection: @map.projection).each do |feature|\n      elevation, modulo, depression = feature.values_at \"elevation\", \"modulo\", \"depression\"\n      category = case\n      when @auxiliary && elevation % (2 * @interval) != 0 then %w[Auxiliary]\n      when modulo.zero? then %w[Index]\n      else %w[Standard]\n      end\n      category << \"Depression\" if depression == 1\n      feature.clear\n      feature[\"elevation\"] = elevation\n      feature[\"category\"] = category\n      feature[\"label\"] = elevation.to_i.to_s if modulo.zero?\n    end\n  end\nend\n"

#margin ⇒ Object 



38
39
40
 
# File 'lib/nswtopo/layer/contour.rb', line 38

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

#to_s ⇒ Object 



227
228
229
230
231
232
233
234
235
236
237
238
 
# File 'lib/nswtopo/layer/contour.rb', line 227

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