Calculating an overlay of day/night for Google Maps

Question

I am trying to find a way to create an overlay for Google Maps API V3 that shows the sunlit areas of the world. This is the basic result I am looking for:

http://www.daylightmap.com/index.php

But want more control over the appearance (ideally just a 10% black overlay with no city lights). I can draw the shape in a canvas element but can not figure out how to calculate the shape based on earth's tilt and rotation etc.

Any help would be appreciated.

EDIT: Javascript

I still don't know where to implement the y-offset variable below. I also need to figure out how to adjust/stretch the y-offset from this (equal distant latitudinal lines) to mercator (closer at poles).

// Get the canvas element
var ctx = document.getElementById('canvas').getContext('2d');
ctx.clearRect( 0, 0, 800, 620 );

// Current time
var map_width = $("#canvas").width();
var map_height = $("#canvas").height();
var now = new Date();
var cur_hour = now.getHours();
var cur_min = now.getMinutes();
var cur_sec = now.getSeconds();
var cur_jul = now.julianDate() - 1;
var equinox_jul = new Date(now.getFullYear(),2,20,24,-now.getTimezoneOffset(),0,0).julianDate() - 1;

var offset_x = Math.round(((cur_hour*3600 + cur_min*60 + cur_sec)/86400) * map_width); // Resulting offset X
var offset_sin = ((365.25 - equinox_jul + cur_jul)%365.25)/365.25; // Day offset, mapped on the equinox offset
var offset_sin_factor = Math.sin(offset_sin * 2 * Math.PI); // Sine wave offset
var offset_y = offset_sin_factor * 23.44; // Map onto angle. Maximum angle is 23.44° in both directions

var degrees_per_radian = 180.0 / Math.PI;
var offset_y_mercator = Math.atan( offset_y.sinh() ) * degrees_per_radian;


// Global wave variables
var period = 1 / 6.28291;   // Original value 2Pi: 6.28291
var amplitude = (map_height/2);

// Draw vertical lines: One for each horizontal pixel on the map
for( var x = 0; x <= map_width; x++ ) {
    ctx.beginPath();

    // Start at the bottom of the map
    ctx.moveTo(x,map_height);

    // Get the y value for the x pixel on the sine wave
    var y = (map_height/2) - (Math.sin( (offset_x / map_width) / period ) * amplitude);

    offset_x++;

    // Draw the line up to the point on the sine wave
    ctx.lineTo(x,y);
    ctx.stroke();
}

Answer 1

If you want it to be physically accurate, you actually need to consider two offsets: a vertical (depending on the current date) and a horizontal (depending on the current time).

The horizontal offset X may be calculated by looking at the current time on some fixed geographic location on earth. The shadow offset will be 0 at midnight and will increase by 1/86400 for each seconds after midnight. So the formular is

offsetX = (curHour*3600 + curMinute*60 + curSeconds)/86400

The vertical offset will change between the Solstices on June 21st and Dec 22nd (if it's not a leap year, where the Solstices are on June 20th and Dec 21st). The maximum angles are 23.44° in both directions. We have 90° per hemisphere and 365/2 = 182.5 days between the two solstices, and we are working with a mapping of a circular motion, so a sin()-function has to be used. The wavelength of a sinus wave is 2pi, so we need pi for half the vertical offset Y of one year.

Please note, that I did not take leap seconds into account, so the calculation might be a bit off in the distant past/future.

// current time
$curHour = date("H");
$curMin  = date("i");
$curSec  = date("s");

// resulting offset X
$offsetX = ($curHour*3600 + $curMin*60 + $curSec)/86400;


echo "======== OFFSET X ==========\n";
echo "curHour:      $curHour\n";
echo "curMin:       $curMin\n";
echo "curSec:       $curSec\n";
echo "offsetX:      $offsetX\n\n";


// spring equinox date as day of year
$equinox = date("z", mktime(0, 0, 0, 3, 20));

// current day of year
    // first line is for testing purposes
//$curDay = date("z", mktime(0, 0, 0, 6, 21));
    $curDay = date("z");

// Day offset, mapped on the equinox offset
$offsetSin = ((365.25 - $equinox + $curDay)%365.25)/365.25;

// sinus wave offset
$offsetSinFactor = sin($offsetSin * 2 * pi());

// map onto angle
$offsetY = $offsetSinFactor * 23.44;

// optional: Mercator projection
$degreesPerRadian = 180.0 / pi();
$offsetYmercator = atan(sinh($offsetY)) * $degreesPerRadian;

// missing: mapping onto canvas height (it's currently
// mapped on $offsetY = 90 as the total height of the
// canvas.


echo "========= OFFSET Y =========\n";
echo "equinox day:  $equinox\n";
echo "curDay:       $curDay\n";
echo "offsetSin:    $offsetSin\n";
echo "offsetSinFac: $offsetSinFactor\n";
echo "offsetY:      $offsetY\n";
echo "offsetYmerc:  $offsetYmercator\n";

You should be able to port this calculation to any language you want.

Answer 2

You asked for more control over appearance

Check out the Geocommons JS api, which may be more suited to your purpose than Google Maps.

If you're going for flexibility, GEOS would be preferable to drawing the shape for a specific projection. I know there are php bindings, but I haven't used them myself.

Helmer Aslaksen has a great writeup about heavenly mathematics that should help you create an alorithm to draw a sunlit-area polygon using geos. You can test your code against measured sunrise/sunset times for accuracy.

Edit 1

Must be Python and Google Maps, you say?

• You can overlay data with KML,
• which will allow you to set a background color as you request.
• I like using the interactive sampler to get a feel for how things will look.
• Check out the Shapely package on pypi,
• and a discussion on how to create kml from it.