[Gdal-dev] RE: shaded relief utility (DIY)

[Charles Wivell]

Oops;

You can't use the the GCA, in this case, that was for a point on the surface of the Earth... You'll have to use the dot product between the normal and a unit vector from the point to the sun's location.

cos(Z) = N dot Vs / |N|

Sorry, its early in the morning here...

Chuck



>>> "Charles Wivell" <charles.wivell at dnr.state.mn.us> 12/13/05 8:54 AM >>>

Hi all;

If you want to have some fun coding. Creating a simple shaded relief image from a DEM is pretty simple to do. I'll outline it below:

Once you have the DEM read in... Create an output space buffer the same size (line,sample, bytes). Step through the postings starting at posting p(0,0). p(line,sample)

Create a 3D vector from posting p(0,0) to p(1,0) = (vecA) and another one from p(0,0) to p(0,1) = (vecB), Taking the cross product of these two vectors vecA x vecB = normal to the plane created by these three points. Now place the sun where ever you want (lat,lon). The brightness of that plane is proportional to the cosine of the zenith angle (cos(Z)), so a quick way to find cos(Z) is to use a form of the law of cosines in spherical trig (which is really just the Great Circle Angle between the sun's nadir and the point p(0,0)). or:

cos(Z) = sin(Latp) * sin(Lats) + cos(Latp) * cos(Lats) * cos(Lonp - Lons);

where Latp = latitude of the p(0,0) and Lonp = longitude of p(0,0) and
Lats = latitude of the sun's nadir and Lons = longitude of the sun's Nadir;

>From there you just have to scale cos(Z) from 0 to 255 or B = 255 * cos(Z)... Or just play with the scaling to get what you like. Place the brightness value in the output space buffer at (0,0). And continue to step through the postings.

If you want to include shadow, it becomes a little more complex.

If you want to get more complex, you can use the points around the point in question to fit a surface and then find the normal of that surface at the point in question... Anyway, I've found that the difference between the fitted surface and using the cross product is small (in most cases) and the cross product gives good qualitative results.


Chuck Wivell
MN DNR Grand Rapids, MN



>>> "Ethan Alpert" <ealpert at digitalglobe.com> 12/12/05 5:04 PM >>>


Thanks a lot! Works great!

-----Original Message-----
From: gdal-dev-bounces at lists.maptools.org 
[mailto:gdal-dev-bounces at lists.maptools.org] On Behalf Of Matthew Perry
Sent: Monday, December 12, 2005 2:16 AM
To: Gdal-Dev
Subject: [Gdal-dev] shaded relief utility


Hey folks,

  If anyone's interested, I created a GDAL-based utility for creating
shaded relief maps from a DEM. This was for a client who was kind enough
to let me keep the rights to all of my code ;-) So I'll release it to
the public domain in case anyone can find a use for it.

  Basically it's a simple GDAL/C++  port of the GRASS algorithm from
r.shaded.relief. Of course this has the advantage of being able to use
any GDAL DEM not just a GRASS map.  I've compiled it against gdal 1.3.1
on linux and I'd like to see if anyone else has any luck on different
platforms. Of course suggestions and changes are very appreciated as my
C++ experience is a bit weak.

--
Matt Perry
perrygeo at gmail.com 
http://www.perrygeo.net 

_______________________________________________
Gdal-dev mailing list
Gdal-dev at lists.maptools.org 
http://lists.maptools.org/mailman/listinfo/gdal-dev 


_______________________________________________
Gdal-dev mailing list
Gdal-dev at lists.maptools.org 
http://lists.maptools.org/mailman/listinfo/gdal-dev