User Tools

Site Tools


r:r-tutorials:shadedterrain

~~DISCUSSION~~

Advanced Relief Shading

Please note there a tons of websites and tutorials online. Most of them are dealing with special techniques or software. This tutorial will just offer a way to simplify and integrate a lot of this “special techniques” for using them in an more or less automated way within R.

Things we cover

Cartographic relief representation is a central aspect of cartography. It deals with a two-dimensional representation of the continuously present three-dimensional real world structure on maps. One of the most famous of this classic cartographic art definitely was Eduard Imhof. So we will follow his approach in our tutorial. Nevertheless in our days digital hillshading as derived by elevation data is still an advanced task unless there are many straightforward solutions that work fine.

Web mapping is technically a bit different from producing a printable map. Furthermore there are fantastic maps served by WMS servers an based on the OpenStreetMap (OSM) project out there. You can mix them render them and even set up your own servers. Geofabrik is providing a toolset which you may use for finding always the appropriate map and informations. So why thinking about and why this hillshading tutorial?

Unfortunately most of the OSM maps are served as bitmap tiles or you can use them only online. Additionally in case of printing it let's say to letter or Din A4 it is hard to determine exactly what you want to see. Furthermore the integration of morphometric data is more ore less pretty basic. It does not integrate the flexibility of using different filter algorithms for smoothen or preparing the elevation data neither it is set up to derive an optimal visualization of morphometric informations according to scale and issues.

It seems obvious if you just want to pimp one the fly a nice but “flat” map of a common tile server or if you want to render your own OSM map using the great OpenStreetMap package or you want to produce your own hillshade design for a leaflet map there are many situation you may choose an comfortable data driven approach to create your personal hillshading model. Additionally you may use the alpha channel technique for all raster data you want to use as web mapping overlays.

In a geographic or let's call it better cartographic meaning hillshading is an complex topic.There is a lot of intuitive (geographic) informations that an observer can derive from a shaded relief, such as steepness, exposition, landforms and so on.

This tutorial will offer you the technical basics to develop your own shading. It uses a simplified approach after Imhof to merge four different insolation situations but it does not use the colorizing of Imhof. Nevertheless if you want to deal with filtering, coloring and alpha channels will have a perfect toolset to derive the best results for automated relief shading meeting your individual demand - without using photoshop etc..

We will cover:

  • getting the best Elevation data that is (automatically) available
  • using advanced filtering this data to derive a better visualization of the structures
  • getting hypsometric isolines
  • Imhof's approach as an example of a more sophisticated hillshading
  • using color tables
  • using alpha channel for transparency

Things we need

  • The tutorial is using Xubuntu 14.04 as operating system (OS). So it may be possible that some of the command line tools will work in a different way if you use Windows. If you are interested in setting up a GIS-shaped XubuntuGIS or if you want to have a all inclusive Windows GIS portable package you may look at XubuntuGIS Live-/Installations-DVD and/or R-WinGIS2Go portable-GIS bundle with fully integrated R. You will find more information at GIS OS distros(German only).

  • The raster package is the major player for dealing with spatial rasterized data. It is for all purposes the agent of choice. It is incredible powerful and comprehensive dealing with “reading, writing, manipulating, analyzing and modeling of gridded spatial data. The package implements basic and high-level functions and processing of very large files is supported”

    if (!require(raster)){install.packages(raster)}
    library(raster)
    
  • The SRTM (Shuttle Radar Topography Mission) sensor was mounted on a Space Shuttle and obtained Earth surface data by remote sensing technology utilizing a synthetic aperture radar (SAR). The data was processed and converted into height data aka Digital Elevation Model (DEM).

    Up to now There are at least 4 versions of SRTM data and version 4.1 available . Additionally since 2012 you can obtain the 30 meter data (only available for the cross referenced parts of the pathes) from the DLR data center . Additionally there are improved “non official” SRTM data providers. To make a long story short, just visit Christoph Hormann's Elevation data search page.

    You can choose almost all available sources for a meta search. By default it tries to offer you the SRTM data of Jonathan de Ferranti and provides you after the search the links to the download links of the data.

    Jonathan is doing an incredible good job. For the actual improvements have a look at Digital Elevation Data.

  • Sun et al, (2007) developed a feature-preserving mesh de-noising algorithm that smooths the surfaces of computer models of three dimensional objects such as those used in computer-aided design and graphics. This algorithm is perfect to reduce unwanted details of an DEM while preserving the main structures as ridges peaks and valleys.

    Please follow the installation instructions as listed under preparation.

  • gdalutils is nothing else then a very comfortable to use R packages that works as wrapper for the Geospatial Data Abstraction Library (GDAL) Utilities. So it not only simplify the GDAL calls but integrates them fully in R. In most cases it follows the naming conventions of the original functions unlike RSAGA allows for more R-like handling.

    if (!require(gdalUtils)){install.packages(gdalUtils)}
    library(gdalUtils)
    
  • GDAL is a translator library for raster and vector geospatial data formats. As a library, it presents a single raster abstract data model and vector abstract data model to the calling application for all supported formats. It also comes with a variety of useful command line utilities for data translation and processing.

    Traditionally GDAL used to design the raster part of the library, and OGR the vector part for Simple Features.

    You need a running binary installation of these Tolls. Please follow the Installation hints on the GDAL binary site to get the appropriate version for your system.

