Working with GIS, terrains and

Or, how to build a precise 3D terrain from any place of the world.

Again not much visually significant FreeCAD development to show this week, so here is another interesting subject, that I started looking at in an earlier post.

We architects should really begin to learn about GIS.

GIS stands for Geographic information system and begins to be more and more used by administrations around the world, specially cities, to manage their physical area and infrastucture in all their complexities.

Also, it is becoming more and more their preferred way to deliver data that you need to use when doing a project on a certain plot of terrain in that city. Before, you would go to the city council, and ask for the data they have on a specific plot. In the old times, they would give you a photocopy of a piece of some plan, and a form filled with some additional information. In the recent years, most of them had switched to CAD systems and gave you DWG files to work with.

Now, DWG files are being abandoned worldwide, and, together with the better integration of the different services that compose a city council, they are adopting GIS systems. So we'd better learn how it works.

It is actually far better than before. GIS is a bit like CAD, but adapted for working with terrains, maps and geographic data. GIS files are usually made of layers, like a CAD file, but layers can contain vectorial data, like CAD, or bitmap (that they call raster) images. The most important feature, though, is that everything is georeferenced, which means that any point has a precise, exact position on earth. So you can join all kinds of data coming from all kinds of sources in a same file, everything will stack exactly on top of each other. In GIS you never move anything.

The amount of data that you can find online, already properly formatted, that you can simply drag and drop in your GIS application, to build beautiful maps with a huge quantity of layers is really impressive. In this article, I'll use data from my hometown, Brussels, and the city I live in now, São Paulo. Both have a public GIS website from where you can download data. I'll also take elevation maps directly from the NASA website. I'll use the most well-known, open-source, multiplatform GIS app available, called QGIS. I recommend you to install it, I bet very soon a GIS tool will be part of any architect's toolbox. Fianlyl, we'll export our map to a DXF file so we can build architecture projects on it, and build a 3D terrain model in FreeCAD. After installing, check the "Plugins" menu of QGIS, it has an amazing collection of useful plugins, a favorite being the openlayer plugin, which lets you work on top of an OpenStretMap layer, which is an excellent way to start a new work in QGIS.

Obtaining GIS data from the São Paulo city council is easy. Head to their GIS portal, zoom to the neighborhood of your choice, click the "Download" button in the left toolbar, and select the type of data you wish to download. Be sure to select "shapefiles" data type, which is the most common GIS format.

Brussels also has a similar geodata portal, a little bit more complex to navigate, but allow to download the same kind of data. However, they also have a much simpler download app which lets you easily specify the dataset you want to download and get the same shapefiles.

There is also one precious type of data available on the NASA website: Elevation maps. Several times in the history, they have covered the world with a satellite able to find map the elevation of any point on earth (apparently with very simple techniques like sending a wave signal and measuring the time it takes to bounce back). The latest data, collected with the ASTER satellite in 2009, has a 30x30m resolution of the whole earth, and is available freely on the NASA Echo Reverb website. The data comes in a special kind of bitmap image format called GeoTIFF which is basically a greyscale TIFF image with floating-point values (unlimited number of grey shades) and it is, of course, georeferenced.

Downloading data from there is however a bit more complex:

  1. egister a new account (it looks commercial, but ASTER data is free)
  2. Log in
  3. Zoom on the map
  4. Draw a rectangle
  5. In the list below, mark "ASTER Global Digital Elevation Model V002"
  6. Click Search for Granules
  7. In the next screen, under Granules you will get one (or more) zip files like Add them to the cart
  8. Click view items in cart. You will get all the zips you added, plus a "ASTER Global Digital Elevation Model V002" that you don't need to worry about
  9. Click "order" (it definitely doesn't look like it, but it's free)
  10. Accept all the requirements and click "submit order"
  11. You will get an email with a download link. Normally it's immediate, but one time I had to wait an hour or two.
  12. Unzip, and unzip another zip inside, there is your geotiff

