BLAINE BELLMANN'S

MOUNTAINOUS TERRAIN INSOLATION SIMULATION

A GEOG 472 PROJECT

Note: Some pics were lost when my UBC geography website was removed to make room for websites for other geography students.

TABLE OF CONTENTS

• INTRODUCTION
• DATA, EQUATIONS AND PARAMETERS USED
• METHODOLOGY
• SETTING UP THE DEM IN IDRISI
• CREATION OF LATITUDE FILE
• CREATION OF LONGITUDE FILE
• CREATION OF SLOPE AND ASPECT FILES
• CALCULATION OF SOLAR DECLINATION ANGLE
• CREATION OF SOLAR ELEVATION FILE
• CREATION OF COSINE OF SOLAR AZIMUTH FILE
• CREATION OF SOLAR AZIMUTH FILE
• CREATION OF INSOLATION INTENSITY FACTOR FILE
• CALCULATION OF SOLAR CONSTANT
• CREATION OF INSOLATION FILE
• THE FINAL PRODUCT!
• PROJECT DISCUSSION
• CONCLUSIONS
• BIBLIOGRAPHY

INTRODUCTION

In today's digital world of computers and Internet links, we are able to accomplish feats that were either impossible or too time consuming to be feasible not too long ago. There is an ever-increasing amount of information being stored in digital form, ready to be processed in a multitude of ways by an ever-increasing number of software products that are getting more powerful every day.