Things to do

Retrieving elevation data

First you need to download the elevation data of our area of interest. Due to the fact that it seems to be more interesting to generate a map for regions with a poor availability of high quality maps we choose the Cordillera area around the Cauca valley in Colombia. For your convenience you may use the prepared download link as generated by the imagico search page. Note the data is taken from the SRTM data distribution of Jonathan de Ferranti. Especially in high mountains he is providing the best SRTM data that is available.

# load raster package
library(raster)
# load gdalUtils
library(gdalUtils)

# define working folder
root.dir <- "~/workspace/tut"      # root folder has to exist!
working.dir <- "r-tut"         # working folder if not exist will be created
# create working directory -> set as working directory ignore warnings if dir exist
dir.create(file.path(root.dir, working.dir),showWarnings = FALSE)
setwd(file.path(root.dir, working.dir))


# get the data as suggested above note this unprojected data i.e. geographic coordinates
download.file('http://www.viewfinderpanoramas.org/dem3/B18.zip','B18.gz')

# unzip the data to a folder named srtm-data -o overwrite -j don't use pathes
system("unzip -j -o B18.gz ")

# generate a list of all srtm files in this folder
file.list<-paste(list.files(getwd(),("hgt")),collapse = ' ')

# merge them with gdal_merge to one file in ASC format
system(paste('gdal_merge.py -init 0 -o dem.tif -of Gtiff', file.list))

# cut specific are of interest (aoi)
gdal_translate('dem.tif','demcut.tif', projwin='-76. 5.5 -75. 4.0')
# alternatively the system call with the command line syntax
# system('gdal_translate -projwin -76. 5.5 -75. 4.0 -of GTiff dem.tif demcut.tif')

Smoothing the elevation data

It is difficult to find optimal balanced rules for relief shading. Nevertheless it is absolutely necessary to derive a good mixture of information and generalization from digital elevation data without to much detailed informations but still keeping the interesting morphometric informations. Up to my knowledge the best algorithm dealing with this opposite goals is the mdenoise algorithm as provided by Sun. You have to deal with the parameters to derive the result you want but you may start with the suggested numbers. Be careful the filtering needs computing time and the area is restricted in size due to hardware limitations.


# mdenoise only works with the ESRI ASCII data format ('asc') so we convert the elevation data to asc
system('gdal_translate -of AAIGrid -ot Int16 dem.tif dem.asc ')

# you may use  mdenoise for benchmarking because it is an extremely cpu time consuming task 
# some numbers for the chosen data set
# Read Model...    17.081 seconds
# Denoising Model...   532.214 seconds
# Saving Model...    17.226 seconds
# you have to play with the -n and -t parameters so you may look at the homepage for further information 
system('~/progs/mdenoise/mdenoise -i dem.asc -e -n 20 -t 0.89  -o dem-denoise.asc')

# project towards EPSG:102033 - South_America_Albers_Equal_Area_Conic
# using cubicspline to smooth it because it is for visualisation only
# and convert it to geotiff note we are using a mixed description for referencing in EPSG Code and PROJ4 string. if you have problems you should use the PROJ4 string
system('gdalwarp -overwrite -s_srs EPSG:4326 -t_srs "+proj=aea +lat_1=-5 +lat_2=-42 +lat_0=-32 +lon_0=-60 +x_0=0 +y_0=0 +ellps=aust_SA +units=m +no_defs" -r cubicspline -of GTiff dem-denoise.asc dem-denoise.tif')

Extracting altitude isolines from the filtered DEM

For more information and to be able to use the map for orientation we extract altitudinal isolines using the gdal_contour tool. Note due to the fact that it is not directly supported by gdalUtils we just use the system call.

