Version4, 01/05/2023
ViewShed Analysis Documentation
Eitam Arnon eitamarnon@[Link]
Orr Spiegel orrspiegel@[Link]
Contents
1. Overview:................................................................................................................................................1
2. Installation ..............................................................................................................................................2
3. Calculation of viewshed layers ...............................................................................................................3
a. Preparing the system and required data............................................................................................3
b. Invoking the calculation .....................................................................................................................4
c. Editing LOS files after calculation .......................................................................................................6
4. Presentation ...........................................................................................................................................8
a. Loading a precalculated layer.............................................................................................................8
b. Cumulative Viewshed .........................................................................................................................9
c. Subtractive Viewshed .......................................................................................................................10
d. LOS analysis ......................................................................................................................................11
5. Obtaining DEM data .............................................................................................................................13
6. Earth curvature .....................................................................................................................................13
1. Overview:
The system is intended to calculate and present a line of sight (LOS) analysis from a set of pre-
defined points (Antennas) to any point in a chosen geographical region (Transmitter). It provides
several analysis tools such as the display of cumulative and intersection of LOS area from several
antennas, as will be detailed below. The results Depend on the height of the antenna, the height
of the transmitter and a Digital Elevation Map (DEM) file that represent the surface height in the
region. This tool can either consider or neglect the effect of earth curvature. The tool was
designed to support both terrestrial antenna arrays and aquatic antenna array. The package
includes several actions:
a) Calculation of a Viewshed layer, a Boolean raster layer for a specific antenna and a
transmitter height, it set a TRUE (visible) or FALSE (invisible) value to all points in a
region (an area predefined by the user). This method includes a parallel calculation to
increase the efficiency of the process.
b) Visualization of the cumulative viewshed. Using a set of precalculated raster layers and a
choice of active antennae (this is the plural form for “antenna”, it is a word of Latin
1
origin), the function sums the number of observant antennae at any point and present
this value on an interactive map.
c) Visualization of intersecting line of sight. Using a set of precalculated raster layers, a
choice of observant antennae and blind antennae, the function identifies the points that
obey these visibility conditions and present these on an interactive map.
d) The function can allow the user to choose an antenna and point on the map and present
the surface height under a line connecting these two points. It also presents the detailed
DEM at the vicinity of the line. This analysis allows the investigation of the obstacles
between any two points to explain the LOS results
2. Installation
To run the program, one must have an installed R program (updated to 4.0.0 or newer) and
RStudio IDE (interactive development environment), instructions for these steps can be found in
[Link]
Using the application can be divided into two steps; a) Calculating viewshed layers, that requires
a list of antennas location and height above ground, and some R programing capabilities and b)
presenting, viewing, and analyzing the precalculated LOS layers, this step does not require any R
programing capabilities. Our recommended approach is that step a can be done once for a
specific array of antennas by a system-manager and shared among a group of users that can
apply step b, without any deep knowledge of the system or R.
It is recommended to save the entire application in a single directory (commonly named
“Viewshed”), containing a R project “Viewshed” (it is recommended to open the RStudio
application using this R project file, this will define the active directory to search locally). The
folder must contain several subdirectories (to allow the correct path for searching files) as
follows (within the distribution, each directory contains some example data from the Harod
valley, northern Israel):
• DEM_Files – a file containing the DEM files needed to compute the LOS analysis
(needed only for the calculation step).
• ANT_Table_Files - a directory containing a csv file with antenna coordinates and
height above ground analysis (needed only for the calculation step).
• LOSData – a directory containing the layers produced in the calculation step (as
many as were created by the system operator). Each calculation setup (for a list of
antennas and target height) produces a single [Link] file that contains 3 items (
'[Link]' - a data frame with antennas details,'LOSLayers' - a raster with Boolean
coverage and ,'DEM' – the area DEM data) . The directory also contains the icons for
antenna ([Link])
• Functions – a directory containing the set of functions required to calculate and
present the analysis.
• Documentation – a directory containing this file (this subdirectory is not
mandatory)
2
3. Calculation of viewshed layers
a. Preparing the system and required data.
This is a preparation step; the end-user can skip it to chapter 4 – Presentation. To initiate the
calculation, one should open the script file “computeViewshed.R” (located in “functions”
subdirectoryThe following paragraph will detail the fundamental functions and procedures a user
should invoke (code is accentuated distinguish font) to complete the calculation, these
commands are already part of the script and the user needs only to change specific input
parameters that will be highlighted in the following lines. The script file contains a reasonable
sequence of commands that will allow the calculation of the viewshed layers, not all of them are
mandatory:
• source('functions/Viewshed.R')
o All functions are in the “Viewshed.R” file and are attached by this
command.
• setEnv()
o This function source external libraries and define several constants, several
warnings usually appear in this process.
o If an error appears: Error in library("XXX") : there is no package called ‘XXX’,
the user is required to install the following package:
[Link]("XXX") (the code is designated to install all
required packages automatically but still these errors sometimes appear in
the first operation of the code).
o It is required to update your R version to 4.0.0 or later!
• box <- [Link](lat=c(32.38,32.68),lon=c(35.27, 35.68))
o The coverage map is calculated only within a predefined region, here the
geographical coordinates of two diagonal corners of the region should be
set.
o Note that larger area requires larger computational effort and therefore
longer calculation times. ( see section b, Invoking the calculation for
calculation time approximations)
o The coordinate system should be WGS84 (EPSG:4326, World Geodetic
System 1984)
• DEM_name<-"DEM_Files/DEM_Harod.tif"
o If the user has the elevation data of the region (digital elevation map- DEM )
on a geographic grid, he can specify it here, other wise the data can be
freely download from the web by specifying type="web" in the
getDEM command ( this usually provide a DEM with 3 arc resolution ( ~ 90
m ) from the SRTM mission).
o The DEM should be given in TIFF format, see Obtaining DEM data section
for more details.
o The DEM coordinate system must be WGS84 (EPSG:4326, World Geodetic
System 1984)
• ANTfilename <- " ANT_Table_Files/[Link]"
o Path to the antenna location csv file. The file should contain antenna
locations and must include the following variables:
3
▪ "ANTName" – name in English letters, numbers and spaces are
allowed.
▪ "ID" – an integer number specifying the antenna.
▪ "LAT", "LON" – geographical coordinates (geographical projection,
wgs84 datum)
▪ "Towerheight" – antenna height above ground in meters
o The file can also include “AntRadius” variable – the reception range of each
antenna (in meters)
• DEM <- getDEM(type="file",box,filename=DEM_name)
o Retrieve and cut DEM from the tiff file, in case the user wants to reduce the
resolution, an argument resoluton = XXX where XXX is the required
resolution in meters can be added to the getDEM function.
o In a case where the user doesn’t have a high-resolution DEM file, the DEM
can be obtained from the web: DEM <- getDEM(type="web",box),
this usually provide a DEM with 3 arc resolution ( ~ 90 m )
o The accuracy and resolution of the DEM depends on its source. The same is
true for question whether the DEM include buildings and trees (the most
widespread source, STRM was create by radar signal to the international
space-station, thus it includes buildings (although not accurately because of
its resolution) and only dense vegetation.
o Note that the provided resolution has a dramatic effect on the calculation
time and file size. Therefore, before dispatching the calculation with
maximal resolution it is strongly recommended to begin with a low-
resolution run.
• [Link] <- setANTS(ANTfilename,DEM)
o Reads and prepares the antenna details file. The file should contain all
antenna locations and must include the following variables:
▪ "ANTName" – name in English letters, numbers and spaces are
allowed.
▪ "ID" – an integer number specifying the antenna.
▪ "LAT", "LON" – geographical coordinates (geographical projection,
wgs84 datum)
▪ "Towerheight" – antenna height above ground in meter
o The file can also include “AntRadius” variable – the reception range of each
antenna (in meters)
o The function returns a SpatialPointsDataFrame object (see sp package for
details) in the geographic, wgs projection, with the added variable
“GroundASL” which is the ground altitude at the antenna position (in
meters)
• viewSetup([Link],DEM)
o presenting the system setup: antennas, region map, and DEM. Click an
antenna to see its details.
b. Invoking the calculation
After setting the required variables the system is ready for the calculation of the
line-of-sight layers (LOSlayer) and saving it. As mentioned before, the provided
4
resolution and region-size has a dramatic effect on the calculation time and file
size. Therefore, before dispatching the calculation with maximal resolution it is
strongly recommended to begin with a low-resolution, small region run.
Gross approximation of calculation time:
• While calculating Harod system we used a laptop (i5 core)
• The size of the region was 40km*40km
• when using a 500m resolution DEM, calculating a layer took 1.5 seconds (for a single
transmitter height and a single antenna)
• when using a 30m resolution DEM, calculating a layer took 8 minutes (for a single
transmitter height and a single antenna)
• Note that the R studio IDE signifies that the command is working by a red sign in the
top of the console window. A user can wait or cancel the command by hitting “esc”.
The calculation can be launched either as serial calculation or parallel calculation
(whenever a large computational power is available) the calculation run for a set
of defined antennae and a set of transmitter heights and saves several required
files for each transmitter height. This is being done once per setup and can take
up to several hours, (but then you do not need to do it again )
The serial code:
• SerialComputeViewShed(layername=Layername,DEM,[Link],transAlts
,ANTlist=ANTlist,includeCurv=F,seaLevel = NA)
The parallel code:
• ParallelComputeViewShed(layername=Layername,DEM,[Link],transAl
ts,ANTlist=ANTlist,includeCurv=F,seaLevel = 0)
The input variables for both these functions are:
• layername="XXX"
o A name for the file produced by this procedure (the same files will be used
to present the LOS analysis). TransAlt0.3m_ XXX_vsf.RData
• DEM
o The DEM raster variable.
• [Link]
o The antenna detailes variable as produced by the setANTS procedure
• ANTlist <- c("1","2","3"),
o A list of antenna ID’s to run the calculation on, The default is running the
calculation for all antennas in the antenna csv file (can be specified explicitly
by: ANTlist <-[Link]([Link]$ID))
• transAlts <- c(-2,0.3, 2, 5,10)
o list of transmitter heights to calculate, each height will be calculated for all
antennas and saved as a single file.
o Since the DEM has a finite accuracy and resolution we do not recommend
to use 0 height.
5
o When working on aquatic systems ( seaLevel parameter is not NA),
negative transAlts can also be specified, this option interprets the
value as meters below sea level see Fig 1 for details
• includeCurv
o Logical parameter (TRUE / FALSE ).
o Determine weather to include earth curvature in the LOS calculation
o Default value is FALSE
• seaLevel
o Numeric value, sets the water surface elevation for aquatic systems
o Affects the transmitter Altitude parameter and the invalid areas parameter,
ignored when set to NA (default value)
o See Figure 1 for detailes
it is recommended to run the following line before invoking the parallel
calculation.
• detectCores()
o detects how many cores are available on your machine.
• registerDoParallel(3)
o assigns X processes for the following calculation.
Figure 1: the transAlts parameter effect for aquatic system. When seaLevel is set to any number (not NA ) the system is
treated as aquatic: positive transAlts values will be interpreted as meters above the seabed (green line) negative values
will be interpreted as meters below water surface ( yellow line). Areas where the LOS calculation is invalid (marked as
red rectangles) will be marked as black patches during any visualization.
c. Editing LOS files after calculation
The calculation of a complete set of antennas might become demanding,
especially at large region sizes and high DEM resolutions. Therefore the ability to
edit a LOSlayers set is required, this ability includes the addition, removal as well
as changing names (antenna name and ID) without recalculation of viewshed.
substituting a single antenna:
6
• AddRmAnts(final_name,base_str_name,add_str_name,ANTID2rm,ANTID2add)
o The function edit the base_str_name file, (note that you need only
to provide the kernel name, without the “ANTS”, “LOSLayers“ or
“DEM” prefix), it removes ANTID2rm antennas and adds ANTID2add
from the file add_str_name
o See code snippet below for example.
o The function will not protect you from combining files that were
calculated for different transmitter height (It will give you a
warning if both strings begin with different name (“TransAlt”).
o The function cannot overwrite the previous files.
• listLOSFiles()
o List of viewshed files available in the LOSDATA directory, it trims
the prefix
• listAnts(str_name)
o List the calculated antennas in a set of LOS files for a given file-
name
• An example editing procedure:
listLOSFiles() print available files
str_name <- "TransAlt3m_Res900a1-6" choosing a file to edit
(an explanation for the name: the “TransAlt3m” represents the height of the
transmitter and “Res900a1_6” was a user specific name that represent the
resolution of the DEM and the active antennas)
str_name_new <- "TransAlt0.3m_Res900a11-16" choosing a source file
for new antennas.
str_name_out <- "TransAlt0.3m_Res900new-a11-16" setting a new name for
the edited file.
listAnts(str_name) listing the antennas in the original file
listAnts(str_name_new) listing the antennas in the additional file
AddRmAnts(final_name=str_name_out,base_str_name=str_name,…
add_str_name=str_name_new,ANTID2rm = c(2,3),ANTID2add = c(15,12))
activating the editing for the previous choices
• UpdateAntnames(final_name,base_str_name,initial_ANTID,
newANTID,newANTname)
o Updates the names and ID’s of antennas.
o Required parameters:
▪ Final name to save the adited LOS files (final_name)
▪ Source file name (base_str_name)
▪ List of initial ID’s (initial_ANTID)
▪ new ID’s, after the editing (newANTID)
▪ New names (newANTname)
▪ newANTID and newANTname can also be omitted if only the ID or
only the name should be changed.
7
o Be careful not to give an ID that already exists!
4. Presentation
After the layers have been calculated, files are automatically saved to the
designated folder. The application is ready for presenting the data.
It is recommended to share the entire directory (including a single R project file)
with the possible users of the application. Thus, any user of the program with
installed R program can use the application by simply executing the command:
source('functions/ShinyViewshed.R').
Basically, upon launching the program a “load data” tab will be opened, the user
should choose a LOSdata set from the drop-menu and hit “Load file” button. The
user can than choose the different analysis tab: “Cumulative viewshed”,
“Subtractive viewshed”, or “Line of sight” as will be detailed below!
To terminate the process a user is just required to press “esc” in the R console,
and to close the web viewer, no saving is required. This is sometimes required
after the system has been ideal for stuck.
a. Loading a precalculated layer
Upon launching this command a web browser window will open:
Figure 2: A screenshot, the user interface after launching the program.
1. The user should choose a set of precalculated antenna file from the
available list and hit “Load file”
2. The system will present the calculated region and antenna, clicking an
antenna will present its details:
8
Figure 3: An overview of the loading data tab, after a LOSdata file was chosen. A red square was added to the top of the screenshot to
signify the tab selection.
3. The “Map vertical size on screen” input allows the user to present the
map on the eniter screen for larger screens
After verifying the basic data is as required, the user can navigate
between the different tabs and present / analyze the data.
b. Cumulative Viewshed
Here the user can choose the antennae that contribute and present the
number of visible antennae at any point in the region (this calculation is
valid for the pre-defined transmitter height)
The map will update when the “Recalc” button is clicked (top red square
on Fig. 3)!
The participating antennae are marked as stars on the map.
The output can be exported (to other application, e.g. google earth) in a
kmz format to “file name” by clicking “Save to kmz” at the bottom of the
page (bottom red square on Fig. 3)!
9
Figure 4: A cumulative viewshed example, after several antennae have been chosen and recalculted. Red squares were added to the
screenshot: top to signify the “Recalc” option and bottom to signify the ”save to kmz” option.
c. Subtractive Viewshed
Here the user can present the intersection between in-site (with LOS) and
blind (without LOS) antennae. The user can choose multiple antennas for
each side, but the list has to be exclusive (none of them can appear in
both lists). Note that this calculation is true for the pre-defined
transmitter height.
The map will update when the “Recalc” button is clicked, area that agree
with the defined condition will be colored yellow!
The in-site antennae are marked as stars ( ) on the map, while the blind
are marked by circled cross ( )
The output can be exported in a kmz format to “file name” by clicking
“Save to kmz” at the bottom of the page!
10
Figure 5: An subtractive viewshed example, after several antennae have been chosen and recalculated, the yellow patches mark the
areas were all conditions agree.
d. LOS analysis
This tab allows some detailed analysis of the area and terrain. It plots the
ground elevation under a line between two points. It also presents the
DEM values in a square around this line (with user-chosen shoulders), this
option allows a high resolution inspection of the input DEM.
• Both points of the line can be inserted manually as a geographic wgs84
coordinates in the “LAT” and “LON” text input fields. Their elevation above
ground can also be inserted manually.
• A second option for Antenna position input is to choose a predefined
antenna from the underlying file.
• A second option for Transmitter position input is via the icon ( ) in the
map. This icon will allow the user to point a location and will copy its
coordinates and ground altitude to the corresponding fields.
• The “DEM shoulders(m)” value defines the margins around the DEM
map.
11
Figure 6: Line of sight analysis tab before the “Calc” button was pressed. Antenna 6 was chosen and its details appear, and the marker
was set on an arbitrary position (blue =marker)
Upon clicking the “Calc” button, the line of sight between the two points
will be displayed in a separate plot (red line signify no LOS and green
signify clear LOS), and a map that includes the zoom-in DEM around the
line.
Figure 7: Line of sight analysis tab after the “Calc” button was pressed. A plot has been added to the screen that shows the altitude
below the chosen line. Another map was added that shows the detailed DEM in the vicinity of the line.
12
5. Obtaining DEM data
1 arc resolution (~ 30 meter) DEM files are available at
[Link] , a user can
download his region of interest ( as geographical anchored TIFF file) and then
provide it to the DEM <- getDEM(type="file",box,filename=DEM_name)
command. Some basic concatenating and saving commands might also be
useful:
1. Read a TIFF file
imported_raster1 <- raster(DEM_name1)
2. Concatenate two DEM files
imported_raster <-
merge(imported_raster1,imported_raster2)
3. save a merged raster file locally for further use:
writeRaster(imported_raster, filename="DEM_Harod.tif",
format="GTiff")
6. Earth curvature
Whenever taking into account earth curvature, we use the following approximation for its effect,
see Figure 8 for the simple vizualization
According to Pythagoras:
𝑅 2 + 𝑑 2 = (𝑅 + ℎ)^2
Simplifying we get
ℎ2 + 2ℎ𝑅 − 𝑑 2 = 0
To find ℎ:
𝑑2
−2𝑅 ± √4𝑅2 + 4𝑑 2 −2𝑅 ± 2𝑅√1 + 𝑅 2
ℎ= =
2 2
Using the Taylor expansion for square root, for 𝑥 ≪ 1 :
𝑥
√1 + 𝑥 ≈ 1 +
2
And the fact that Earth radius is far greater than the typical distance, and ignoring the negative
solution, we can simplify to
𝑑2 𝑑2
ℎ ≈ −𝑅 + 𝑅 (1 + ) ≈
2𝑅2 2𝑅
13
Within the code, whenever taking the curvature into consideration we use the value 𝑅 =
6371000𝑚, and subtract ℎ from the DEM value
Figure 8: Earth curvature effect on LOS, R signifies Earth radius, d, the distance between tower and target, and h, earth curvuture effect
at distance d
14