Sunday, 6 October 2013

Converting terrain data to a WebGL-friendly format

Three.js is a very promising tool if you want to add a third dimension to your web maps. In my last blog post, I showed how easy it was to create a WebGL Earth in the browser. Today, I'm starting a new blog series about terrain building with three.js.


Last year, I wrote a blog series about all the fun you can do with digital terrain data:
  1. Digital terrain modelling and mapping
  2. Creating hillshades with gdaldem
  3. Creating color relief and slope shading with gdaldem
  4. Terrain mapping with Mapnik
  5. Land cover mapping with Mapnik
  6. Using custom projections with TileCache, Mapnik and Leaflet
  7. Creating contour lines with GDAL and Mapnik
  8. Showing GPS tracks with Leaflet
I will continue using map data from Jotunheimen, a mountainous area of Norway. The 29 highest mountains in Norway are all in Jotunheimen, as well as the deepest valley, Utladalen. It's a great spot for 3D terrain mapping. But the same techniques applies to all terrains, - you only need some terrain data. Instead of Leaflet, we're going to use WebGL and three.js to render the maps. 

Last friday, 27th September 2013, was a milestone in the mapping history of Norway. The Norwegian Mapping Authority released its topographic datasets to the public, free of charge. Included was also a digital elevation model (DEM) of the Norwegian mainland, at 10 meters resolution. You can download the data on this page (unfortunately only in Norwegian), under a CC BY 3.0 licence. The terrain files created in this blog post are also available on GitHub.

Norway is divided into 50 x 50 km tiles, and you can select the areas you want to download on a map.
We need 4 tiles to cover Jotunheimen.  

The DEM files are using the USGS file format for terrain data, and I've selected the UTM 32N projection for my data. I'm using gdalbuildvrt to create a combined virtual dataset (jotunheimen.vrt) of the DEM files: 

gdalbuildvrt jotunheimen.vrt 6704_1_10m_z32.dem 6704_4_10m_z32.dem 6804_2_10m_z32.dem 6804_3_10m_z32.dem

Then we can use gdalwarp to clip the DEM to the area of interest and convert to GeoTIFF

gdalwarp -te 432000 6790000 492000 6850000 jotunheimen.vrt jotunheimen.tif

Use gdalinfo to see the properties of this image:

gdalinfo -mm jotunheimen.tif

Driver: GTiff/GeoTIFF
Files: jotunheimen.tif
Size is 6000, 6000
Coordinate System is:
PROJCS["WGS 84 / UTM zone 32N",
    GEOGCS["WGS 84",
        DATUM["WGS_1984",
            SPHEROID["WGS 84",6378137,298.257223563,
                AUTHORITY["EPSG","7030"]],
            AUTHORITY["EPSG","6326"]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433],
        AUTHORITY["EPSG","4326"]],
    PROJECTION["Transverse_Mercator"],
    PARAMETER["latitude_of_origin",0],
    PARAMETER["central_meridian",9],
    PARAMETER["scale_factor",0.9996],
    PARAMETER["false_easting",500000],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]],
    AUTHORITY["EPSG","32632"]]
