Last weekend I did some code tests related to path planning and trayectory control for UAVs. All the tests required working in different coordinate systems and access to geographic information systems.

In detail this post contains my notes on some Python libraries and tools that I found useful related to the WGS84 and UTM coordinate systems, the Digital Elevation Model (DEM), elevation profiles and the Open Source Geographic Information System QGIS.

**Universal Transverse Mercator coordinate system**

In our flight platform the vehicle's latitude and longitude belong to the WGS84 coordinate system although, at the level of waypoints calculation and precision, we will be interested on the Universal Transverse Mercator (UTM) conformal projection.

The UTM uses a 2-dimensional Cartesian coordinate system to give locations on the surface of the Earth independently of vertical position. Like the traditional method of latitude and longitude, it is a horizontal position representation.

The UTM system is not a single map projection. The system divides the Earth into sixty zones, each being a six-degree band of longitude, and uses a secant transverse Mercator projection in each zone.

We can use the 'utm' library as a simple and bidirectional UTM-WGS84 converter.

So, to convert a (latitude, longitude) tuple into an UTM coordinate:

>>> import utm >>> utm.from_latlon(43.38593855, -8.40656065) (548066.4735342383, 4803845.057732525, 29, 'T')

In the same way, to convert an UTM coordinate into a (latitude, longitude) tuple:

>>> import utm >>> utm.to_latlon(548066.4735342383, 4803845.057732525, 29, 'T') (43.3859385507919, -8.406560635013417)

We can verify that the library works as expected by calculating a new coordinate at a known distance (e.g. 1000 meters = 1 km) in the first coordinate system and verifying that the same distance is maintained in the second coordinate system. We can use the distance tool of Google Maps.

>>> import utm >>> (x, y, z, l) = utm.from_latlon(43.3710452, -8.396409) >>> (x, y, z, l) (548900.6904711355, 4802196.961608924, 29, 'T') >>> utm.to_latlon(x, y+1000, z, l) (43.380049197415424, -8.39631967000107)

The equations that solve the conversion of WGS84 to UTM and vice versa published by Coticchia & Surace, can be found in this resource.

The original source from which the equations are extracted seems to be the following document located in the library of Florence, Italy:

Coticchia, Alberto, Surace, Luciano. Monografia: Risoluzione di problemi geodetici con le minicalcolatrici elettroniche programmabili, Bollettino di geodesia e scienze affini, Anno 37. No. 1. P. 108 - 136. Editoriale: Istituto geografico militare, Firenze Italia. 1978. BIBLIOTECA NAZIONALE CENTRALE DI FIRENZE.

**Digital Elevation Model**

As mentioned in the previous point, the UTM does not consider the vertical position. To solve this drawback we can use a matrix of terrain elevation values, i.e., a Digital Elevation Model (DEM)

The elevation matrix can be indexed at the coordinate level taking into account that the vehicle altitude is relative to mean sea-level (MSL) in the case of LocationGlobal(), and the home position in the case of LocationGlobalRelative()

In our case we will use a Python script called Elevation to download global terrain digital elevation models, SRTM 30m DEM and SRTM 90m DEM.

As commented in the documentation, Elevation provides easy download, cache and access of the global datasets SRTM 30m Global 1 arc second V003 elaborated by NASA and NGA hosted on Amazon S3 and SRTM 90m Digital Elevation Database v4.1 elaborated by CGIAR-CSI.

Elevation will be very useful since the user experience covers both the command line interface and a simple API.

$ eio clip -o Rome-30m-DEM.tif --bounds 12.35 41.8 12.65 42 [...] $ ls Rome-30m-DEM.tif

At this point we can use QGIS to load and view the downloaded DEM.

It is also immediate to generate a contour line map for the area in QGIS.

It is interesting to note that if we use a version of QGIS higher than 3.0 we will have 3D view support, so we will be able to visualizate the 3D DEM easily.

