Seasonality displays a IR composite satellite image that covers the whole globe. Unfortunately, during the past few weeks, the previous source of this satellite image ceased to update. It was time to write my own program to create a composite image.

The toughest part of this process (at least for us non-GIS types) is handling the geometry involved in converting map projections. It was extremely difficult for me to piece all of this together into something that works, as no one has posted any kind of solution online. I thought I might save some other people the trouble, and post the solution.

To create a composite, the basic algorithm is to iterate over a map of your chosen projection and to query the various satellite image sources for a color value for the latitude/longitude point you are currently generating. Once you have a color, you can simply set that color for the pixel that it applies to in your new map. I use a mercator/cylindrical style projection in the final map, with equal pixel distance for a change in latitude or longitude. This results in an image that is twice as wide as it is tall, with 360 degrees all the way around the equator along the X axis, and 180 degrees from the north pole to the south as the Y axis.

The tough part is that query to the source satellite image. The projection of these satellite images has a few different names. Officially it’s a Vertical Perspective Near-Side projection, but it’s also known as the GEOS projection. I boiled it down into a request:

float getPixelForLatLon(float latitude, float longitude)

Satellite images are 1-channel grayscale, so just returning a single float with the gray value works fine here. Now before we can understand what is going on, here are some of the rules and assumptions taking place here. First, the satellite is in geosynchronous orbit, at an elevation of 35,786 kilometers above the surface of the earth. The satellite is above the equator (so latitude=0) and the longitude varies based on the portion of the earth being covered.

With that in mind, let’s look at the pseudocode for this function.

// First convert from degrees to radians.
float latitudeRad = degrees2Radians(latitude);
float longitudeRad = degrees2Radians(longitude);
float satelliteLonRad = degrees2Radians(satelliteLongitude);

// Now for a bunch of trig equations.

// c_lat is the corrected angle the given latitude is above the equator
// after taking into account the ellipsoid shape of the earth.
// Here, 0.993243 is the polar radius (6356.5838 km)
// divided by the equator radius (6378.1690 km).
float c_lat = atan(0.993243 * tan(latitudeRad));

// r_L is the radius from the center of the earth to the given point.
float r_L = 6356.5838 / sqrt(1 - (0.00675701 * cos(c_lat) * cos(c_lat)));

// r_[123] are xyz components of the vector from the satellite
// to the point on earth.
float r_1 = 42164 - (r_L * cos(c_lat) * cos(longitudeRad - satelliteLonRad));
float r_2 = -1 * r_L * cos(c_lat) * sin(longitudeRad - satelliteLonRad);
float r_3 = r_L * sin(c_lat);

// r_N is the distance from the satellite to the point on earth.
float r_N = sqrt(r_1 ** 2 + r_2 ** 2 + r_3 ** 2);

// Finally we get to an intermediate x, y coordinate.
// These are NOT x, y coordinates in the image.  Rather,
// they are angles (in radians) that the satellite camera
// would have to rotate off-center to see that lat/lon
// location.  If x = y = 0, then the camera is looking directly
// at the earth point underneath it.  Positive X rotates to the
// right, positive Y rotates down towards the south.
float int_x = atan2(-r_2, r_1);  // Equivalent: atan(-r_2 / r_1)
float int_y = asin(-r_3 / r_N);

// We got to this point without knowing anything about our
// actual satellite image, now it's time to look at the
// full disc satellite image we are working with.
// Note that coordinates x = y = 0 in the upper-left corner
// of the image.

// offset_x is the number of pixels that pad the disc on the left.
// offset_y is the number of pixels that pad the disc on the top.
// x_radius is the radius of the disc in the horizontal direction.
// y_radius is the radius of the disc in the vertical direction.
float x = offset_x + x_radius + (int_x * 6.6073212476 * x_radius);
float y = offset_y + y_radius + (int_y * 6.6073212476 * y_radius);

// Here, get the pixel in the satellite at position x, y, and return it.

That’s it! You now have a color from the x, y position in the full disc satellite image, and you can paint that in the final satellite composite. Note that the satellite composite will not have data on the poles, because the satellite can only see to a latitude of 81.3 degrees. You can see sample output of a composite image using this method (overlaid on a standard earth map). To find more information about how this works, check out the WMO document number C1-SP-810-004.