# calculate contour lines every 20m and every 50m
# the parameter -3d extract the altitude data from the DEM and assigns it to each node
# calculate contour lines every 20m and every 50m you can not overwrite existing files!
system('gdal_contour -a height dem-denoise.tif contour20.shp -3d -i 20.0')
system('gdal_contour -a height dem-denoise.tif contour100.shp -3d -i 100.0')
system('gdal_contour -a height dem-denoise.tif contour500.shp -3d -i 500.0')

Calculating hillshade models

Hillshading algorithm are very popular. In the code below we are taking the Imhof approach into account that we need different views on the morphometric structures. So we use for different azimuth angels. Later we will merge them. Additionally we are inverting the grey scales and add an so called alpha channel for opacity. While the resulting hillshade file is used as an overlay (in any webbrowser or in QGIS…) it will be multiplied with the other layers. Due to the inverted grey scale and the opacity it will nicely merge. You can play around with the color ramp file. Please look also at the further readings to get an idea how and why to do so.

# generating a hillshade base map

# generate grey scaled alpha channels to the single hillshade
# we are using 4 different azimuth for a more structured 3D relief
gdaldem('hillshade', 'dem-denoise.tif','hs340_45.tif', az=340, alt=45,  of='GTiff')
gdaldem('hillshade', 'dem-denoise.tif','hs160_45.tif', az=160, alt=45,  of='GTiff')
gdaldem('hillshade', 'dem-denoise.tif','hs295_45.tif', az=295, alt=45,  of='GTiff')
gdaldem('hillshade', 'dem-denoise.tif','hs340_45.tif', az=250, alt=45,  of='GTiff')

# generate grey scaled hillshades to demonstrate the use of:
# - alpha channels and color-tables
# - output_raster flag
# note we put alpha and output_raster both to false because do not need it here 
# if you need/want it put it on TRUE
gdaldem('color-relief',input_dem='hs16045.tif', output='a16045.tif', of='GTiff',output_Raster=FALSE, col='grey2.txt', alpha=FALSE, nearest_color_entry=TRUE)
gdaldem('color-relief',input_dem='hs34045.tif', output='a34045.tif', of='GTiff',output_Raster=FALSE, col='grey2.txt', alpha=FALSE, nearest_color_entry=TRUE)
gdaldem('color-relief',input_dem='hs29545.tif', output='a29545.tif', of='GTiff',output_Raster=FALSE, col='grey2.txt', alpha=FALSE, nearest_color_entry=TRUE)
gdaldem('color-relief',input_dem='hs25045.tif', output='a25045.tif', of='GTiff',output_Raster=FALSE col='grey2.txt', alpha=FALSE, nearest_color_entry=TRUE)

# read resulting GeoTiff output files as raster object
# NOTE we are not using the gdalUtils option  output_Raster=TRUE  due to the fact that it produces a rasterBrick
a160<-raster('a16045.tif')
a340<-raster('a34045.tif')
a295<-raster('a29545.tif')  
a250<-raster('a25045.tif')

# merge them equally weighted
dem.a<-(a160+a340+a295+a250)/4

# mean filter for smoothing
dem.a.f<-focal(dem.a, w=matrix(1/9,nrow=3,ncol=3))

# write the derived hillshading rasters to a geotiff file
# without alpha channel
writeRaster(dem.a, filename='dem.tif',overwrite=TRUE)
writeRaster(dem.a.f, filename='dem_meanfilter.tif',overwrite=TRUE)

# assign alpha channels via gdaldem
# alpha value are in the forth col of the grey2.txt file
gdaldem('color-relief',input_dem='newdem.tif', output='dem_alpha.tif', of='GTiff',output_Raster=TRUE, col='grey1.txt', alpha=FALSE, nearest_color_entry=TRUE)
gdaldem('color-relief',input_dem='newdem_f.tif', output='dem_alpha_meanfilter.tif', of='GTiff',output_Raster=TRUE, col='grey1.txt', alpha=FALSE, nearest_color_entry=TRUE)



The grey scale alpha channel ramp file

For the above example you may use the following grey2.txt file. You can alter it in any way. The structure is simple and documented on the GDAL site. First value is the grid value of the file (index, altitude, whatever), the next three values are the RGB code and the last one is the value for the alpha channel (scaled from 0-255).

