This week we launched a new API providing access to our high resolution Digital Surface Model (DSM) at 7.5cm! This data has been available to our customers over the past years via download of large GeoTiffs from our cloud storage, ideal for offline use. But now with the DSM API a whole world of applications opens up with on-demand access to the Surface Model on an individual property basis. With the surface model, you know the elevation in meters above the Earth’s surface of every pixel. Whether you simply need to know the height of an exhaust pipe on a building or you need to accurately measure the north facing facet of a roof, our DSM provides you with the necessary data. Detecting the footprint of a building, calculating the volume of the building, or identifying and measuring individual roof slopes and sizes are all straightforward operations with the our DSM.
And it couldn’t be any easier to use. The API returns standard 256X256 tiles in Geotiff format. As these tiles conform to the WMTS spec, they align perfectly with other Vexcel data like our True Orthomosaic tiles and Near-Infrared tiles, as well as tiled data from other providers in the industry. Here is an excellent overview on the standard tiling system from Microsoft, which also includes some helpful sample code for converting between lat/long and tile coordinates. If you know how to request an X/Y/Z tile from another provider, you already know how to request data from our DSM API. From street map tiles from Open Street Map to Land Parcels from Landgrid, our DSM aligns perfectly without any painful up/down sampling or transform of the data. and of course it aligns with our own high quality aerial imagery.
Here are two tiles at location Tile X = 359441, Tile Y = 836948, Zoom=21. The first tile is from the Vexcel Blue Sky Ultra Orthomosaic. The second image is a simple rendering of the SAME geography’s DSM. Note the clean alignment between them. the black areas are flat to the Earth’s Surface, with the brightest red areas representing the highest points. the Delta here is 11.87 meters from the ground to the flat roof. the elevation of EVERY pixel is easily read from the Geotiff.
Here is a series of tiles stitched together showing the orthomosaic and DSM together. use the slider to compare.
Getting started with the DSM API is easy. Begin with the documentation which you can access here. If all you need to do is display the DSM data in an application like QGIS that knows how to render DSM, simply download some tiles and drop them into a supported application. But if you want to access the elevation data programmatically in your own applications, you’ll likely want to use one of the many code libraries out there that making working with Tiff images (or the GeoTiff extension) quite easy. Here is a link to LibTiff, a very popular library for opening Tiff Files: http://libtiff.org/ And for .NET developers, here is a port of Libtiff for .NET: https://bitmiracle.com/libtiff/ LibTiff is so popular that many of the Geotiff libraries available rely on it under the hood. So you can’t go wrong by looking into LibTiff to speed your development process with our DSM.
Below is some sample code written in C# that shows how to open a DSM Geotiff tile and read the matrix of elevation values out of it. It relies on the aforementioned LibTiff libarary, but you should be able to adapt this code to any language and library of your choice. There are two important extended tags in the Geotiff that you will want to be aware of. Extended Tag 33550 contains the Pixelsize information. Tag 33922 contains the X/Y location on earth of the tile, in Meters. These tags are highlighted in the comments in the code below.
This should get you off to a good start working with the DSM API, but as always reach out to our support team with any questions. support@vexcelgroup.com In future articles here, we’ll build on the code below to do more interesting things like building footprint detection and 3D extrusion.
using System;
using BitMiracle.LibTiff.Classic;
using System.Net;
namespace DSMtileConsole
{
class Program
{
private const double EarthRadius = 6378137;
static void Main(string[] args)
{
string outputFolder = System.IO.Path.GetTempPath();
string VexcelToken = "YOURTOKENHERE";
//Request a DSM Geotiff tile from the service and save it as a local file
string filename = outputFolder + "TestDSMTile.tiff";
string imageURL_DSM = "https://api.gic.org/images/GetDSMTile/21/359441/836948/?layer=bluesky-ultra" +
"&token=" + VexcelToken;
fetchImageFromURL(imageURL_DSM, filename);
//Open and Parse the Geotiff file
using (Tiff tiff = Tiff.Open(filename, "r"))
{
if (tiff == null)
{
// Error - could not open
Console.WriteLine("ERROR: Couldn't open the file. ");
return;
}
Console.WriteLine("Input File: " + filename);
//Read the image Height and width from standard tags
int imageWidth = tiff.GetField(TiffTag.IMAGEWIDTH)[0].ToInt();
int imageHeight = tiff.GetField(TiffTag.IMAGELENGTH)[0].ToInt();
Console.WriteLine("Image width: " + imageWidth + " Image height: " + imageHeight);
//Parse the Pixel Width and Height. this information is contained in Extended Tag 33550
//These values are each 8 byte doubles. the first double is PixelsSizeX and the second is PixelsSizeY
byte[] destin = new byte[8];
FieldValue[] values = tiff.GetField((TiffTag)33550);
byte[] bytesGeo = values[1].ToByteArray();
Array.Copy(bytesGeo, 0, destin, 0, 8);
double PixelSizeX = BitConverter.ToDouble(destin, 0);
Array.Copy(bytesGeo, 8, destin, 0, 8);
double PixelSizeY = BitConverter.ToDouble(destin, 0);
Console.WriteLine("Pixel Size X: " + PixelSizeX + " Pixel Size Y: " + PixelSizeY);
//Parse the X and Y origin. this represents the coordinate of the upper left pixel of the geotiff
//these values are in meters and found in Extended Tag 33922
//These values are each 8 byte doubles. the 4th double is OffsetX and the fifth is OffsetY
values = tiff.GetField((TiffTag)33922);
bytesGeo = values[1].ToByteArray();
Array.Copy(bytesGeo, 24, destin, 0, 8);
double OriginMetersX = BitConverter.ToDouble(destin, 0);
Array.Copy(bytesGeo, 32, destin, 0, 8);
double OriginMetersY = BitConverter.ToDouble(destin, 0);
Console.WriteLine("X Origin: " + OriginMetersX + " Y Origin: " + OriginMetersY);
//Convert the X,Y from meters to decimal degrees
MetersXYToLatLong(OriginMetersX, OriginMetersY, out double originLatitude, out double originLongitude);
Console.WriteLine("Origin Longitude: " + originLongitude + " Origin Latitude: " + originLatitude);
//Now we can start reading the body of the geotiff containing the elevation data
//the body is organized as a series of 4 byte Singles
Console.WriteLine("Reading the body of the file...");
int stripSize = tiff.StripSize();
int stripMax = tiff.NumberOfStrips();
int bufferSize = stripMax * stripSize;
byte[] buf = new byte[bufferSize];
int result = tiff.ReadEncodedTile(0, buf, 0, stripSize);
//Scan the file. Report each pixel's Elevation and find the min/max elevation along the way
Single minElevation = 999999;
Single maxElevation = 0;
Single curEvelation = 0;
FieldValue[] value = tiff.GetField(TiffTag.IMAGELENGTH);
int imageLength = value[0].ToInt();
string outline = "";
for (int row = 0; row < imageLength; row++)
{
outline = "Row " + row.ToString().PadRight(4,' ');
for (int col = 0; col < (buf.Length / (imageWidth * 4)); col++)
{
int offset = (row * (imageWidth * 4)) + (col * 4);
byte[] bytes = { buf[offset], buf[offset + 1], buf[offset + 2], buf[offset + 3] };
curEvelation = BitConverter.ToSingle(bytes, 0);
minElevation = Math.Min(minElevation, curEvelation);
maxElevation = Math.Max(maxElevation, curEvelation);
outline += Truncate( curEvelation.ToString(),6).PadLeft(8, ' ');
}
//Console.WriteLine(outline);
}
Single elevationDelta = maxElevation - minElevation;
Console.WriteLine("Elevation range: " + minElevation + " to " + maxElevation);
Console.WriteLine("Elevation delta: " + elevationDelta );
}
}
public static void fetchImageFromURL(string imageURL, string outputFilePath)
{
using (WebClient webClient = new WebClient()) {
try {
webClient.DownloadFile(imageURL, outputFilePath);
}
catch (Exception ex)
{
}
}
return;
}
public static void MetersXYToLatLong(double x, double y, out double latitude, out double longitude)
{
//Adapted from https://docs.microsoft.com/en-us/bingmaps/articles/bing-maps-tile-system Thanks Joe!!
latitude = 90 - 360 * Math.Atan(Math.Exp(-y / EarthRadius)) / Math.PI;
longitude = 360 * x / (EarthRadius * 2 * Math.PI);
}
public static string Truncate(string source, int length)
{
if (source.Length > length)
{
source = source.Substring(0, length);
}
return source;
}
}
}
The Ground Truth 