Origin = (432000.000000000000000,6850000.000000000000000)
Pixel Size = (10.000000000000000,-10.000000000000000)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (  432000.000, 6850000.000) (  7d42'39.90"E, 61d46'36.81"N)
Lower Left  (  432000.000, 6790000.000) (  7d43'59.45"E, 61d14'18.21"N)
Upper Right (  492000.000, 6850000.000) (  8d50'54.02"E, 61d46'58.29"N)
Lower Right (  492000.000, 6790000.000) (  8d51' 3.39"E, 61d14'39.21"N)
Center      (  462000.000, 6820000.000) (  8d17' 9.19"E, 61d30'42.34"N)
Band 1 Block=6000x1 Type=Float32, ColorInterp=Gray
    Computed Min/Max=1.900,2467.800


What does it tell us?
  • The image or grid size is 6000 x 6000 pixels or cells.
  • The coordinate system is UTM zone 32N using the WGS 84 spheroid.
  • The coordinate unit is meters.
  • Each pixel or grid cell covers an area of 10 x 10 meters.
  • The extent of the area is defined by four coordinate pairs, covering an area of 3,600 km² (60 x 60 km).
  • The image is single band and the elevation data is stored as 16-bit integers. The minimum elevation is 1.9 meters and the maximum is 2,467.8 meters.

We got the data and the information we need to start terrain modelling. Our first task is to transfer the terrain data to the browser. As most web browsers don't support TIFF-files, it's common to convert the DEM into a heightmap. It can be thought of as a grayscale image where the intensity of each pixel represents the height at that position. Black indicates the minimum height and white the maximum height. It's easy to create a heightmap with gdal_translate:

gdal_translate -scale 0 2470 0 255 -outsize 200 200 -of PNG jotunheimen.tif jotunheimen.png

This command will create a PNG where the height values are reduced to 256 shades of gray. I'm also reducing the size from 6000 x 6000 px to only 200 x 200 px to save our GPU. Each pixel or elevation value is now covering an area of 300 meters.  


We can then transfer the image to the browser, and read the height values directly from the image. This can be sufficient for many uses, but I'm afraid my terrain will look "blocky" with only 256 height values. It's possible to use colour channels to get a wider height span with PNGs, but I couldn't find an easy way to achieve this with GDAL. By reading these notes, I saw that I could use a format called ENVI. In this format the height values are provided as a sequence of raw bytes, and I'm storing the values as 16-bit unsigned integers:

gdal_translate -scale 0 2470 0 65535 -ot UInt16 -outsize 200 200 -of ENVI jotunheimen.tif jotunheimen.bin

To preserve as much detail as possible, I'm scaling the floating point elevation values to the full range of a 16-bit unsigned integer (0-65535). I'm also creating a heightmap in full 10 m resolution of the Besseggen mountain ridge (2 x 2 km): 

gdalwarp -te 484500 6818000 486500 6820000 jotunheimen.vrt besseggen.tif

gdal_translate -scale 982 1742 0 255 -of PNG besseggen.tif besseggen.png

gdal_translate -scale 982 1905 0 65535 -ot UInt16 -of ENVI besseggen.tif besseggen.bin



Lastly, I'm creating a tiny heightmap of only 10 x 10 px for testing purposes: 

gdal_translate -scale 982 1905 0 255 -outsize 10 10 -of PNG besseggen.tif besseggen10.png

gdal_translate -scale 982 1905 0 65535 -outsize 10 10 -ot UInt16 -of ENVI besseggen.tif besseggen10.bin




In the next blog post, we'll start playing with our terrain data in three.js.

Autumn in Nordmarka in Oslo. Photo: Bjørn Sandvik, 5th October 2013. 

7 comments:

  1. For the gdalwarp, where do your values come from for: 432000 6790000 492000 6850000?

    ReplyDelete
  2. Are you familiar with any technique for smoothly transistioning from a two dimensional view to a three dimensional view?

    Consider a 2D square for example. If I push upwards on the lower most edge of the square I would like the image to transistion into a cube.

    If there is a 2D square and a 2D rectangle on a canvas if I push upwards on the lower edge of the canvas the image transistions into a 3D cube and 3D rectangle block.

    Google Earth has this similar effect. As you fly into the map the image smoothly transistions from a 2D map image to a 3D real world terrain image just as happens in real life. Any suggestions would be welcomed. Thanks.

    ReplyDelete
  3. Hello Bjorn,
    Love your work:) everything is explained in detail and final results are always astonishing.

    I am having a problem with converting my ASCII data to binary format as you described it. I have the data in WGS 84 format:
    50.7076,16.5009,640.00
    50.7073,16.4998,639.73

    and I try to use gdalbuildvrt I get the following:
    gdalbuildvrt cos.vrt M-33-45-D-d-4_p.asc
    0...10...20...30...40...50...60...70...80...90...100 - done.
    ERROR 1: Ungridded dataset: At line 3, X is 50.707300, where as 50.699700 was expected (1)
    Warning 1: Can't open M-33-45-D-d-4_p.asc. Skipping it

    I tried to use gdal_grid as they suggest on their page and I have received an empty .tiff file with following stats:

    $ gdalinfo p.tiff -stats
    Driver: GTiff/GeoTIFF
    Files: p.tiff
    p.tiff.aux.xml
    Size is 256, 256
    Coordinate System is `'
    Origin = (50.663200000000003,16.431699999999999)
    Pixel Size = (0.000188671875000,0.000284765625000)
    Image Structure Metadata:
    INTERLEAVE=BAND
    Corner Coordinates:
    Upper Left ( 50.6632000, 16.4317000)
    Lower Left ( 50.6632000, 16.5046000)
    Upper Right ( 50.7115000, 16.4317000)
    Lower Right ( 50.7115000, 16.5046000)
    Center ( 50.6873500, 16.4681500)
    Band 1 Block=256x4 Type=Float64, ColorInterp=Gray
    Min=474.683 Max=1008.174
    Minimum=474.683, Maximum=1008.174, Mean=720.570, StdDev=106.421
    Metadata:
    STATISTICS_MAXIMUM=1008.1740722656
    STATISTICS_MEAN=720.56987645384
    STATISTICS_MINIMUM=474.68276977539
    STATISTICS_STDDEV=106.42117154646

    What should I do to receive the correct tiff file (just as you did)?

    ReplyDelete
  4. Ok I have managed to sort that out by using gdal_grid to make my data gridded with the following command:

    $ gdal_grid -a_srs EPSG:4326 -of GTiff -ot Float64 -outsize 1000 1000 -l all all.vrt all.tif

    ReplyDelete
  5. I don`t understand were you get thee coordinates for the gdalwarp.

    ReplyDelete
  6. Thanks for this walk through the gdal tools, they’re not the most accessible. Quick question — how did you figure out that “the elevation data is stored as 16-bit integers” from the gdalinfo output? Based on this line, I’d have guessed they were stored as 32-bit floating point values:

    Band 1 Block=6000x1 Type=Float32, ColorInterp=Gray


    ReplyDelete