Page 1 of 1 (5 posts)

  • talks about »
  • dem


Last update:
Sun Sep 14 21:00:08 2014

A Django site.

QGIS Planet

Selective data removal in an elevation map by means of floodfilling

Do you also sometimes get maps which contain zero (0) rather than NULL (no data) in some parts of the map? This can be easily solved with “floodfilling”, even in a GIS.

My original map looks like this (here, Trentino elevation model):

The light blue parts should be no data (NULL) rather than zero (0)...

Now what? In a paint software we would simply use bucket fill but what about GIS data? Well, we can do something similar using “clumping”. It requires a bit of computational time but works perfectly, even for large DEMs, e.g., all Italy at 20m resolution. Using the open source software GRASS GIS 7, we can compute all “clumps” (that are many for a floating point DEM!):

# first we set the computational region to the raster map:
g.region rast=pat_DTM_2008_derived_2m -p
r.clump pat_DTM_2008_derived_2m out=pat_DTM_2008_derived_2m_clump

The resulting clump map produced by r.clump is nicely colorized:

Clumped map derived from DEM (generated with r.clump)

As we can see, the area of interest (province) is now surrounded by three clumps. With a simple map algebra statement (r.mapcalc or GUI calculator) we can create a MASK by assigning these outer boundary clumps to NULL and the other “good” clumps to 1:

r.mapcalc "no_data_mask = if(pat_DTM_2008_derived_2m_clump == 264485050 || \
   pat_DTM_2008_derived_2m_clump == 197926480 || \
   pat_DTM_2008_derived_2m_clump == 3, null(), 1)"

This mask map looks like this:

Mask map from all clumps except for the large outer clumps

We now activate this MASK and generate a copy of the original map into a new map name by using map algebra again (this just keeps the data matched by the MASK). Eventually we remove the MASK and verify the result:

# apply the mask
r.mask no_data_mask
# generate a copy of the DEM, filter on the fly
r.mapcalc "pat_DTM_2008_derived_2m_fixed = pat_DTM_2008_derived_2m"
# assign a nice color table
r.colors pat_DTM_2008_derived_2m_fixed color=srtmplus
# remove the MASK
r.mask -r

And the final DEM is now properly cleaned up in terms of NULL values (no data):

DEM cleaned up for no data


The post Selective data removal in an elevation map by means of floodfilling appeared first on GFOSS Blog | GRASS GIS Courses.

Will the sun shine on us?

Recently, in order to nicely plan ahead for a birthday lunch at the agritur Malga Brigolina (a farm-restaurant near Sopramonte di Trento, Italy, at 1,000m a.s.l. in the Southern Alps), friends of mine asked me the day before:

Will the place be sunny at lunch time for a nice walk?

[well, the weather was close to clear sky conditions but mountains are high here and casting long shadows in the winter time].


A rather easy task I thought, so I got my tools ready since that was an occasion to verify the predictions with some photos! Thanks to the new EU-DEM at 25m I was able to perform the computations right away in a metric system rather than dealing with degree in LatLong.

Direct sunlight can be assessed from the beam radiation map of GRASS GIS’ r.sun when running it for a specific day and time. But even easier, there is a new Addon for GRASS GIS 7 which calculates right away time series of insolation maps given start/stop timestamps and a time step: r.sun.hourly. This shrinks the overall effort to almost nothing.

1. Creating the direct sunlight maps

The first step is to calculate where direct sunlight reaches the ground. Here the input elevation map is the European “eu_dem_25″, while the output is the beam radiation for a certain day (15 Dec 2013 is DOY 349).

Important hint: the computational region must be large enough to east/south/west to capure the cast shadow effects of all relevant surrounding mountains.

I let calculations start at 8am and finish at 5pm which an hourly time step. The authors have kindly parallelized r.sun.hourly, so I let it run on four processors simultaneously to speed up:

# calculate DOY (day-of-year) from given date:
date -d 2013-12-15 +%j

# calculate beam radiation maps for a given time period
# (note: minutes are to be given in decimal time: 30min = 0.5)
r.sun.hourly elev_in=eu_dem_25 beam_rad=trento_beam_doy349 \
      start_time=8 end_time=17 time_step=1 day=349 year=2013 \