255	255	255	255	128
253	254	254	254	128
252	253	253	253	128
251	252	252	252	128
250	251	251	251	128
249	250	250	250	128
248	249	249	249	128
247	248	248	248	128
246	247	247	247	128
245	246	246	246	128
244	245	245	245	128
243	244	244	244	128
242	243	243	243	128
241	242	242	242	128
240	241	241	241	128
239	240	240	240	128
238	239	239	239	128
237	238	238	238	128
236	237	237	237	128
235	236	236	236	128
234	235	235	235	128
233	234	234	234	128
232	233	233	233	128
231	232	232	232	128
230	231	231	231	128
229 	230	230	230	128
228	229	229	229	128
227	228	228	228	128
226	227	227	227	128
225	226	226	226	128
224	225	225	225	128
223	224	224	224	128
222	223	223	223	128
221	222	222	222	128
220	221	221	221	128
219	220	220	220	128
218	219	219	219	128
217	218	218	218	128
216	217	217	217	128
215	216	216	216	128
214	215	215	215	128
213	214	214	214	128
212	213	213	213	128
211	212	212	212	128
210	211	211	211	128
209	210	210	210	128
208	209	209	209	128
207	208	208	208	128
206	207	207	207	128
205	206	206	206	128
204	205	205	205	128
203	204	204	204	128
202	203	203	203	128
201	202	202	202	128
200	201	201	201	128
199	200	200	200	128
198	199	199	199	128
197	198	198	198	128
196	197	197	197	128
195	196	196	196	128
194	195	195	195	128
193	194	194	194	128
192	193	193	193	128
191	192	192	192	128
190	191	191	191	128
189	190	190	190	128
188	189	189	189	128
187	188	188	188	128
186	187	187	187	128
185	186	186	186	128
184	185	185	185	128
183	184	184	184	128
182	183	183	183	128
181	182	182	182	128
180	181	181	181	128
179	180	180	180	128
178	179	179	179	128
177	178	178	178	128
176	177	177	177	128
175	176	176	176	128
174	175	175	175	128
173	174	174	174	128
172	173	173	173	128
171	172	172	172	128
170	171	171	171	128
169	170	170	170	128
168	169	169	169	128
167	168	168	168	128
166	167	167	167	128
165	166	166	166	128
164	165	165	165	128
163	164	164	164	128
162	163	163	163	128
161	162	162	162	128
160	161	161	161	128
159	160	160	160	128
158	159	159	159	128
157	158	158	158	128
156	157	157	157	128
155	156	156	156	128
154	155	155	155	128
153	154	154	154	128
152	153	153	153	128
151	152	152	152	128
150	151	151	151	128
149	150	150	150	128
148	149	149	149	128
147	148	148	148	128
146	147	147	147	128
145	146	146	146	128
144	145	145	145	128
143	144	144	144	128
142	143	143	143	128
141	142	142	142	128
140	141	141	141	128
139	140	140	140	128
138	139	139	139	128
137	138	138	138	128
136	137	137	137	128
135	136	136	136	128
134	135	135	135	128
133	134	134	134	128
132	133	133	133	128
131	132	132	132	128
130	131	131	131	128
129	130	130	130	128
128	129	129	129	128
127	128	128	128	128
126	127	127	127	128
125	126	126	126	128
124	125	125	125	128
123	124	124	124	128
122	123	123	123	128
121	122	122	122	128
120	121	121	121	128
119	120	120	120	128
118	119	119	119	128
117	118	118	118	128
116	117	117	117	128
115	116	116	116	128
114	115	115	115	128
113	114	114	114	128
112	113	113	113	128
111	112	112	112	128
110	111	111	111	128
109	110	110	110	128
108	109	109	109	128
107	108	108	108	128
106	107	107	107	128
105	106	106	106	128
104	105	105	105	128
103	104	104	104	128
102	103	103	103	128
101	102	102	102	128
100	101	101	101	128
99	100	100	100	128
98	99	99	99	128
97	98	98	98	128
96	97	97	97	128
95	96	96	96	128
94	95	95	95	128
93	94	94	94	128
92	93	93	93	128
91	92	92	92	128
90	91	91	91	128
89	90	90	90	128
88	89	89	89	128
87	88	88	88	128
86	87	87	87	128
85	86	86	86	128
84	85	85	85	128
83	84	84	84	128
82	83	83	83	128
81	82	82	82	128
80	81	81	81	128
79	80	80	80	128
78	79	79	79	128
77	78	78	78	128
76	77	77	77	128
75	76	76	76	128
74	75	75	75	128
73	74	74	74	128
72	73	73	73	128
71	72	72	72	128
70	71	71	71	128
69	70	70	70	128
68	69	69	69	128
67	68	68	68	128
66	67	67	67	128
65	66	66	66	128
64	65	65	65	128
63	64	64	64	128
62	63	63	63	128
61	62	62	62	128
60	61	61	61	128
59	60	60	60	128
58	59	59	59	128
57 	58	58	58	128
56	57	57	57	128
55	56	56	56	128
54	55	55	55	128
53	54	54	54	128
52	53	53	53	128
51	52	52	52	128
50	51	51	51	128
49	50	50	50	128
48	49	49	49	128
47	48	48	48	128
46	47	47	47	128
45	46	46	46	128
44	45	45	45	128
43	44	44	44	128
42	43	43	43	128
41	42	42	42	128
40	41	41	41	128
39	40	40	40	128
38	39	39	39	128
37	38	38	38	128
36	37	37	37	128
35	36	36	36	128
34	35	35	35	128
33	34	34	34	128
32	33	33	33	128
31	32	32	32	128
30	31	31	31	128
29	30	30	30	128
28	29	29	29	128
27	28	28	28	128
26	27	27	27	128
25	26	26	26	128
24	25	25	25	128
23	24	24	24	128
22	23	23	23	128
21	22	22	22	128
20	21	21	21	128
19	20	20	20	128
18	19	19	19	128
17	18	18	18	128
16	17	17	17	128
15	16	16	16	128
14	15	15	15	128
13	14	14	14	128
12	13	13	13	128
11	12	12	12	128
10	11	11	11	128
9	10	10	10	128
8	9	9	9	128
7	8	8	8	128
6	7	7	7	128
5	6	6	6	128
4	5	5	5	128
3	4	4	4	128
2	3	3	3	128
1	2	2	2	128
0	0	0	0	128