In our case we will use the stable version 3.2.1 'Bonn'. The last releases are available here.

Take into account we can add also a world satellite imagery layer from ESRI web map server.

QGIS will automatically overlay the satellite imagery above the DEM data.

**Elevation profile**

Storing a matrix of elevation in memory and computing matrix operations can be a challenge in terms of temporal/spatial complexity and consumption for an UAV.

A lighter abstraction that we can easily derive from a DEM and use to program missions, paths and so on is an elevation profile. QGIS allows to obtain an elevation profile through a plugin called Profile tool.

Once activated the plugin, the tool plots profile lines from raster layers or point vector layer with elevation field. It supports multiple lines as well as graph export to svg, pdf, png or csv file.

The elevation profile indicates, starting from an origin point and reaching a destination point, the value of the corresponding height for each of the intermediate points.

In our example we draw a polyline with an arbitrary origin towards the center of the Roman Colosseum.

The plugin generates the elevation chart.

The plugin also allows to visually inspect the values of the profile and thus obtain the maximums and minimums easily.

The tool allows exporting the data for the values of the elevation chart. In the example, the units for the variables X and Y are meters.

0.0 29.0 54.648293915284334 25.0 109.29658783056867 24.0 ... 808.201193834435 25.0 1840.9052211536393 26.0 1840.9052211536393 26.0

There is also the option to include the coordinates (lon-lat) for the variable X and thus have the location, the distance between points not interpolated and the different elevations.

0.0 12.482835396791444 41.8881108248405 29.0 54.648293915284334 12.483117227466614 41.888411174560886 25.0 109.29658783056867 12.483399058141787 41.88871152428126 24.0 ... 808.201193834435 12.491999843922018 41.89025616876877 25.0 1840.9052211536393 12.49228365467219 41.89019945403496 26.0 1840.9052211536393 12.49228365467219 41.89019945403496 26.0

**Wrapping up**

Although it is practical and useful to operate in the domain of the UTM coordinates for the calculation of new waypoints, it is necessary to bear in mind that the conversion between coordinate systems is expensive.

The utm library is implemented in Python, and although it is notably faster than the more generic pyproj library, it works in the order of 0.4 - 0.5 seconds according to the tests documented in the project.

This conversion speed may be too high for missions with cruising speeds between 25 m/s and 35 m/s.

From an implementation point of view we can consider the rewriting of the library in a less heavy language as well as the use of a coordinate cache, or indexed and pre-computed coordinate tables.

It also seems logical to think that in a fleet of UAVs that can communicate with a common backend, this may be the point where these computations are carried out, either through an API or a set of higher level path planning services.

For fully distributed systems it is also interesting to think about how the computational load can be distributed among the different participating vehicles, as well as the most appropriate instants to carry it out.

In relation to the Data Elevation Model (DEM), the open source tools available to download and explore it are flexible and powerful.

The DEM can become a heavy structure to be transmitted from a remote backend by network, or be manipulated locally in the UAV. It seems more practical to work with a data structure that can be derived in turn from the DEM as is the elevation profile.

From a practical point of view, we are interested in the fact that an elevation profile is associated with a mission.

The elevation profile can be assigned to the mission once the coordinates of the waypoints are known and can be transmitted to the UAV at the same time as the coordinates.

As a first impression, it can be considered that the obtaining and manipulation of a DEM, although carried out in an automated way, seems more appropriate to perform it in a backend leaving the elevation profiles as lighter structures to be sent and manipulated by the UAVs.

In any case, the impact of the coordinate conversion algorithms and z-information must be considered in the design of the path planner of the solution.

- Open Source UAV, SITL in Docker, Mission Planner and MAV tools
- Open Source UAV and survival analysis with R
- Open Source UAV and mobile cellular networks
- Open Source UAV API, DroneKit-Python and Geopy
- Open Source UAV Autopilot with Ardupilot and Pixhawk