Project

General

Profile

Download (2.33 KB) Statistics
| Branch: | Revision:
1
# GDAL commands to generate hillshade images using the different fused
2
# DEMs we've produced in the 60N Canada boundary region.
3
#
4
# Also included below is a bit of R code to generate a PNG of quick
5
# hillshade visuals of two 1-degree wide swaths (one in the more
6
# mountainous west, and one in the flatter east) spanning the 60N
7
# boundary in Canada, for several different fused layers.
8
#
9
# Jim Regetz
10
# NCEAS
11

    
12
export DEMDIR=~/media/temp/terrain/dem
13
export HSDIR=~/media/temp/terrain/hillshade
14

    
15
gdaldem hillshade -s 111120 $DEMDIR/fused_300straddle.tif \
16
  $HSDIR/fused_300straddle_h.tif
17
gdaldem hillshade -s 111120 $DEMDIR/fused_300straddle_rampexp.tif \
18
  $HSDIR/fused_300straddle_rampexp_h.tif
19
gdaldem hillshade -s 111120 $DEMDIR/fused_300straddle_predexp.tif \
20
  $HSDIR/fused_300straddle_predexp_h.tif
21
gdaldem hillshade -s 111120 $DEMDIR/fused_300straddle_predgau.tif \
22
  $HSDIR/fused_300straddle_predgau_h.tif
23
gdaldem hillshade -s 111120 $DEMDIR/fused_300straddle_blendgau.tif \
24
  $HSDIR/fused_300straddle_blendgau_h.tif
25
gdaldem hillshade -s 111120 $DEMDIR/fused_300straddle_enblend.tif \
26
  $HSDIR/fused_300straddle_enblend_h.tif
27

    
28
# Use R to generate a PDF with some visuals 
29

    
30
echo '
31
  library(raster)
32
  hsdir <- "/home/regetz/media/temp/terrain/hillshade"
33
  
34
  h.uncor <- raster(file.path(hsdir, "fused_300straddle_h.tif"))
35
  h.re <- raster(file.path(hsdir, "fused_300straddle_rampexp_h.tif"))
36
  h.enblend <- raster(file.path(hsdir, "fused_300straddle_enblend_h.tif"))
37
  h.bg <- raster(file.path(hsdir, "fused_300straddle_blendgau_h.tif"))
38
  
39
  window1 <- extent(-135, -134, 59.875, 60.125)
40
  window2 <- extent(-118, -117, 59.875, 60.125)
41
  
42
  png("boundary-hillshade.png", height=10, width=7.5, units="in", res=600)
43
  par(mfrow=c(4,2), mar=c(2,2,3,0))
44
  plot(crop(h.uncor, window1), legend=FALSE)
45
  mtext("Simple fuse (hillshade)", adj=0, cex=0.8)
46
  plot(crop(h.uncor, window2), legend=FALSE)
47
  plot(crop(h.re, window1), legend=FALSE)
48
  mtext("North exponential ramp (hillshade)", adj=0, cex=0.8)
49
  plot(crop(h.re, window2), legend=FALSE)
50
  plot(crop(h.enblend, window1), legend=FALSE)
51
  mtext("Multiresolution spline (hillshade)", adj=0, cex=0.8)
52
  plot(crop(h.enblend, window2), legend=FALSE)
53
  plot(crop(h.bg, window1), legend=FALSE)
54
  mtext("Gaussian weighted average (hillshade)", adj=0, cex=0.8)
55
  plot(crop(h.bg, window2), legend=FALSE)
56
  dev.off()
57
' | Rscript --vanilla -
(20-20/23)