The ten resulting maps contain the beam radiation for each pixel considering the cast shadow effects of the pixel-surrounding mountains. However, the question was not to calculate irradiance raster maps in Wh/m^2 but simply “sun-yes” or “sun-no”. So a subsequent filtering had to be applied. Each beam radiation map was filtered: if pixel value equal to 0 then “sun-no”, otherwise “sun-yes” (what my friends wanted to achieve; effectively a conversion into a binary map).
Best done in a simple shell script loop:

[Edit 30 Dec 2013: thanks to Anna you can simplify below loop to the r.colors call thanks to the new -b flag for binary output in r.sun.hourly!]

for map in `g.mlist rast pattern="trento_beam_doy349*"` ; do
    # rename current map to tmp
    g.rename rast=$map,$map.tmp
    # filter and save with original name
    r.mapcalc "$map = if($map.tmp == 0, null(), 1)"
    # colorize the binary map
    echo "1 yellow" | r.colors $map rules=-
    # remove cruft
    g.remove rast=$map.tmp

As a result we got ten binary maps, ideal for using them as overlay with shaded DEMs or OpenStreetMap layers. The areas exposed to direct sunlight are shown in yellow.

Trento direct sunlight 15 Dec 2013 Animation

Trento, direct sunlight, 15 Dec 2013 between 10am and 5pm (See here for creating an animated GIF). Quality reduced for this blog.

2. Time to look at some details: So, is Malga Brigolina in sunlight at lunch time?

Situation at 12pm (noon): predicted that the restaurant is still in shadow – confirmed in the photo:


(click to enlarge)

Situation at 1:30/2:00pm: sun is getting closer to the Malga, as confirmed in photo (note that the left photo is 20min ahead of the map). The small street in the right photo is still in the shadow as predicted in the map):

malga_brigolina_direct_sunlight_15dec2013_14pm_foto(click to enlarge)

Situation at 3:00pm: sun here and there at Malga Brigolina:

malga_brigolina_direct_sunlight_15dec2013_15pm_3DSunlight map blended with OpenStreetmap layer (r.blend + r.composite) and draped over DEM in wxNVIZ of GRASS GIS 7 (click to enlarge). The sunday walk path around Malga Brigolina is the blue/red vector line shown in the view center.

Situation at 4:00pm: we are close to sunset in Trentino… view towards the Rotaliana (Mezzocorona, S. Michele all’Adige), last sunlit summits also seen in photo:

rotaliana_direct_sunlight_15dec2013_16pm_foto(click to enlarge)

3. Outcome

The resulting sunlight/shadow map appear to match nicely realty. Perhaps r.sun.hourly should get an additional flag to generate the binary “sun-yes” – “sun-no” maps directly.

trento_direct_sunlight_15dec2013_15pm_3DDirect sunlight zones (yellow, 15 Dec 2013, 3pm): Trento with Monte Bondone, Paganella, Marzolan, Lago di Caldonazzo, Lago di Levico and surroundings (click to enlarge)

4. GRASS GIS usage note

The wxGUI settings were as simple as this (note the transparency values for the various layers):


Data sources:

The post Will the sun shine on us? appeared first on GFOSS Blog | GRASS GIS Courses.

Using the 25m EU-DEM for shading OpenStreetMap layers

Inspired by Václav Petráš posting about “Did you know that you can see streets of downtown Raleigh in elevation data from NC sample dataset?” I wanted to try the new GRASS GIS 7 Addon r.shaded.pca which creates shades from various directions and combines then into RGB composites just to see what happens when using the new EU-DEM at 25m.

To warm up, I registered the “normally” shaded DEM (previously generated with gdaldem) with r.external in a GRASS GIS 7 location (EPSG 3035, LAEA) and overlayed the OpenStreetMap layer using WMS with GRASS 7′s An easy task thanks to University of Heidelberg’s Indeed, they offer a similar shading via WMS, however, in the screenshot below you see the new EU data being used for controlling the light on our own:

OpenStreetMap shaded with EU DEM 25m

OpenStreetMap shaded with EU DEM 25m (click to enlarge)

Next item: trying r.shaded.pca… It supports multi-core calculation and the possibility to strengthen the effects through z-rescaling. In my example, I used:

r.shaded.pca input=eu_dem_25 output=eu_dem_25_shaded_pca nproc=3 zmult=50

The leads to a colorized hillshading map, again with the OSM data on top (50% transparency):


OpenStreetMap shaded with r.shaded.pca using EU DEM 25m (click to enlarge)

Yes, fun – I like it :-)

Data sources:

The post Using the 25m EU-DEM for shading OpenStreetMap layers appeared first on GFOSS Blog | GRASS GIS Courses.

