CALCULATE LINE OF SIGHT FROM POINT TO POINT USING GDAL ON LEAFLET This blog illustrates the methods and procedures involved in the calculation of Line of Sight (LOS). Line of sight helps in the calculation of intervisibility between two points and is defined as the path followed by a signal from an emitter at position A to reach an observer at position B which are a certain distance apart. If this line is blocked by the geographical relief/terrain, then intervisibility is lost. The following graphs show line of sight with both possibilities. Line of sight with intervisibility Line of sight with no intervisibility For this calculation, we need DEM files which are available online for the whole world. This files have the AMSL values for every x meters depending on its spatial resolution. How to retrieve these values on a web page has already been explained in the previous blog on displaying altitude on map. In this blog, I will be using Java language as scripting language and GDAL java library for determining the status of LOS between two given points. If you do not know Java language, you can go through Guru99 for an overview of Java flavour. Following pre-requisites are to be fulfilled before proceeding to the coding part. Set up the GDAL environment on you PC/Server depending on its architecture (x64/x86). Download the ASTER GDEM/SRTM DEM files of your area of interest (AOI) and generate its VRT file as explained in the previous blog on displaying altitude on map. We require an IDE for coding in Java. The most popular are Netbeans and Eclipse. We require Geographiclib-Java-1.47.jar for splitting line of sight apart from gdal.jar. These files would be set in Java built path library location of the project later. LINE OF SIGHT CALCULATION METHODOLOGY Consider two points on earth as A (emitter) and B (observer) which are x meter apart. We need to find the intervibility of these two points. We would draw a line connecting these two points as shown below. This line is called the Line of Sight. This Line of Sight is split into n equidistant number of points depending on the spatial resolution of DEM or user requirement (using GeographicLib-Java-1.47.jar) and at each point the value of elevation is determined using the same process explain in the previous blog on displaying altitude on map. Simulteneously, the height of the line of sight is calculated using the simple mathematics. Both these elevation and height of LOS is compared and at any point if the value of elevation is more than the height of LOS, then A and B are said to be invisible or else they achieve intervisibility. The right side diagrams explains it. Let us try to achieve this concept through GDAL and Geographiclib in Java. Open a new Web application project in your IDE say losCalc. Create a servlet say CalcLOS.java in your xyz package. We will be using only the doget method here to handle the incoming request for LOS status between asked locations. We will create a Java Class file say GDALLOS.java which will evaluate the line of sight status i.e, visible or invisible. You need to add gdal.jar and GeographicLib-Java-1.47.jar files in the JAVA BUILD PATH Library location. For Netbeans Right-Click project --> Properties --> Libraries --> Compile tab For Eclipse Right-Click project --> Properties --> Java Build Path --> Libraries tab. Create a java class called Coordinate.java and create class variables like latitude, longitude and altitude along with constructors and setters/getters as follows: ```package xyz; public class Coordinate { private double lat; private double lon; private double height = 0; //initialize to 0; public Coordinate(){ } public Coordinate(double height){ this.height = height; } public Coordinate(int height){ this.height = height; } public Coordinate(double lat, double lon){ this.lat = lat; this.lon = lon; } public Coordinate(double lat, double lon, double height){ this.lat = lat; this.lon = lon; this.height = height; } public double getLat() { return lat; } public void setLat(double lat) { this.lat = lat; } public double getLon() { return lon; } public void setLon(double lon) { this.lon = lon; } public double getHeight() { return height; } public void setHeight(double height) { this.height = height; } } ``` Copy the following Code in your GDALLOS.java file and create setters/getters as shown below: ```package xyz; import java.util.ArrayList; import java.util.Iterator; import net.sf.geographiclib.*; import org.gdal.gdal.*; import org.gdal.gdalconst.gdalconstConstants; public class GDALLOS { private Coordinate coordTransmitter; // Transmitter Location private Coordinate coordObs; // Observer Location private int obsHeight = 8000; // Default altitude private double transOffset = 0.0d; // Transmitter height above ground private int accuracyDem = 30; // DEM resolution protected final static String GDEM_location = "D:\\GDEMv2\\TIF\\GDemV2.vrt"; // DEM VRT location private static Dataset hDataset = null; // DEM GDAL dataset initialised to null private final static Geodesic geod = Geodesic.WGS84;//Default datum private static void regDataset(){ //method to initialise static dataset only once. if(hDataset==null){ gdal.AllRegister(); //Registers GDAL hDataset = gdal.Open(GDEM_location, gdalconstConstants.GA_ReadOnly); } } public Viewshed(){ // Default constructor regDataset(); } public Coordinate getCoordTransmitter() { return coordTransmitter; } public void setCoordTransmitter(Coordinate coordTransmitter) { this.coordTransmitter = coordTransmitter; } public Coordinate getCoordObs() { return coordObs; } public void setCoordObs(Coordinate coordObs) { this.coordObs = coordObs; } public int getObsHeight() { return obsHeight; } public void setObsHeight(int obsHeight) { this.obsHeight = obsHeight; } public double getTransOffset() { return transOffset; } public void setTransOffset(double transOffset) { this.transOffset = transOffset; } public int getAccuracyDem() { return accuracyDem; } public void setAccuracyDem(int accuracyDem) { this.accuracyDem = accuracyDem; } } ``` The following methods are added to the above GDALLOS.java files. splitCoord: This will split the line joining Transmitter Coordinate and Observer Coordinate into n number of equidistant points and save them in an ArrayList. ```public ArrayList splitCoord(){ ArrayList cd = new ArrayList(); GeodesicLine line = geod.InverseLine(coordTgt.getLat(), coordTgt.getLon(), coordObs.getLat(), coordObs.getLon(), GeodesicMask.DISTANCE_IN | GeodesicMask.LATITUDE | GeodesicMask.LONGITUDE); int num = (int)(Math.ceil(line.Distance() / accuracyDem)); { // Use intervals of equal length double ds = line.Distance() / num; for (int i = 0; i <= num; ++i) { GeodesicData g1 = line.Position(i * ds, GeodesicMask.LATITUDE | GeodesicMask.LONGITUDE); Coordinate c = new Coordinate(g1.lat2, g1.lon2); cd.add(c); } } return cd; } ``` splitCoord(Coordinate tran, Coordinate obs): This method splits the provided Transmitter Coordinate and Observer Coordinate into n number of equidistant points and save them in an ArrayList. In java, this method is an example of polymorphism in Java. ```public ArrayList splitCoord(Coordinate tran, Coordinate obs){ ArrayList cd = new ArrayList(); GeodesicLine line = geod.InverseLine(tran.getLat(), tran.getLon(), obs.getLat(), obs.getLon(), GeodesicMask.DISTANCE_IN | GeodesicMask.LATITUDE | GeodesicMask.LONGITUDE); int num = (int)(Math.ceil(line.Distance() / accuracyDem)); { // Use intervals of equal length double ds = line.Distance() / num; for (int i = 0; i <= num; ++i) { GeodesicData g1 = line.Position(i * ds, GeodesicMask.LATITUDE | GeodesicMask.LONGITUDE); Coordinate c = new Coordinate(g1.lat2, g1.lon2); cd.add(c); } } return cd; } ``` getElevation: This method calculates the value of ground elevation at the provided coordinate using the same technique explained in the previous blog. ```private static String getEle(Coordinate A){ regDataset(); Band hBand; int iBand; double[] adfGeoTransform = new double[6]; int iPixel; int iLine; double[] value = new double[1]; String elevationval = null; int X = hDataset.getRasterXSize(); int Y = hDataset.getRasterYSize(); if (hDataset == null) { return "NODATA"; } iBand = hDataset.getRasterCount(); if( iBand == 1 ){ hBand = hDataset.GetRasterBand(iBand); hDataset.GetGeoTransform(adfGeoTransform); if (adfGeoTransform[2]==0 && adfGeoTransform[4]==0){ iPixel = (int) Math.floor((A.getLon() - adfGeoTransform[0]) / adfGeoTransform[1]); iLine = (int) Math.floor((A.getLat() - adfGeoTransform[3]) / adfGeoTransform[5] ); if (iPixel<0 || iLine<0 || iPixel>X||iLine>Y){ return "0"; }else{ try{ hBand.ReadRaster(iPixel, iLine, 1, 1, value); elevationval = String.valueOf((int)value[0]); hBand.delete(); }catch(Exception e){ elevationval = "0"; } } }else{elevationval = "INVALIDDATA";} }else{elevationval = "MULTIBANDDEM";} return elevationval ; } ``` getGHeight(Coordinate A): This method calculates the value of ground elevation using the above method at the provided coordinate and returns a double (datatype) value. ```private double getGHeight(Coordinate A){ double height = 0; //default Elevation try{ height = double.parseDouble(getElevation(A)); }catch(NumberFormatException e){} return height; } ``` getLHeight(Coordinate tran, Coordinate obs): This method calculates the LOS height at any location between Transmitter Coordinate and Observer Coordinate. Notice the class variable coordObs and transOffset below. ```private double getLHeight(Coordinate tran, Coordinate obs){ double effG = getGHeight(tran); double effHeight = ((obsHeight+transOffset-effG)/getDistance(tran, coordObs))*getDistance(tran, obs) + effG; return effHeight; } ``` ViewStatus(): This method finally check the visibility status between both the cordinates and returns 1 for true visiblility and 0 for no visibility. ```public boolean viewStatus(){ boolean stat = true; //default value true means transmitter is visible ArrayList c = splitCoord(coordTgt, coordObs); //Split for (Coordinate t : c) { double gHeight = getGHeight(t); double losHeight = getLHeight(coordTransmitter, t); if((gHeight - losHeight)>0){ //if elevation is greater than los height? stat = false; break; } } return stat; } ``` If at any point the value of ground height (AMSL) is more than the line of sight height, the intervisibility is said to be false and the value of stat becomes 0. Finally create a servlet say CalcLOS.java to handle the request from leaflet map in doget method to check if intervisibility is achieved between two points. ```protected void doGet(HttpServletRequest request, HttpServletResponse response) throws ServletException, IOException { processRequest(request, response); PrintWriter out = response.getWriter(); double latA = 0; double lonA = 0; double latB = 0; double lonB = 0; int obsHeight = 8000; //default double transHeight = 0.0d; //default Coordinate tgt = new Coordinate(); Coordinate obs = new Coordinate(); String lata = request.getParameter("lata"); String lona = request.getParameter("lona"); String latb = request.getParameter("latb"); String lonb = request.getParameter("lonb"); String obsheight = request.getParameter("obsHeight"); String transheight = request.getParameter("transHeight"); try{ latA = Double.parseDouble(lata); lonA = Double.parseDouble(lona); latB = Double.parseDouble(latb); lonB = Double.parseDouble(lonb); if(obsheight!=null){ obsHeight = Integer.parseInt(obsheight); } if(transheight!=null){ transHeight = Double.parseDouble(transheight); } tgt.setLat(latA); tgt.setLon(lonA); obs.setLat(latB); obs.setLon(lonB); GDALLOS vs = new GDALLOS(); vs.setCoordTgt(tgt); vs.setCoordObs(obs); vs.setObsHeight(obsHeight); vs.setTransOffset(TransHeight); if(vs.ViewStatus()){ out.println("1"); }else{ out.println("0"); } }catch(Exception e){ out.println(e.toString()); } } ``` Run this web application in Java supported server like Tomcat, Glassfish, JBoss etc., and pass the values of latitude and longitude as parameter in url. For example if your server runs on port 8080 on http on localhost then the url will be `http://localhost:8080/losCalc/CalcLOS?lata=23.5687&lona=45.9862&latb=23.3258&lonb=45.7862&obsHeight=9000&transHeight=25` This will print 1 if the intervisibility is achieved and 0 if the transmitter is obscured by the geographical feature. In case of any error/exception, it will be printed on screen. You call pass parameters through Ajax calls using EventListener either click or mousemove on your GIS web clients. Remember to check for cross-origin CORS if web client and the above codes run on different servers. Remember that while working with java code along with GDAL on web aplications, you need to restart the server every time you make some changes with Java codes else it would show UnsatisfiedLinkError from GDAL. The above mentioned procedure for calculation of LOS will give good accuracy for upto 20-30 km only. If the requirement exists for LOS calculation beyond 30 Km then the effect of earth curvature and atmospheric refrection need to be accounted for.
 QUICK LINKS Coordinates Earth Shape Ellipsoid & Geoid Datum Projection UTM Rasters Vectors GDAL GDAL Setup GDAL Rasters GDAL Vectors GDAL Tools WMS MS4W Spatial DB Oracle Spatial SDO_GEOMETRY Geo-Crawler Web Clients OL3+ Closest Distance Map Openlayers v2.x Openlayers+Cesiumjs ThreeJS Terrain Viewer Bluemarble (250m tiff) JS Tools Geographiclib JS2Shapefile CSV2Shape