Module: NSWTopo::Relief

Includes:
ArcGISServer, DEM, Log, Raster, Shapefile
Defined in:
lib/nswtopo/layer/relief.rb

Constant Summary collapse

CREATE = 
%w[altitude azimuth factor sources yellow smooth median bilateral contours]
DEFAULTS = 
YAML.load <<~YAML
  altitude: 45
  azimuth: 315
  factor: 2.0
  sources: 3
  yellow: 0.2
  smooth: 4
  resolution: 5.0
  opacity: 0.3
YAML

Constants included from ArcGISServer

ArcGISServer::ERRORS, ArcGISServer::Error, ArcGISServer::SERVICE

Constants included from Shapefile

Shapefile::Error

Constants included from Log

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

Instance Method Summary collapse

Methods included from Raster

#create, #empty?, #filename, #render, #size_resolution, #to_s

Methods included from ArcGISServer

===, #arcgis_layer, check_uri, start

Methods included from Shapefile

===, #shapefile_layer

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

#get_raster(temp_dir) ⇒ Object 



20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
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
 
# File 'lib/nswtopo/layer/relief.rb', line 20

def get_raster(temp_dir)
  dem_path = temp_dir / "dem.tif"
  flat_relief = (Math::sin(@altitude * Math::PI / 180) * 255).to_i

  case
  when @path
    get_dem temp_dir, dem_path

  when @contours
    bounds = @map.bounds(margin: margin)
    txe, tye, spat = bounds[0], bounds[1].reverse, bounds.transpose.flatten
    outsize = (bounds.transpose.difference / @resolution).map(&:ceil)

    collection = @contours.map do |url_or_path, attribute_or_hash|
      raise "no elevation attribute specified for #{url_or_path}" unless attribute_or_hash
      options   = Hash == attribute_or_hash ? attribute_or_hash.transform_keys(&:to_sym).slice(:where, :layer) : {}
      attribute = Hash == attribute_or_hash ? attribute_or_hash["attribute"] : attribute_or_hash
      case url_or_path
      when ArcGISServer
        arcgis_layer url_or_path, margin: margin, **options do |index, total|
          log_update "%s: retrieved %i of %i contours" % [@name, index, total]
        end
      when Shapefile
        shapefile_layer source_path, margin: margin, **options
      else
        raise "unrecognised elevation data source: #{url_or_path}"
      end.each do |feature|
        feature.properties.replace "elevation" => feature.fetch(attribute, attribute).to_f
      end.reproject_to(@map.projection)
    end.inject(&:merge)

    log_update "%s: calculating DEM" % @name
    OS.gdal_grid "-a", "linear:radius=0:nodata=-9999", "-zfield", "elevation", "-ot", "Float32", "-txe", *txe, "-tye", *tye, "-spat", *spat, "-outsize", *outsize, "/vsistdin/", dem_path do |stdin|
      stdin.puts collection.to_json
    end

  else
    raise "no elevation data specified for relief layer #{@name}"
  end

  log_update "%s: generating shaded relief" % @name
  reliefs = -90.step(90, 90.0 / @sources).select.with_index do |offset, index|
    index.odd?
  end.map do |offset|
    (@azimuth + offset) % 360
  end.map do |azimuth|
    relief_path = temp_dir / "relief.#{azimuth}.bil"
    OS.gdaldem "hillshade", "-of", "EHdr", "-compute_edges", "-s", 1, "-alt", @altitude, "-z", @factor, "-az", azimuth, dem_path, relief_path
    [azimuth, ESRIHdr.new(relief_path, 0)]
  rescue OS::Error
    raise "invalid elevation data"
  end.to_h

  bil_path = temp_dir / "relief.bil"
  if reliefs.one?
    reliefs.values.first.write bil_path
  else
    blur_path = temp_dir / "dem.blurred.tif"
    blur_dem dem_path, blur_path

    aspect_path = temp_dir / "aspect.bil"
    OS.gdaldem "aspect", "-zero_for_flat", "-of", "EHdr", blur_path, aspect_path
    aspect = ESRIHdr.new aspect_path, 0.0

    log_update "%s: combining shaded relief" % @name
    reliefs.map do |azimuth, relief|
      [relief.values, aspect.values].transpose.map do |relief, aspect|
        relief ? aspect ? 2 * relief * Math::sin((aspect - azimuth) * Math::PI / 180)**2 : relief : flat_relief
      end
    end.transpose.map do |values|
      values.inject(&:+) / @sources
    end.map do |value|
      [255, value.ceil].min
    end.tap do |values|
      ESRIHdr.new(reliefs.values.first, values).write bil_path
    end
  end

  tif_path = temp_dir / "relief.tif"
  OS.gdalwarp "-co", "TFW=YES", "-s_srs", @map.projection, "-dstnodata", "None", bil_path, tif_path

  filters = []
  if @median
    pixels = (2 * @median + 1).to_i
    filters += %W[-channel RGBA -statistic median #{pixels}x#{pixels}]
  end
  if @bilateral
    threshold, sigma = *@bilateral, (60.0 / @resolution).round
    filters += %W[-channel RGB -selective-blur 0x#{sigma}+#{threshold}%]
  end
  if filters.any?
    log_update "%s: applying filters" % @name
    OS.mogrify "-virtual-pixel", "edge", *filters, tif_path
  end

  log_update "%s: rendering shaded relief" % @name
  vrt_path = temp_dir / "coloured.vrt"
  OS.gdalbuildvrt vrt_path, tif_path

  xml = REXML::Document.new vrt_path.read
  vrt_raster_band = xml.elements["VRTDataset/VRTRasterBand[ColorInterp[text()='Gray']]"]
  vrt_raster_band.elements["ColorInterp[text()='Gray']"].text = "Palette"
  color_table = vrt_raster_band.add_element "ColorTable"

  shade, sun = 90 * flat_relief / 100, (10 + 90 * flat_relief) / 100
  256.times do |index|
    case
    when index < shade
      color_table.add_element "Entry", "c1" => 0, "c2" => 0, "c3" => 0, "c4" => (shade - index) * 255 / shade
    when index > sun
      color_table.add_element "Entry", "c1" => 255, "c2" => 255, "c3" => 0, "c4" => ((index - sun) * 255 * @yellow / (255 - sun)).to_i
    else
      color_table.add_element "Entry", "c1" => 0, "c2" => 0, "c3" => 0, "c4" => 0
    end
  end

  vrt_path.write xml
  coloured_path = temp_dir / "coloured.tif"
  OS.gdal_translate "-expand", "rgba", vrt_path, coloured_path
  FileUtils.mv coloured_path, tif_path
  return @resolution, tif_path
end

#margin ⇒ Object 



16
17
18
 
# File 'lib/nswtopo/layer/relief.rb', line 16

def margin
  { mm: 3 * @smooth }
end