One type of useful data readily available online free of charge is Digital elevation models (DEM's) covering large areas of the earth (30 arc-second resolution), and the United States (3 arc-second resolution). In this project, I examine one use of the 30 arc-second resolution Digital Elevation Models in the field of atmospheric science by making use of a Geographic Information System software package called Idrisi for Windows and a number of equations to create a realistic simulation of the solar illumination of mountainous terrain. All of the relevant equations and Idrisi commands used are described in this report, followed by the resulting image from several perspectives, some sample insolation measurements and a discussion of the results, including how the resulting insolation simulation could be useful for studies in other fields than meteorology.

Since my home is British Columbia, Canada, I decided to work with 30 arc-second digital elevation data covering all of British Columbia and some surrounding areas. The data set I used came from the U.S. Geological Survey's FTP site, which is now edcwww.cr.usgs.gov/pub/data/gtopo30/global/ . The data used had the following attributes:

• Horizontal Resolution: 30 Arc-Seconds Or Roughly 1 km
• Horizontal Datum: WGS84
• Elevation Units: Meters Above Mean Sea Level

Since I was interested in the solar illumination of the topography, I had to select a time and date for the simulation. My choice was 1:00 PM Pacific Daylight Savings Time on June 22'nd, which was the middle of the summer solstice. What follows is a description of the equations I made use of for my mountainous terrain insolation simulation.

These equations come from the text "Meteorology Today For Scientists and Engineers" by Dr. Roland B. Stull.

This equation calculates the solar declination angle, which is the latitude where the sun is directly overhead. It assumes a circular orbit and is therefore only approximate:

Where:

• is solar declination angle in radians.
• is the angle between the plane of the earth's orbit (ecliptic) and the plane of the earth's equator. It is also the maximum latitude that the sun reaches overhead (23.45or 0.409 radians).
• C is the circle circumference (2 in radians).
• d is the number of days into the year (Julian Day).
• is the Julian day of the summer solstice (173 for non-leap years).
• is the number of days in a year (365 for non-leap year, 366 for leap year).

The next equation calculates the solar constant for any day of the year in Watts per square meter ():



Where:

• S is the calculated solar constant for that day.
• is the average solar constant (1368 ).
• is the average earth-sun distance, which is where is measured at (149.6 Gm or 149,600,000 km).
• R is the distance between the sun and the earth on a given day, calculated by the following equation:
  
  
  
  Where:
  • a is half the distance of the longest axis (semi-major axis) in the earth's elliptical orbit (149.457 Gm).
  • e is the eccentricity of the earth's orbit, which is a measure of how the earth's orbit deviates from a perfect circle. It is given by the equation:
    e=c/a=0.0167
    Where c is half of the distance between the two foci of the earth's elliptical orbit (2.5 Gm).
  • v is the true anomaly, which is the angle at the sun between the perihelion and the location of the earth-moon system's center of gravity (barycenter) and is given by this equation:
    
    v=M+0.0333988sin(M)+0.0003486sin(2M)+0.000005sin(3M)
    
    Where M is the mean anomaly, which is given by yet another equation:
    
    
    
    Where:
    • t is the Julian day.
    • is the perihelion day, which is the day the earth is closest to the sun (Jan 3).
    • P is the period of the earth's orbit (365.25463 days).

The next equation calculates the angle at which the sun is above the horizon (solar elevation):

Where:

• is the solar elevation in radians.
• is the latitude.
• is the longitude.
• is the Universal time (20 for 1:00 pm PDT).
• is the number of hours in the day (24).

The next equation calculates the cosine of the solar azimuth, where the azimuth is the direction in degrees (or radians) where:

• 0 (0 radians) is north.
• 90 (/2 radians) is east.
• 180 ( radians) is south.
• 270 (3/2 radians) is west:
  

Where is the solar azimuth angle in radians.

Where:

• d is the difference between the solar azimuth and the surface aspect in radians
• is the solar zenith angle in radians.
• is the surface slope in radians.
• is the insolation zenith angle in radians.

However, this simple model neglected a number of factors, including the effects of the atmosphere, the shading of other flatter surfaces by steeper topography and the terrestrial reflection on the incoming solar radiation. My choice of the summer solstice allowed me to avoid problems with the shading due to the high elevation of the sun.

METHODOLOGY

Here are the steps I took to create the insolation simulation model for BC using the equations defined above:

SETTING UP THE DEM IN IDRISI

The Idrisi format digital elevation model data was provided for me by my Geography 472 instructor, Brian Klinkenberg. He accessed the digital elevation model data from the U.S. Geological Survey's FTP site , which is now edcwww.cr.usgs.gov/pub/data/gtopo30/global/. He used ARC/INFO, a much more sophisticated GIS software package to import the large raw data file, windowing out the portion containing BC, and exporting it to Idrisi format, which was much more convenient since he had downloaded the data into a UNIX machine which could handle extremely LARGE data files. Memory limitations on the PC's would have created problems for processing the data using Idrisi's BILIDRIS import module and WINDOW sub-image extraction module.

The water was reassigned from -9999 to 0 meters above sea level:

Here are a couple of versions of the highly revered Digital Elevation Model:
Plain Greyscale DEM	Color Enhanced DEM
Click here for a medium size greyscale DEM!	Click here for a medium size color-enhanced DEM!
Click here for the FULL-SIZE (1046 K) greyscale DEM!	Click here for the FULL-SIZE (529 K) color-enhanced DEM!

CALCULATION OF SOLAR DECLINATION ANGLE

Making use of the previously defined equation , and plugging in the following variables for June 22'nd on a non-leap year:

• = 0.409 radians
• C = 2 radians
• = 173 days
• = 365 days
• d = 173 days (June 22'nd)

The solar elevation was calculated using the previously defined equation:



Where:

• C=2 radians
• =20 hours (for 1:00 Pacific Daylight Savings Time)
• =24 hours
• =5.235987757
• =0.40927971 radians
• =0.3979486315
• =0.9174077001

FUNCTION	INPUT FILE	OPERATION	OUTPUT FILE	SYMBOLIZED BY:
TRANSFOR	BCLA	SINE(BCLA)	SINLAT
TRANSFOR	BCLO	COSINE(BCLO)	COSLAT
SCALAR	BCLO	BCLO-5.235987757 Value units=RADIANS	BCLODIFF	-()

The sign was irrelevant due to cosine's symmetry about 0.

FUNCTION	INPUT FILE	OPERATION	OUTPUT FILE	SYMBOLIZED BY:
TRANSFOR	BCLODIFF	COSINE(BCLODIFF)	BCOSDIFF
SCALAR	BCOSDIFF	0.9174077001*BCOSDIFF	BCCOS2
SCALAR	SINLAT	SINLAT*0.3979486315	BCSIN1
OVERLAY	COSLAT BCCOS2	COSLAT*BCCOS2	BCCOS1
OVERLAY	BCSIN1 BCCOS1	BCSIN1-BCCOS1	SINELEV
TRANSFOR	SINELEV	ARCSINE(SINELEV) VALUE UNITS=RADIANS	SOLELEV

Here is what the solar elevation file looked like:

• =0.40927971 radians
• =0.3979486315
• C=2 radians
• C/4=1.570796327 radians

FUNCTION	INPUT FILE	OPERATION	OUTPUT FILE	SYMBOLIZED BY:
SCALAR	SOLELEV	SOLELEV-1.570796327	TEMPZEN
SCALAR	TEMPZEN	TEMPZEN*-1 Value Units=RADIANS	SOLZEN
TRANSFOR	SOLZEN	SINE(SOLZEN)	SINZEN
TRANSFOR	SOLZEN	COSINE(SOLZEN)	COSZEN
OVERLAY	SINZEN COSLAT	SINZEN*COSLAT	TMPDENOM
OVERLAY	SINLAT COSZEN	SINLAT*COSZEN	TMPNUM1
SCALAR	TMPNUM1	TMPNUM1-0.3979486315	TMPNUM2
SCALAR	TMPNUM2	TMPNUM2*-1	TMPNUM
OVERLAY	TMPNUM TMPDENOM	TMPNUM/TMPDENOM	COSAZIM

CREATION OF SOLAR AZIMUTH FILE

Before calculating the solar azimuth, a correction file had to be made to prevent the azimuth from doubling back for regions east of 120 west, the latitude of the sun's highest point at 1:00 PM PDT. The correction file shown below, had a value of 2 for regions to the east of 120 west, and 0 everywhere else. The solar azimuth file was then subtracted from the correction file to set things straight:

FUNCTION	INPUT FILE	OPERATION	OUTPUT FILE
RECLASS	BCLONG	0 to all values from 120 to just less than 150 2 to all values from 110 to just less than 120	CORRECTN
SCALAR	CORRECTN	CORRECTN*3.141592654	CORRECT2

Here is what the correction file looked like:

FUNCTION	INPUT FILE	OPERATION	OUTPUT FILE
TRANSFOR	COSAZIM	ARCCOSINE(COSAZIM)	TEMPAZIM
OVERLAY	CORRECT2 TEMPAZIM	CORRECT2-TEMPAZIM	TEMPAZM2
TRANSFOR	TEMPAZM2	ABS(TEMPAZM2) Value units: RADIANS	SOLAZIM

Here is what the solar azimuth file looked like:

---

FUNCTION	INPUT FILE	OPERATION	OUTPUT FILE	SYMBOLIZED BY
TRANSFOR	BCSLOPRD	COSINE(BCSLOPRD)	COSSLOPE
OVERLAY	COSZEN COSSLOPE	COSZEN*COSSLOPE	COSLOPZN
DOCUMENT	BCASPECT	Value units: DEGREES	BCASPECT
TRANSFOR	BCASPECT	RADIANS(BCASPECT)	BCASPCTR
OVERLAY	SOLAZIM BCASPCTR	SOLAZIM-BCASPCTR Value units: RADIANS	DIFFAZIM
TRANSFOR	DIFFAZIM	COSINE(DIFFAZIM)	COSDFAZM
TRANSFOR	BCSLOPRD	SINE(BCSLOPRD)	SINSLOPE
OVERLAY	COSDFAZM SINSLOPE	COSDFAZM*SINSLOPE	SNSLCSDA
OVERLAY	SNSLCSDA SINZEN	SNSLCSDA*SINZEN	TEMP
OVERLAY	COSLOPZN TEMP	COSLOPZN+TEMP	COSOLANG

The calculation of the solar constant for June 22'nd involved several equations. The first equation to be evaluated was the mean anomaly equation:

Where:

• t=173 days (June 22'nd)
• =3 (Jan 3)
• P=365.25463 days
• C=2 radians

M=2.924374983

M was then plugged into the equation for the true anomaly:

v=M+0.0333988sin(M)+0.0003486sin(2M)+0.000005sin(3M)=2.931429183

Which, in turn, was inserted into the equation for the sun to earth-moon barycenter distance for June 22'nd:



Where:

• a=149.457 Gm
• e=c/a=2.5 Gm/149.457 Gm=0.016727219
  

• =1368
• =149.6 Gm

THE FINAL PRODUCT!

Here it is! The product of all of those calculations. This is a full map of British Columbia showing the simulated solar illumination. The gradations seen on the ocean are due to the varying angle of the sun over different latitudes and longitudes.

Click here for the FULL-SIZE (1170 K) image!

Here is a closer view of southwestern BC cut out using Idrisi's sub-image extraction module WINDOW:

Sorry, but I was unable to recover the imagery below after setting up this page here after it was removed from the UBC server.

Here is a side view of southern BC (same as at the top of this document) using Idrisi's ORTHO module on a WINDOW'ed out portion of the insolation image:

The Idrisi format of the insolation image consisted of the actual values of insolation as calculated by the above procedure. Using Idrisi's "Cursor Inquiry Mode", the values of the direct insolation (assuming no atmosphere) could be read off any point on the map in . Here is an example of a few points in a zoomed section of the Sumas region in the lower Fraser valley:

Insolation at point A on north side of Sumas Mountain: 951.9475
Insolation at point B on south side of Sumas Mountain: 1293.42
Insolation at point C on valley bottom: 1195.418

PROJECT DISCUSSION

As previously mentioned, this model neglected a number of factors, such as atmospheric absorption and scattering, and terrestrial reflection of the sunlight as well as the potential shading of terrain whose slope and azimuth would have otherwise allowed for exposure to the sun. However, the combination of my choice of "high noon" on the summer solstice and the coarseness of the data (~1 km) "smoothing" out the steepness of the slopes resulted in no total shading. If it did occur, then the feature blocking out the sun would have had a negative insolation intensity factor on the shaded side, making the reclassification of negative values to 0 mentioned earlier necessary. If I had more time, I would have created a simulated insolation map for the winter solstice, December 22'nd. However, the shading problem would have been quite pronounced due to the sun's low elevation. Idrisi's SURFACE module has an "analytical hillshading" feature which has an effect similar to my results, except that it doesn't account for variations of the solar elevation and azimuth with latitude and longitude. Therefore, it is not an accurate representation of solar illumination of the region, limiting its uses to just illustration of topography.

The above procedure was extremely cumbersome, with each operation taking roughly 10 minutes to process 16 megabyte files. I never had time to investigate the use of Idrisi's macros, which would have "packaged" groups of operations quite nicely. However, the execution time would not have been shortened. Other GIS packages have the feature that allows for flexible user-specified equations to operate on several rasters to create a new raster. ArcGIS's ArcView's Map Calculator, and public domain GRASS's r.mapcalc program are two examples. This feature has since then been incorporated into all the newer versions of Idrisi as the Image Calculator, which has increased its raster image processing flexibility greatly. This feature would have significantly reduced the amount of processing time and work for the above equations. My work was done using Idrisi because its low-cost, easy-to-use features made it the software of choice for learning the concepts of GIS at UBC.

Creating realistic shading and being able to read the insolation off the map is only the tip of the iceberg of possible uses. The results could be the first step in the process of modelling the radiation budget in complex terrain, incorporating many other factors such as atmospheric and terrestrial absorption, emission, reflection and scattering, taking into account their variation with elevation. The 30 arc-second (~1 km) resolution of the DEM used is compatible with the resolution of the Advanced Very High Resolution Radiometer (AVHRR), allowing for the comparison between calculated insolation and thermal satellite imagery to obtain information about the varying albedo (reflectivity) of the surface as well as the atmospheric and terrestrial effects not accounted for.

The use of this information can extend well beyond the field of meteorology. An accurate map of varying solar intensity on the ground could be used to cancel out the shading on satellite imagery, resulting in a more uniformly illuminated surface, eliminating one problem encountered in image classification. Relating the spatial insolation intensity variation with other factors such as elevation, latitude, proximity to climate-moderating bodies of water, general climate type and soil type could result in maps showing the most suitable regions for growing certain crops.

The inspiration for this project came from the United Nations Institute for Training and Research (UNITAR) publication "Explorations in Geographic Information Systems Technology" Workbook Series, Volume 5: "GIS and Mountain Environments" Exercise 4: "Modeling the Impact of Mountains on Regional Climates". The exercise was part of a project in the Puna Atacama province in Chile, whose aim was to examine availability of water and its sources. It made use of digital elevation models of the area's mountainous terrain to model the development of heat islands responsible for heavy summer convective precipitation in the mountains to the northeast. This information was then used, along with drainage patterns, to help learn more about potential water sources.

CONCLUSIONS

There are a lot of potential uses for insolation models in a GIS. This project looked at how a simple insolation model could be created using GIS software, a DEM and some equations. Some useful applications were also examined. Idrisi is a relatively small GIS software package, compared to others such as ARC/INFO, and therefore required more work than that using the more sophisticated GIS software packages. With the continuing increase in the amount of available data and the ongoing development of more powerful computers and software packages, models that were once too intimidating due to complexity are presently not only within reach, but are now being implemented.

Stull, Roland B. 1995. "Meteorology Today For Scientists and Engineers", West Publishing Company, Minneapolis/St. Paul, New York, Los Angeles, San Francisco. pp. 42 - 45, ISBN 0-314-06471-0

Clark Labs 1995. "Explorations in Geographic Information Systems Technology" Workbook Series Volume 5: "GIS and Mountain Environments", Exercise 4: "Modeling the Impact of Mountains on Regional Climates". Kristin Schneider and Paul Robbins, editors. Published by the United Nations Institute for Training and Research (UNITAR). pp. 4-1  -  4-15.

Head back to "Blaine's World" Home Page!