module NSWTopo::Relief
Constants
- CREATE
- DEFAULTS
Public Instance Methods
get_raster(temp_dir)
click to toggle source
# 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()
click to toggle source
# File lib/nswtopo/layer/relief.rb, line 16 def margin { mm: 3 * @smooth } end