Geographic Information Systems Stack Exchange is a question and answer site for cartographers, geographers and GIS professionals. Join them; it only takes a minute:

Sign up
Here's how it works:
  1. Anybody can ask a question
  2. Anybody can answer
  3. The best answers are voted up and rise to the top

I'm looking to plot day/night on a Google map, for an arbitrary point in time. I'm familiar with generating map tiles; I'm just looking for an algorithm to tell me whether a particular point on the globe is currently in daylight or darkness, or to otherwise plot the curve of the day/night interface onto the map.

I've done some searching, but it's possible I don't even know enough about the problem domain here to know what terms to search for!

Any ideas? Doesn't have to be perfect -- basically, I'm comparing Flickr geolocation data of sunrise and sunset photos (and their "date taken" timestamps) with reality, and this is to help me visualise it.

share|improve this question
Solutions also appear in a very closely related question at…. – whuber Oct 29 '14 at 23:44
up vote 5 down vote accepted

This page gives equations good to 1 degree. It looks like this code calculates it too but I didn't actually check.

share|improve this answer
Fantastic, just what I was looking for. And yes, looks like the projillum() function in vplanet.c in that code matches up pretty well with that algorithm, so that should definitely get me heading down the right track, thanks. – Matt Gibson Mar 21 '11 at 23:54
That link is broken. Please include relevant details instead of just a link. – captncraig Oct 26 '12 at 3:33
second link works though – iant Oct 26 '12 at 7:14

Also an example

    Sunrise/Sunset Algorithm

    Almanac for Computers, 1990
    published by Nautical Almanac Office
    United States Naval Observatory
    Washington, DC 20392

    day, month, year:      date of sunrise/sunset
    latitude, longitude:   location for sunrise/sunset
    zenith:                Sun's zenith for sunrise/sunset
      offical      = 90 degrees 50'
      civil        = 96 degrees
      nautical     = 102 degrees
      astronomical = 108 degrees

    NOTE: longitude is positive for East and negative for West
        NOTE: the algorithm assumes the use of a calculator with the
        trig functions in "degree" (rather than "radian") mode. Most
        programming languages assume radian arguments, requiring back
        and forth convertions. The factor is 180/pi. So, for instance,
        the equation RA = atan(0.91764 * tan(L)) would be coded as RA
        = (180/pi)*atan(0.91764 * tan((pi/180)*L)) to give a degree
        answer with a degree input for L.

1. first calculate the day of the year

    N1 = floor(275 * month / 9)
    N2 = floor((month + 9) / 12)
    N3 = (1 + floor((year - 4 * floor(year / 4) + 2) / 3))
    N = N1 - (N2 * N3) + day - 30

2. convert the longitude to hour value and calculate an approximate time

    lngHour = longitude / 15

    if rising time is desired:
      t = N + ((6 - lngHour) / 24)
    if setting time is desired:
      t = N + ((18 - lngHour) / 24)

3. calculate the Sun's mean anomaly

    M = (0.9856 * t) - 3.289

4. calculate the Sun's true longitude

    L = M + (1.916 * sin(M)) + (0.020 * sin(2 * M)) + 282.634
    NOTE: L potentially needs to be adjusted into the range [0,360) by adding/subtracting 360

5a. calculate the Sun's right ascension

    RA = atan(0.91764 * tan(L))
    NOTE: RA potentially needs to be adjusted into the range [0,360) by adding/subtracting 360

5b. right ascension value needs to be in the same quadrant as L

    Lquadrant  = (floor( L/90)) * 90
    RAquadrant = (floor(RA/90)) * 90
    RA = RA + (Lquadrant - RAquadrant)

5c. right ascension value needs to be converted into hours

    RA = RA / 15

6. calculate the Sun's declination

    sinDec = 0.39782 * sin(L)
    cosDec = cos(asin(sinDec))

7a. calculate the Sun's local hour angle

    cosH = (cos(zenith) - (sinDec * sin(latitude))) / (cosDec * cos(latitude))

    if (cosH >  1) 
      the sun never rises on this location (on the specified date)
    if (cosH < -1)
      the sun never sets on this location (on the specified date)

7b. finish calculating H and convert into hours

    if if rising time is desired:
      H = 360 - acos(cosH)
    if setting time is desired:
      H = acos(cosH)

    H = H / 15

8. calculate local mean time of rising/setting

    T = H + RA - (0.06571 * t) - 6.622

9. adjust back to UTC

    UT = T - lngHour
    NOTE: UT potentially needs to be adjusted into the range [0,24) by adding/subtracting 24

10. convert UT value to local time zone of latitude/longitude

    localT = UT + localOffset

using the algorithm

share|improve this answer
While this link may answer the question, it is better to include the essential parts of the answer here and provide the link for reference. Link-only answers can become invalid if the linked page changes. – iant Aug 14 '12 at 13:22

Your Answer


By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.