1 /**********************************************************
2 * Version $Id$
3 *********************************************************/
4
5 ///////////////////////////////////////////////////////////
6 // //
7 // SAGA //
8 // //
9 // System for Automated Geoscientific Analyses //
10 // //
11 // Tool Library //
12 // ta_lighting //
13 // //
14 //-------------------------------------------------------//
15 // //
16 // HillShade.cpp //
17 // //
18 // Copyright (C) 2003 by //
19 // Olaf Conrad //
20 // //
21 //-------------------------------------------------------//
22 // //
23 // This file is part of 'SAGA - System for Automated //
24 // Geoscientific Analyses'. SAGA is free software; you //
25 // can redistribute it and/or modify it under the terms //
26 // of the GNU General Public License as published by the //
27 // Free Software Foundation, either version 2 of the //
28 // License, or (at your option) any later version. //
29 // //
30 // SAGA is distributed in the hope that it will be //
31 // useful, but WITHOUT ANY WARRANTY; without even the //
32 // implied warranty of MERCHANTABILITY or FITNESS FOR A //
33 // PARTICULAR PURPOSE. See the GNU General Public //
34 // License for more details. //
35 // //
36 // You should have received a copy of the GNU General //
37 // Public License along with this program; if not, see //
38 // <http://www.gnu.org/licenses/>. //
39 // //
40 //-------------------------------------------------------//
41 // //
42 // e-mail: oconrad@saga-gis.org //
43 // //
44 // contact: Olaf Conrad //
45 // Institute of Geography //
46 // University of Goettingen //
47 // Goldschmidtstr. 5 //
48 // 37077 Goettingen //
49 // Germany //
50 // //
51 ///////////////////////////////////////////////////////////
52
53 //---------------------------------------------------------
54
55
56 ///////////////////////////////////////////////////////////
57 // //
58 // //
59 // //
60 ///////////////////////////////////////////////////////////
61
62 //---------------------------------------------------------
63 #include "HillShade.h"
64
65
66 ///////////////////////////////////////////////////////////
67 // //
68 // //
69 // //
70 ///////////////////////////////////////////////////////////
71
72 //---------------------------------------------------------
CHillShade(void)73 CHillShade::CHillShade(void)
74 {
75 //-----------------------------------------------------
76 Set_Name (_TL("Analytical Hillshading"));
77
78 Set_Author ("O.Conrad, V.Wichmann (c) 2003-2013");
79
80 Set_Description(_TW(
81 "This tool performs an analytical hillshade computation for an elevation grid. "
82 "The 'Standard' method simply calculates the angle at which light coming from the "
83 "position of the light source would hit the surface. This method can produce angles "
84 "greater than 90 degree. With the second method all values are kept within the "
85 "range of 0-90 degree. It can be enhanced with shadowing effects, where shadowed "
86 "cells will be marked with a value of exactly 90 degree. 'Shadows Only' creates "
87 "a mask for the shadowed areas and sets all other cells to no-data. 'Combined Shading' "
88 "takes the values of the standard method and multiplies these with the normalized slope. "
89 "'Ambient Occlusion' is based on the concepts of Tarini et al. (2006), but only "
90 "the northern half-space is considered here. "
91 ));
92
93 Add_Reference(
94 "Tarini, M. / Cignoni, P. / Montani, C.", "2006",
95 "Ambient Occlusion and Edge Cueing to Enhance Real Time Molecular Visualization",
96 "IEEE Transactions on Visualization and Computer Graphics, Vol. 12, No. 5, pp. 1237-1244."
97 );
98
99 //-----------------------------------------------------
100 Parameters.Add_Grid("",
101 "ELEVATION" , _TL("Elevation"),
102 _TL(""),
103 PARAMETER_INPUT
104 );
105
106 Parameters.Add_Grid("",
107 "SHADE" , _TL("Analytical Hillshading"),
108 _TL("The angle between the surface and the incoming light beams, measured in radians."),
109 PARAMETER_OUTPUT
110 );
111
112 //-----------------------------------------------------
113 Parameters.Add_Choice("",
114 "METHOD" , _TL("Shading Method"),
115 _TL(""),
116 CSG_String::Format("%s|%s|%s|%s|%s|%s|",
117 _TL("Standard"),
118 _TL("Limited Maximum"),
119 _TL("With Shadows"),
120 _TL("Shadows Only"),
121 _TL("Ambient Occlusion"),
122 _TL("Combined Shading")
123 ), 0
124 );
125
126 Parameters.Add_Choice("",
127 "POSITION" , _TL("Sun's Position"),
128 _TL(""),
129 CSG_String::Format("%s|%s|",
130 _TL("azimuth and height"),
131 _TL("date and time")
132 ), 0
133 );
134
135 Parameters.Add_Double("POSITION",
136 "AZIMUTH" , _TL("Azimuth"),
137 _TL("Direction of the light source, measured in degree clockwise from the North direction."),
138 315.0, 0.0, true, 360.0, true
139 );
140
141 Parameters.Add_Double("POSITION",
142 "DECLINATION" , _TL("Height"),
143 _TL("Height of the light source, measured in degree above the horizon."),
144 45.0, 0.0, true, 90.0, true
145 );
146
147 Parameters.Add_Date("POSITION",
148 "DATE" , _TL("Day"),
149 _TL(""),
150 CSG_DateTime::Now().Get_JDN()
151 );
152
153 Parameters.Add_Double("POSITION",
154 "TIME" , _TL("Hour"),
155 _TL(""),
156 12.0, 0.0, true, 24.0, true
157 );
158
159 Parameters.Add_Double("",
160 "EXAGGERATION" , _TL("Exaggeration"),
161 _TL("The terrain exaggeration factor allows one to increase the shading contrasts in flat areas."),
162 1.0
163 );
164
165 Parameters.Add_Choice("",
166 "UNIT" , _TL("Unit"),
167 _TL(""),
168 CSG_String::Format("%s|%s|",
169 _TL("radians"),
170 _TL("degree")
171 ), 0
172 );
173
174 Parameters.Add_Choice("",
175 "SHADOW" , _TL("Shadow"),
176 _TL("Choose 'slim' to trace grid node's shadow, 'fat' to trace the whole cell's shadow. The first is slightly faster but might show some artifacts."),
177 CSG_String::Format("%s|%s|",
178 _TL("slim"),
179 _TL("fat")
180 ), 0
181 );
182
183 Parameters.Add_Int("",
184 "NDIRS" , _TL("Number of Directions"),
185 _TW("Number of sample directions for ambient occlusion. Divides azimuth range (270 to 0 to 90 degree) into sectors. "
186 "Declination (0 to 90 degree) is divided in (Number of Directions / 4) sectors."),
187 8, 2, true
188 );
189
190 Parameters.Add_Double("",
191 "RADIUS" , _TL("Search Radius"),
192 _TL("Radius used to trace for shadows (ambient occlusion) [map units]."),
193 10.0, 0.001, true
194 );
195 }
196
197
198 ///////////////////////////////////////////////////////////
199 // //
200 ///////////////////////////////////////////////////////////
201
202 //---------------------------------------------------------
On_Parameters_Enable(CSG_Parameters * pParameters,CSG_Parameter * pParameter)203 int CHillShade::On_Parameters_Enable(CSG_Parameters *pParameters, CSG_Parameter *pParameter)
204 {
205 if( pParameter->Cmp_Identifier("METHOD") )
206 {
207 pParameters->Set_Enabled("POSITION" , pParameter->asInt() != 4);
208 pParameters->Set_Enabled("EXAGGERATION", pParameter->asInt() != 4 && pParameter->asInt() != 3);
209 pParameters->Set_Enabled("UNIT" , pParameter->asInt() < 3);
210 pParameters->Set_Enabled("SHADOW" , pParameter->asInt() == 2 || pParameter->asInt() == 3);
211 pParameters->Set_Enabled("NDIRS" , pParameter->asInt() == 4);
212 pParameters->Set_Enabled("RADIUS" , pParameter->asInt() == 4);
213 }
214
215 if( pParameter->Cmp_Identifier("POSITION") )
216 {
217 pParameters->Set_Enabled("AZIMUTH" , pParameter->asInt() == 0);
218 pParameters->Set_Enabled("DECLINATION" , pParameter->asInt() == 0);
219 pParameters->Set_Enabled("DATE" , pParameter->asInt() == 1);
220 pParameters->Set_Enabled("TIME" , pParameter->asInt() == 1);
221 }
222
223 return( CSG_Tool_Grid::On_Parameters_Enable(pParameters, pParameter) );
224 }
225
226
227 ///////////////////////////////////////////////////////////
228 // //
229 ///////////////////////////////////////////////////////////
230
231 //---------------------------------------------------------
On_Execute(void)232 bool CHillShade::On_Execute(void)
233 {
234 //-----------------------------------------------------
235 m_pDEM = Parameters("ELEVATION")->asGrid();
236 m_pShade = Parameters("SHADE" )->asGrid();
237
238 //-----------------------------------------------------
239 bool bResult;
240
241 switch( Parameters("METHOD")->asInt() )
242 {
243 default: bResult = Get_Shading(false, false); break; // Standard
244 case 1: bResult = Get_Shading( true, false); break; // Limited Maximum
245 case 5: bResult = Get_Shading(false, true); break; // Combined Shading
246 case 2: bResult = Get_Shadows(false ); break; // With Shadows
247 case 3: bResult = Get_Shadows( true ); break; // Shadows Only
248 case 4: bResult = AmbientOcclusion( ); break; // Ambient Occlusion
249 }
250
251 if( !bResult )
252 {
253 return( false );
254 }
255
256 //-----------------------------------------------------
257 if( Parameters("METHOD")->asInt() >= 3 )
258 {
259 m_pShade->Set_Unit(_TL(""));
260 }
261 else if( Parameters("UNIT")->asInt() == 0 )
262 {
263 m_pShade->Set_Unit(_TL("radians"));
264 }
265 else
266 {
267 m_pShade->Set_Unit(_TL("degree"));
268
269 m_pShade->Multiply(M_RAD_TO_DEG);
270 }
271
272 //-----------------------------------------------------
273 if( Parameters("METHOD")->asInt() == 3 ) // Shadows Only
274 {
275 DataObject_Set_Parameter(m_pShade, "UNISYMBOL_COLOR", (int)SG_COLOR_BLACK);
276 DataObject_Set_Parameter(m_pShade, "COLORS_TYPE", 0); // Single Symbol
277 }
278 else
279 {
280 DataObject_Set_Colors (m_pShade, 11, SG_COLORS_BLACK_WHITE, true);
281 DataObject_Set_Parameter(m_pShade, "COLORS_TYPE", 3); // Graduated Colors
282 }
283
284 //-----------------------------------------------------
285 return( bResult );
286 }
287
288
289 ///////////////////////////////////////////////////////////
290 // //
291 ///////////////////////////////////////////////////////////
292
293 //---------------------------------------------------------
Get_Position(double & Azimuth,double & Decline)294 bool CHillShade::Get_Position(double &Azimuth, double &Decline)
295 {
296 if( Parameters("POSITION")->asInt() == 0 )
297 {
298 Azimuth = Parameters("AZIMUTH" )->asDouble() * M_DEG_TO_RAD;
299 Decline = Parameters("DECLINATION")->asDouble() * M_DEG_TO_RAD;
300 }
301 else
302 {
303 CSG_Shapes Origin(SHAPE_TYPE_Point), GCS;
304
305 Origin.Get_Projection().Create(m_pDEM->Get_Projection());
306 Origin.Add_Shape()->Add_Point(m_pDEM->Get_System().Get_Extent().Get_Center());
307
308 SG_RUN_TOOL_ExitOnError("pj_proj4", 2, // Coordinate Transformation (Shapes)
309 SG_TOOL_PARAMETER_SET("SOURCE", &Origin)
310 && SG_TOOL_PARAMETER_SET("TARGET", &GCS)
311 )
312
313 TSG_Point Center = GCS.Get_Shape(0)->Get_Point(0);
314
315 CSG_DateTime Date(Parameters("DATE")->asDate()->Get_Date());
316
317 double Hour = Parameters("TIME")->asDouble();
318
319 Date.Set_Hour(Hour);
320 Date.Set(floor(Date.Get_JDN()) - 0.5 + Hour / 24.0); // relate to UTC, avoid problems with daylight saving time
321
322 SG_Get_Sun_Position(Date, 0.0, Center.y * M_DEG_TO_RAD, Decline, Azimuth);
323
324 Message_Fmt("\n%s: %f", _TL("Longitude"), Center.x);
325 Message_Fmt("\n%s: %f", _TL("Latitude" ), Center.y);
326 Message_Fmt("\n%s: %f", _TL("Azimuth" ), Azimuth * M_RAD_TO_DEG);
327 Message_Fmt("\n%s: %f", _TL("Height" ), Decline * M_RAD_TO_DEG);
328 }
329
330 return( Decline >= 0.0 );
331 }
332
333 //---------------------------------------------------------
Get_Shading(bool bDelimit,bool bCombine)334 bool CHillShade::Get_Shading(bool bDelimit, bool bCombine)
335 {
336 double Azimuth, Decline;
337
338 if( !Get_Position(Azimuth, Decline) )
339 {
340 return( false );
341 }
342
343 double sinDecline = sin(Decline);
344 double cosDecline = cos(Decline);
345
346 double Scale = Parameters("EXAGGERATION")->asDouble();
347
348 //-----------------------------------------------------
349 for(int y=0; y<Get_NY() && Set_Progress(y); y++)
350 {
351 #pragma omp parallel for
352 for(int x=0; x<Get_NX(); x++)
353 {
354 double Slope, Aspect;
355
356 if( !m_pDEM->Get_Gradient(x, y, Slope, Aspect) )
357 {
358 m_pShade->Set_NoData(x, y);
359 }
360 else
361 {
362 if( Scale != 1.0 )
363 {
364 Slope = atan(Scale * tan(Slope));
365 }
366
367 double d = M_PI_090 - Slope;
368
369 d = acos(sin(d) * sinDecline + cos(d) * cosDecline * cos(Aspect - Azimuth));
370
371 if( bDelimit && d > M_PI_090 )
372 {
373 d = M_PI_090;
374 }
375
376 if( bCombine )
377 {
378 d *= Slope / M_PI_090;
379 }
380
381 m_pShade->Set_Value(x, y, d);
382 }
383 }
384 }
385
386 return( true );
387 }
388
389
390 ///////////////////////////////////////////////////////////
391 // //
392 ///////////////////////////////////////////////////////////
393
394 //---------------------------------------------------------
395 #define SHADOW_SLIM 0x0
396 #define SHADOW_FAT_X 0x1
397 #define SHADOW_FAT_Y 0x2
398
399 #define EPSILON 0.0001
400
401 //---------------------------------------------------------
Get_Shadows(bool bMask)402 bool CHillShade::Get_Shadows(bool bMask)
403 {
404 double Azimuth, Decline;
405
406 if( !Get_Position(Azimuth, Decline) )
407 {
408 return( false );
409 }
410
411 //-----------------------------------------------------
412 int Shadow = SHADOW_SLIM;
413 double dx = sin(Azimuth + M_PI_180);
414 double dy = cos(Azimuth + M_PI_180);
415
416 if( fabs(dx) - fabs(dy) > EPSILON )
417 {
418 dy /= fabs(dx);
419 dx = dx < 0 ? -1 : 1;
420
421 if( Parameters("SHADOW")->asInt() && fabs(dy) > EPSILON )
422 Shadow = SHADOW_FAT_X;
423 }
424 else if( fabs(dy) - fabs(dx) > EPSILON )
425 {
426 dx /= fabs(dy);
427 dy = dy < 0 ? -1 : 1;
428
429 if( Parameters("SHADOW")->asInt() && fabs(dx) > EPSILON )
430 Shadow = SHADOW_FAT_Y;
431 }
432 else // diagonal
433 {
434 dx = dx < 0 ? -1 : 1;
435 dy = dy < 0 ? -1 : 1;
436
437 if( Parameters("SHADOW")->asInt() )
438 Shadow = SHADOW_FAT_X|SHADOW_FAT_Y;
439 }
440
441 double dz = tan(Decline) * sqrt(dx*dx + dy*dy) * Get_Cellsize();
442
443 //-----------------------------------------------------
444 if( bMask )
445 {
446 m_pShade->Assign_NoData();
447 }
448 else
449 {
450 Get_Shading(true, false);
451 }
452
453 //-----------------------------------------------------
454 for(int y=0; y<Get_NY() && Set_Progress(y); y++)
455 {
456 for(int x=0; x<Get_NX(); x++)
457 {
458 if( !m_pDEM->is_NoData(x, y) )
459 {
460 Set_Shadow_Trace(x, y, m_pDEM->asDouble(x, y), dx, dy, dz, Shadow);
461 }
462 }
463 }
464
465 return( true );
466 }
467
468 //---------------------------------------------------------
Set_Shadow_Trace(double x,double y,double z,double dx,double dy,double dz,int Shadow)469 bool CHillShade::Set_Shadow_Trace(double x, double y, double z, double dx, double dy, double dz, int Shadow)
470 {
471 for(x+=dx+0.5, y+=dy+0.5, z-=dz; ; x+=dx, y+=dy, z-=dz)
472 {
473 int ix = (int)x, iy = (int)y;
474
475 if( !is_InGrid(ix, iy) )
476 {
477 return( true );
478 }
479
480 if( !m_pDEM->is_NoData(ix, iy) )
481 {
482 if( z < m_pDEM->asDouble(ix, iy) )
483 {
484 return( true );
485 }
486
487 m_pShade->Set_Value(ix, iy, M_PI_090);
488
489 if( Shadow & SHADOW_FAT_X )
490 {
491 int xx = x - ix < 0.5 ? ix - 1 : ix + 1;
492
493 if( m_pDEM->is_InGrid(xx, iy) && z < m_pDEM->asDouble(xx, iy) )
494 {
495 m_pShade->Set_Value(xx, iy, M_PI_090);
496 }
497 }
498
499 if( Shadow & SHADOW_FAT_Y )
500 {
501 int yy = y - iy < 0.5 ? iy - 1 : iy + 1;
502
503 if( m_pDEM->is_InGrid(ix, yy) && z < m_pDEM->asDouble(ix, yy) )
504 {
505 m_pShade->Set_Value(ix, yy, M_PI_090);
506 }
507 }
508 }
509 }
510
511 return( false );
512 }
513
514
515 ///////////////////////////////////////////////////////////
516 // //
517 ///////////////////////////////////////////////////////////
518
519 //---------------------------------------------------------
AmbientOcclusion(void)520 bool CHillShade::AmbientOcclusion(void)
521 {
522 double dRadius = Parameters("RADIUS")->asDouble();
523 int iDirs = Parameters("NDIRS" )->asInt();
524
525 CSG_Points_Z Directions; Directions.Set_Count(iDirs);
526
527 for(int i=0; i<iDirs; i++)
528 {
529 Directions[i].z = (M_PI_180 * i) / iDirs; // only northern half-space, otherwise: (M_PI_360 * i) / iDirs;
530 Directions[i].x = sin(Directions[i].z - M_PI_090);
531 Directions[i].y = cos(Directions[i].z - M_PI_090);
532 }
533
534 m_pShade->Assign(0.0);
535
536 for(int y=0; y<Get_NY() && Set_Progress(y); y++)
537 {
538 #pragma omp parallel for
539 for(int x=0; x<Get_NX(); x++)
540 {
541 double s, a;
542
543 if( !m_pDEM->Get_Gradient(x, y, s, a) )
544 {
545 m_pShade->Set_NoData(x, y);
546 continue;
547 }
548
549 CSG_Point_Z Normal(sin(s) * sin(a), sin(s) * cos(a), cos(s));
550
551 //-----------------------------------------------------
552 for(int i=0; i<Directions.Get_Count(); i++)
553 {
554 for(int j=0; j<Directions.Get_Count(); j++)
555 {
556 Directions[i].z = sin((M_PI_090 * j) / (iDirs / 4.0));
557
558 double dAngle = (Normal.Get_X() * Directions[i].x + Normal.Get_Y() * Directions[i].y + Normal.Get_Z() * Directions[i].z);
559
560 if (dAngle <= 0.0)
561 continue;
562
563 if( AmbientOcclusion_Trace(x, y, Directions[i], dRadius) )
564 {
565 m_pShade->Add_Value(x, y, dAngle);
566 }
567 }
568 }
569
570 if( !m_pShade->is_NoData(x, y) )
571 m_pShade->Set_Value(x, y, M_PI - m_pShade->asDouble(x, y) / (iDirs * (iDirs / 4.0)));
572 }
573 }
574
575 return( true );
576 }
577
578 //---------------------------------------------------------
AmbientOcclusion_Trace(int x,int y,CSG_Point_Z Direction,double dRadius)579 bool CHillShade::AmbientOcclusion_Trace(int x, int y, CSG_Point_Z Direction, double dRadius)
580 {
581 double dDistance, iDistance, dx, dy, dz, ix, iy, iz, z;
582
583 z = m_pDEM->asDouble(x, y);
584 dx = Direction.Get_X();
585 dy = Direction.Get_Y();
586 dz = tan( asin(Direction.Get_Z()) ) * sqrt(dx*dx + dy*dy) * Get_Cellsize();
587 ix = x;
588 iy = y;
589 iz = m_pDEM->asDouble(x, y);
590 dDistance = 0.0;
591 iDistance = Get_Cellsize() * M_GET_LENGTH(dx, dy);
592
593 while( is_InGrid(x, y) && dDistance <= dRadius )
594 {
595 ix += dx; x = (int)(0.5 + ix);
596 iy += dy; y = (int)(0.5 + iy);
597 iz += dz;
598
599 if( m_pDEM->is_InGrid(x, y) && m_pDEM->asDouble(x, y) > iz )
600 return( false );
601
602 dDistance += iDistance;
603 }
604
605 return( true );
606 }
607
608
609 ///////////////////////////////////////////////////////////
610 // //
611 // //
612 // //
613 ///////////////////////////////////////////////////////////
614
615 //---------------------------------------------------------
616