EU-DEM: new Digital Surface Model at 25m


EU DEM 25m upper Garda Lake area with Riva del Garda and Arco (Italy). 3D view in wxNVIZ – GRASS GIS 7

The 25m European Digital Elevation Model (EU-DEM, Version 1) is a Digital Surface Model (DSM) representing the first surface as illuminated by the sensors:


EU DEM Rotaliana with Mezzocorona and S. Michele (Italy). Produced using Copernicus data and information funded by the European Union – EU-DEM layers.

Its elevations were captured at 1 arc second postings (2.78E-4 degrees). The tiles are provided at 25m resolution in EU-LAEA (EPSG. 3035) projection, temporal coverage: 2000, published in Oct 2013. It is a realisation of the Copernicus programme, managed by the European Commission, DG Enterprise and Industry. Metadata are provided here. According to their “Methodology” page it is a hybrid product based on SRTM and ASTER GDEM data fused by a weighted averaging approach and it has been generated as a contiguous dataset divided into 1 degree by 1 degree tiles, corresponding to the SRTM naming convention. In addition to the DEM data, a colour shaded relief image over Europe is provided.

From the metadata page: “The EU-DEM data are provided as is, i.e. without a formal validation yet. An independent statistical validation is scheduled as part of the GIO land monitoring service activities, and will be made available in the course of 2014.

1. Data download

Note that the GeoTIFF files are of major size, up to 5 GB:

Hint for GRASS GIS users: instead of importing the data, you can use r.external to simply register the GeoTIFF file within a EU LAEA projected location.


The post EU-DEM: new Digital Surface Model at 25m appeared first on GFOSS Blog | GRASS GIS Courses.

Generating contours using GDAL ( via shell or QGIS)

I tell you. It always amazes me how much cool stuff you can do with great open source GIS software these days. GDAL is one of those great open source projects that I have just found a great use for (apart from just opening every raster type under the sun in QGIS).

GDAL has the ability to generate contours from a DEM, something that I have always wanted to try for my town but have never been able due to lack of a good DEMs.

Recently we purchased a set of DEMs that cover a large area as part of a study. Each DEM uses a grid size of 1mx1m. Prefect for generating contours.

Using the GdalTools QGIS plugin.

First make sure that you have the latest version of the GdalTools plugin installed (GdalTools should be installed by default with QGIS. If it’s not, search “Gdal” in the plugin installer). Enable the plugin once it’s installed.

Load the DEM into QGIS using the Load Raster icon.

DEM loaded in QGIS

Now head up to the menu Raster->Extraction->Contour

Raster menu in QGIS

Select the settings that you need. For this DEM I am going to generate 250mm contours.

Contour dialog.

Take note of the text area at the bottom of the dialog as that is the exact command sent to GDAL in order to generate the contours. If you take a copy of that you can run it on the command line for batch processing later.

Hit ok.

250mm contours from the DEM

BAM! :)

Using the command line/shell.

Using QGIS for a one off DEM is fine and dandy but what if you have 3000 DEMs that you need to process. To hell with doing that by hand!

Remember the contour tool in QGIS told us the exact command line args to use, so creating a shell script for automating the process is pretty easy.

for f in *.asc
  echo "Processing $f"
 gdal_contour -a ELEV -i 0.25 $f $f-250mm.shp

The above code will loop though the current directory and process all the DEMs generating 250mm contours for each one. It saves each new contour file as {filename}-250mm.shp. You will need to change *.asc to whatever format your DEM is in.

Copy the above code into a file somewhere and call it This can then be called from the command line using


Running sh on Windows

If you’re a windows user you will need to run OsGeo4W Shell in order to use sh.

Loading OsGeo4w shell.

Once the shell is loaded you can just call:


Output from running

Final remarks

I ran the above sh script on a folder with about 2500 DEMs and it took around 4 hours to complete the whole folder. Of course performance will vary but it seems pretty quick considering.

Once again the possibly to use open source GIS tools to get my work done is amazing.  No expensive software here.

So far I’m not aware of any line smoothing/generalizing abilities using GDAL/OGR.  Although you can import the contours into GRASS and use that:

You can also generate the contours using GDAL via Python but that is a topic for another day.

Filed under: Open Source, qgis Tagged: contours, DEM, FOSSGIS, gdal, gis, grid, ogr, Open Source, osgeo, qgis, Quantum GIS, raster

  • Page 1 of 1 ( 5 posts )
  • dem

Back to Top