Once we have all this data unzipped, we can just drag-and-drop the geotiffs and the shapefiles from the Browser Panel of QGIS into a fresh, blank document. Right-clicking on each new layer in the Layers Panel will allow to change visual properties such as color and transparency.

One very important note: Projecting coordinates onto the earth is a very complex business (think that the earth is not even completely spherical!). As a result, there are hundreds of different projection methods. You will quickly see that a shapefile is actually several files. One of these files defines the projection method. On the GIS portals above they usually inform which projection system they use for their data. In this case, Brussels uses Belgian Lambert 1972 projection, and São Paulo uses SAD69/96 UTM 23S. So make sure, when you insert a new layer, that it uses the correct projection. You can check or change by right-clicking a layer and "Set Layer CRS".

Above is Brussels central area buildings layer with the corresponding GeoTIFF image in the background. Generating contour lines from a GeoTIFF image involves two steps: Cropping (the area of the GeoTIFF is usually way too large for what we need) and generating the contour lines. Both operations are available from menu Raster -> Extraction. When using the Contour tool, be sure to turn "Attribute generation" on. This will, together with the contour lines, build a table and store the elevation of each line. This will be useful when we'll recreate a 3D model.

Our map is now complete, and if we right-click the contours layer, and check the attributes, we see our elevations list:

Before exporting to DXF, we need to perform one last step. The DXF exporter of QGIS will only export vector layers with the same projection system. However, the contours have taken the same projection system as the GeoTIFF, which are not Belgian Lambert. We then need to convert first. This is done by saving the contours layer (right-click it -> Save as) with the Belgian Lambert 72 projection system. Once that is done, we can discard the old one, and save to DXF via File -> DXF export. Make sure you mark the "Export features intersecting the current map extent" option, otherwise you'll be exporting the whole city, which will be quite huge.

We can then open our file in our vavorite DXF application(such as QCAD and do the necessary cleaning, remove stuff we don't need, etc. Both DXF files exported with Lambert or SAD69 will produce files in meters, I don't know if it is the case with any projection system, you'd better verify if using another system.

Unfortunately this DXF contains only 2D data. If you check the "Elevation" property of the terrain polylines, they are all zero. However, GDAL, the processing engine behind QGIS, is able to produce such a 3D DXF file from a contour shapefile with an "ELEV" attribute as we've taken care to do. Open a shell window in the folder that contains your contours file, and run a command like:

ogr2ogr -f "dxf" contours-lambert.dxf contours-lambert.shp -zfield ELEV

Changing "contours-lambert" by the name you used to save the contour file. You will see that a DXF file has been added. I wish someone would produce a plugin to do this from inside QGIS!

The new DXF file has the exact same coordinates as the previous one, so you can just import it in the former one, everything should just click into place (Any problem, open it in a separate tab and use the "Copy with basepoint" function of your CAD app and copy/paste the curves from/to the same point such as 0,0). Make sure to leave the 3D curves in their separate layer, so we can isolate them later to build the terrain.

It's now time to open the DXF file in FreeCAD. Make sure the "group objects by layer" option is turned on in the DXF preferences so we get all our curves as one object:

To build a nice BSpline surface from our curves can be easily done from the python console (at the moment there is no GUI tool for it). First we need to turn our curved into a points cloud by activating the "Points" workbench and using the "Points -> Convert to Points..." menu item, with a rather large value, like 5 or 10. After that, issuing these commands in the python console will create the BSpline surface (change the getObject("Points") with yours if your points cloud object is not named "Points"):

import Part,ReverseEngineeringobj = FreeCAD.ActiveDocument.getObject("Points"),NbUPoles=16,NbVPoles=16).toShape())

This might take a little while. Use lower values than 16 to get a coarser but faster result, or higher to get a finer result. But you'll get a nice 3D surface built from the curves. This surface can then be measured, cropped, turned into a solid, etc.. with precise calculations.