Now you can use it for instance to extract data from huge files:

 chmod +x ECWJP2SDKSetup_5.2.1.bin
 
 # and start the setup NOTICE desktop-read-only is the license you need
 ./​ECWJP2SDKSetup_5.2.1.bin ​
 # it will be installed at ~/​hexagon/​ERDAS-ECW_JPEG_2000_SDK-5.2.1 copy it to the correct /usr/local folder
 
 sudo cp -r ~/​hexagon/​ERDAS-ECW_JPEG_2000_SDK-5.2.1 /usr/local/
 sudo ln -s /​usr/​local/​ERDAS-ECW_JPEG_2000_SDK-5.2.1/​Desktop_Read-Only/​lib/​x64/​release/​libNCSEcw.so /​usr/​local/​lib/​libNCSEcw.so
 sudo ldconfig
 
 # now get the libgdal plugin
 wget https://​launchpad.net/​~ubuntugis/​+archive/​ubuntugis-unstable/​+files/​libgdal-ecw-src_1.10.0-1~precise4_all.deb
 # unpack it
 ar vx libgdal-ecw-src_1.10.0-1~precise4_all.deb
 # and untar the data archive
 tar -xzf data.tar.gz
 # copy to /usr/local
 sudo cp usr/​src/​libgdal-ecw-1.10.0.tar.gz /usr/src/
 # and usr/bin
 sudo cp    usr/​bin/​gdal-ecw-build /usr/bin/
 # now build it
 sudo gdal-ecw-build ~/​hexagon/​ERDAS-ECW_JPEG_2000_SDK-5.2.1/​Desktop_Read-Only
 
 # make a copy  of the folder according to the expected version number
 
  sudo mkdir /​usr/​lib/​gdalplugins/1.11
  cd /​usr/​lib/​gdalplugins/1.10
  cp gdal_ECW_JP2ECW.so /​usr/​lib/​gdalplugins/1.11 ​
 
 # check it
  ​gdalinfo --formats | grep -i ECW
  ECW (rw+): ERDAS Compressed Wavelets (SDK 5.1)
  ​JP2ECW (rw+v): ERDAS JPEG2000 (SDK 5.1)
 
 #!/bin/bash
 # get all files in current directory with the ecw extension
 for filename in  *ecw
 do
 # generate outputfilname
   outfilename=$filename'​.tif'echo "​Processing $filename file..."#  extract a certain area an write it to jpeg embedded tif
   gdal_translate -a_nodata 0 -projwin -21543.0 252606.0 -5846.0 235465.0 -of GTiff -co COMPRESS=JPEG filenameoutfilename
 done

Things of further interest

If you want to get an starting point for cartographic hillshading have a look at the the technically outdated but in a cartographic manner comprehensive Relief Shading page. For a more technical approach have a look at Ideas and Techniques about Relief Presentation on Maps.

According to R and OpenStreetMap at least the following sites should be acknowledged:

Christoph Reudenbach 2014/12/31 15:21

r/r-tutorials/shadedterrain.txt · Last modified: 2018/12/23 19:46 (external edit)