5

I'm trying to use GrassGIS to extract the shadow from a DEM (elevation data).

The DEM I'm using shows Geneva (plus a part of the Alps mountains), and was captured by the SRTM. The elevation data was downloaded from the USGS website, and was already projected to the WGS84 coordinate system:

srtm

The command I've used in GrassGIS to calculate the shadow is r.sunmask, for which the date, time, and timezone has to be provided:

r.sunmask elevation=n46_e006_1arc_v3@PERMANENT output=alps-sunmask.tif year=2014 month=9 day=23 hour=11 minute=0 timezone=1

The problem is that the image produced with r.sunmask is completely empty. Is something wrong with the command above?

P.S: I've checked the answer suggested here, but I'm unable to make r.sun work since I can't find which parameter correspond to the output.

Community
  • 1
Hakim
  • 956
  • 12
  • 33
  • have you tried reprojecting to a local projection in meters (e.g. epsg:21781) first so that the x/y/elevation are all in meters? otherwise the vertical exaggeration will be way out. – Steven Kay Aug 21 '17 at 16:40
  • @StevenKay The elevation seems to be already in meter. – Hakim Aug 21 '17 at 17:27

3 Answers3

7

You need x, y and z units in meters. Or at least, consistent units. If your raster is in WGS84 and elevations are in meters, they're not consistent. This can give weird results for a lot of algorithms.

If the units aren't consistent, the mountains are thousands of times steeper than you think they are - nearly vertical - and almost everything is in shadow... that would explain things if your raster is mostly or all values of 1

I picked a small area in the Swiss Alps from SRTM and used warp to convert it to epsg:21781 (using Cubic interpolation to avoid 'artifacts').

Then used r.sunmask.datetime . This example was 23rd September but I chose 9am to give more obvious shadows.

In my case I used it via processing/QGIS, which insists on easting/northings. Those two parameters should be optional, so I manually chose a point near the centre, using the coordinate capture plugin. If you're using this from the GRASS command line you don't need to do that.

EDIT note that if doing this outside of qgis/processing, you'll need to remember to set the region using g.region

Here's the result... enter image description here

And with the same settings, but done as WGS84 - almost everything in shadow

enter image description here

Steven Kay
  • 20,405
  • 5
  • 31
  • 82
  • You're right the elevation data I have is in the WGS84 projection. How long did take for you to extract the shadow (it's taking ages for me)? – Hakim Aug 28 '17 at 14:50
  • 1
    can take quite a long time, I tried this on an old laptop! It's why I only ran it on a small part of the SRTM raster ;-) I reckon it would have taken an hour or more, and the raster's only 6000x4000 or so. That small subsection took about 10 minutes I think. – Steven Kay Aug 28 '17 at 15:37
  • 1
    Could you just add in your response that the the region had to be set to the DEM's region using g.region (because without it the resulting shadow mask would be empty, i.e. No Data)? – Hakim Sep 04 '17 at 16:21
3

Please use r.sun and not r.sunmask for this. The underlying algorithm is even more precise.

The manual page of r.sun offers a dedicated section on "Extraction of shadow maps", explaining that a map of shadows can be extracted from the solar incidence angle map ("incidout" output parameter) in which areas with zero values are shadowed.

markusN
  • 12,976
  • 31
  • 48
  • Can I get the solar incidence angle map from the date, time of the day, and the timezone? I'm confused about what should be provided in that field (apparently it's required to make the command r.sun work). – Hakim Aug 21 '17 at 18:00
2

Maybe your method is correct but there is no shadow...

You are close to midday, so your sun elevation will be close to 40° (check with the "-g" option when you run r.sunmask ). This means that you need a slope of more than 40° in the direction of the sun in order to have a shadow. This is possible in mountain like the alps, but unlikely (especially with a 90m DEM, as it would mean about 75 m difference between two adjacent pixels).

Please check close to sunrise/sunset to verify that the application is working properly, because it seems to be OK at a first glance as long as your XYZ units are the same.

WARNING : for correct shadow location on a study area, you need to take into account the mountains that are outside of your study areaon the side of the sun.

radouxju
  • 49,636
  • 2
  • 71
  • 144