The Finite Nearly-Conformal Normal Cylindrical Projection (Aasgaard's map projection)

The properties are:

• It maps the whole earth into a rectangle of finite size.
• Shape is approximately preserved over allmost the whole earth.
• The scale increases to a very large value towards the poles.
• The polar points (the north and south pole) are mapped to a line (the upper and lower edge of the rectangle).
• The meridians are mapped as vertical lines.
• The paralells are mapped as horizontal lines.
• The equations are very simple.

The projection was presented at the 10th International Symposium on Spatial Data Handling in Ottawa July 9-12 2002. It is publiched as: "Projecting a Regular Grid onto a Sphere or Ellipsoid", Rune Aasgaard, in: "Advances in Spatial Data Handling", Dianne Richardson and Peter van Oosterom (eds.), Springer-Verlag 2002, pp 339-350.

A pdf version of the article can be found here. A PowerPoint presentation used at the conference is found here (Warning: 16MB!).

Equations and program code

For most of my work I usually uses the parameters x_max = y_max = 2^30. When storing coordinates as integers that gives a minimum resolution of ca. 4 cm at the equator, better towards the poles.

//===========================================================================
// This goes in a .h file
//===========================================================================

#include <cmath>

/*
 * The Aasgaard Projection class definitions.
 */
class AasgaardProjection {
protected:
  double x_max_;
  double y_max_;

  double a_;
  double b_;
  double c_;
  double d_;

  double r_;
  double e2_;

  void adj_cyclic(double& x, const double& maxval) const {
    while (x < 0)
      x += maxval;
    while (x >= maxval)
      x -= maxval;
  }

  double diff_x(double x1, double x2) const {
    adj_cyclic(x1, x_max_);
    adj_cyclic(x2, x_max_);
    double dx = x1-x2;
    if (dx > x_max_/2.)
      dx = x_max_ - dx;
    if (dx < -x_max_/2.)
      dx = x_max_ + dx;
    return dx;
  }

  double W(const double& lat) const {
    double sin_lat = sin(lat);
    return sqrt(1.-e2_*sin_lat*sin_lat);
  }
  // Ellipsoidal normal radius of curvature
  double N(const double& lat) const {
    return r_/W(lat);
  }
  // Ellipsoidal meridial radius of curvature
  double M(const double& lat) const {
    double w = W(lat);
    return r_*(1.-e2_)/(w*w*w);
  }


  // Transformation of latitude
  double toLat(const double& y) const {
    return b_ == 0. ? y/a_ : tanh(y*b_/a_)/b_;
  }
  double toY(const double& lat) const {
    return b_ == 0. ? lat*a_ : log((1.+b_*lat)/(1.-b_*lat))*a_*.5/b_;
  }

public:
  AasgaardProjection(const double& x_max, const double& y_max,
		     const double& r = 6370997.0, const double& e2 = 0.);

  bool forward(double  lon, const double& lat,
	       double& x,         double& y) const;
  bool reverse(double  x,   const double& y,
	       double& lon,       double& lat) const;

  double dist(const double& x1, const double& y1,
	      const double& x2, const double& y2) const;
  double area(const double& x1, const double& y1,
	      const double& x2, const double& y2) const;


  double invAspect(const double& lat) const {
    double sin_lat = sin(lat);
    return
      (a_*(1.-e2_*sin_lat*sin_lat)*cos(lat)) /
      (c_*(1.-e2_)*(1.-(b_*b_*lat*lat)));
  }

  double getR()  const { return r_;  }
  double getE2() const { return e2_; }
};


//===========================================================================
// The rest goes in the .cc file
//===========================================================================
/* Constructor, takes the domain of the projection, 
 * 0 <= x < x_max and -y_max <= y <= y_max and 
 * an earth ellipsoid definition. 
 */
AasgaardProjection::AasgaardProjection(const double& x_max,
				       const double& y_max,
				       const double& r,
				       const double& e2)
{
  x_max_ = x_max;
  y_max_ = y_max;
  r_     = r;
  e2_    = e2;
  
  c_ = x_max/(2*M_PI);

  if (y_max*4. <= x_max) {
    // If we should use uniform latlon scaling instead.
    a_ = y_max/M_PI_2;
    d_ = 1.;
    b_ = 0.;
  }
  else {
    a_ = c_*(1.-e2_);
    // Compute d by fix-point iteration
    d_ = 0.;
    while (2+2==4) {
      double prev_d = d_;
      d_ = 2./(1.+exp(4.*(1.-d_)*y_max/(a_*M_PI)));
      if (fabs(d_-prev_d) < d_*1.0e-9)
	break;
    }
    b_ = 2*(1.-d_)/M_PI;
  }
}


//===========================================================================
/* Forward projection, from lat,lon to x,y */
bool AasgaardProjection::forward(double lon, const double& lat,
				 double& x, double& y) const
{
  if (lat < -M_PI/2. || lat > M_PI/2.)
    return false;

  adj_cyclic(lon, 2*M_PI);
  x = c_ * lon;
  y = toY(lat);

  return true;
}


//===========================================================================
/* Reverse projection, from x,y to lat,lon */
bool AasgaardProjection::reverse(double x, const double& y,
				 double& lon, double& lat) const
{
  if (y < -y_max_ || y > y_max_)
    return false;
  adj_cyclic(x, x_max_);
  
  lon = x/c_;
  lat = toLat(y);

  return true;
}


//===========================================================================
/* Simple computation of distance. */
double AasgaardProjection::dist(const double& x1, const double& y1,
				const double& x2, const double& y2) const
{
  if (y1 < -y_max_ || y1 > y_max_ ||
      y2 < -y_max_ || y2 > y_max_)
    return HUGE_VAL;
  double dy = y1-y2;
  double y0 = (y1+y2)/2.;
  double lat0 = toLat(y0);

  double inv_asp = invAspect(lat0);
  double m = M(lat0);
  double coshy0 = cosh(y0*b_/a_);

  double dx = diff_x(x1, x2);
  dx = dx*inv_asp;

  return m*sqrt(dy*dy+dx*dx)/(a_*coshy0*coshy0);
}


//===========================================================================
/* Simple computation of area */
double AasgaardProjection::area(const double& x1, const double& y1,
				const double& x2, const double& y2) const
{
  if (y1 < -y_max_ || y1 > y_max_ ||
      y2 < -y_max_ || y2 > y_max_)
    return HUGE_VAL;
  double dy = y1-y2;
  double y0 = (y1+y2)/2.;
  double lat0 = toLat(y0);

  double inv_asp = invAspect(lat0);
  double m = M(lat0);
  double coshy0 = cosh(y0*b_/a_);
  coshy0 *= coshy0;
  coshy0 *= coshy0;

  double dx = diff_x(x1, x2);
  dx = dx*inv_asp;

  return m*m*fabs(dx)*fabs(dy)/(a_*a_*coshy0);
}

Examples

The first level of quad tree adresses:

2	3
0	1

The second level of quad tree adresses:

22	23	32	33
20	21	30	31
02	03	12	13
00	01	10	11

The coordinate system uses the domain 0≤x<2^30 and -2^30≤y≤2^30. As the root of the quad tree covers -2^31≤x<2^31 and -2^31≤y<2^31 the data here starts at quad tree index 12 (southern hemisphere) and 30 (northern hemisphere.

Example code for generating a similar texture pyramid

#include <string>

//===========================================================================
/*
 * A pixel color (or an elevation value?)
 */
class Color {
public:

  // Adding up values
  Color& operator+=(const Color& c);

  // Scalar division
  Color  operator/(double v) const;

  // Not equal test
  bool operator!= (const Color& c) const;

  // The "undefined" value
  static Color undef();

  // The additive identity
  static Color zero();
};

//===========================================================================
/*
 * An image (or elevation grid?)
 */
class Image {
public:
  Image();
  Image(int width, int height);
  ~Image();

  Color getPixel(int col, int row) const;
  void  setPixel(int col, int row, const Color& c);

  bool load(std::string filename);
  bool save(std::string filename) const;
};

//===========================================================================
/*
 * A geo referenced image mosaic
 * This is really the difficult part!
 * In this example I use geographical coordinates in radians.
 */
class GeoImageMosaic {
public:
  // Find resolution of "best" image in the area
  // Return "false" if no image found
  bool getMaxResolution(const double& source_min_x,
			const double& source_min_y,
			const double& source_max_x,
			const double& source_max_y,
			double &rx, double &ry) const;

  // Compute pixel value for a geographical position
  Color calcPixel(const double& x, const double& y) const;
};

/*
 * This should go in a .h file
 */

//===========================================================================
/*
 * The real workhorse
 */
class Pyramid {
protected:
  AasgaardProjection proj_;
  
  // The source image (in lat/lon coordinates)
  GeoImageMosaic* source_;

  // The destination root
  std::string     dirname_;

  int             image_dim_;
  float           min_pixel_size_;
  bool            reuse_;

  Image* create(unsigned long long q_code, int q_level);

public:
  Pyramid(int image_dim = 128,
	  float min_pixel_size = 0.,
	  bool reuse = true);
  
  void   createWorld(GeoImageMosaic* source, std::string dirname);
};

//===========================================================================
/*
 * This should go in a .cc file
 */

#include <iostream>
#include <sstream>
#include <iomanip>
#include <unistd.h>
#include <sys/stat.h>
#include <sys/types.h>

using namespace std;

//===========================================================================
Pyramid::Pyramid(int image_dim, float min_pixel_size, bool reuse)
  : proj_((1<<30), (1<<30)), dirname_("./")
{
  source_         = NULL;
  min_pixel_size_ = min_pixel_size;
  image_dim_      = image_dim;
  reuse_          = reuse;
}

//===========================================================================
// Do the real work
void Pyramid::createWorld(GeoImageMosaic* source, std::string dirname)
{
  source_    = source;
  dirname_   = dirname;
  
  // Whole world
  create( 6ULL << 60, 2);
  create(12ULL << 60, 2);
}

//===========================================================================
Image* Pyramid::create(unsigned long long q_code, int q_level)
{
  // Create filename
  stringstream filename;
  filename << dirname_;
  for (int i = 1; i <= q_level; ++i) {
    filename << ((q_code >> (64-2*i)) & 0x3);
    if (i % 3 == 0 && i < q_level)
      filename << "/";
  }
  filename << ".jpg";
  string filestr = filename.str();

  // Check if file exists and is readable
  if (reuse_ && access(filestr.c_str(), R_OK) != -1) {
    // Read the file and return image
    cerr << "Reading : " << filestr << "\n";
    Image* image = new Image();
    if (image->load(filestr))
      return image;
  }

  // Compute image position
  unsigned int x = 0; 
  unsigned int y = 0;
  unsigned int width  = (1UL << (32 - q_level));
  unsigned int height = (1UL << (32 - q_level));
  
  for (int i = 1; i <= q_level; ++i) {
    x |= ((q_code >> (64-2*i)) & 0x1) << (32-i);
    y |= ((q_code >> (65-2*i)) & 0x1) << (32-i);
  }

  double min_x = (double)x - (double)(1UL << 31);
  double min_y = (double)y - (double)(1UL << 31);
  double source_min_x, source_min_y;
  proj_.reverse(min_x, min_y, source_min_x, source_min_y);
  double source_max_x, source_max_y;
  proj_.reverse(min_x+width, min_y+height, source_max_x, source_max_y);

  float pixel_size = (source_max_y - source_min_y) * proj_.getR() / image_dim_;

  if (source_max_x <= 0.000001)
    source_max_x = 2*M_PI;
  
  // Find resolution
  double rx, ry;
  if (!source_->getMaxResolution(source_min_x, source_min_y,
				 source_max_x, source_max_y,
				 rx, ry))
    return NULL;
  
  Image* image = new Image(image_dim_, image_dim_);

  if (rx > (source_max_x - source_min_x) / image_dim_ ||
      ry > (source_max_y - source_min_y) / image_dim_) {
    // Extract image from mosaic
    for (int row = 0; row < image_dim_; ++row)
      for (int col = 0; col < image_dim_; ++col) {
	double src_x, src_y, dst_x, dst_y;
	// Compute destination coordinates (in projected coordinate values)
	// Pixel surface centered...
	dst_x = min_x + (col+.5)*width /image_dim_;
	dst_y = min_y + (row+.5)*height/image_dim_;

	// Compute source coordinates (in radians)
	proj_.reverse(dst_x, dst_y, src_x, src_y);

	// Compute source "color" (or elevation value?)
	Color c = source_->calcPixel(src_x, src_y);
	image->setPixel(col, row, c);
      }
  }
  else {
    // Assemble images from children
    Image* img[2][2];
    img[0][0] = create(q_code, q_level+1);
    img[1][0] = create(q_code + (1ULL << (62-2*q_level)), q_level+1);
    img[0][1] = create(q_code + (2ULL << (62-2*q_level)), q_level+1);
    img[1][1] = create(q_code + (3ULL << (62-2*q_level)), q_level+1);
    
    for (int imr = 0; imr < 2; ++imr)
      for (int imc = 0; imc < 2; ++imc) 
	if (img[imc][imr] != NULL)
	  for (int row = 0; row < image_dim_; row += 2)
	    for (int col = 0; col < image_dim_; col += 2) {
	      Color pixel, sum_pixel = Color::zero();
	      int sum = 0;
	      Image* im = img[imc][imr];
	      pixel = im->getPixel(col,   row);
	      if (pixel != Color::undef()) {
		sum_pixel += pixel;
		++sum;
	      }
	      pixel = im->getPixel(col+1, row);
	      if (pixel != Color::undef()) {
		sum_pixel += pixel;
		++sum;
	      }
	      pixel = im->getPixel(col,   row+1);
	      if (pixel != Color::undef()) {
		sum_pixel += pixel;
		++sum;
	      }
	      pixel = im->getPixel(col+1, row+1);
	      if (pixel != Color::undef()) {
		sum_pixel += pixel;
		++sum;
	      }
	      image->setPixel((image_dim_*imc + col)/2,
			      (image_dim_*imr + row)/2,
			      sum > 0 ? pixel/sum : Color::undef());
	    }
    delete img[0][0];
    delete img[0][1];
    delete img[1][0];
    delete img[1][1];
  }
  
  // Save only if we are above the accepted minimum detail level
  if (pixel_size > min_pixel_size_) {
    // Save image
    cout << "Saving : " << filestr << "\n";
    
    string::size_type dir_ix = filestr.find('/');
    while (dir_ix != string::npos) {
      string dir_str(filestr, 0, dir_ix);
      mkdir(dir_str.c_str(),
	    S_IRUSR | S_IWUSR | S_IXUSR |
	    S_IRGRP | S_IXGRP | S_IWGRP |
	  S_IROTH | S_IXOTH);
      dir_ix = filestr.find('/', dir_ix+1);
    }
    
    image->save(filestr);
  }
  return image;
}

int main(int argc, char* argv[])
{
  GeoImageMosaic source;

  // Do all of the loading of source files
  // source->addLoad("World-data-set");
  
  // Create image hierarchy
  Pyramid pyramid(256);
  pyramid.createWorld(&source, "textures/");

  exit(0);
}