• Derivation of the Conversion:
  • \(\cos \phi = d / R \Rightarrow d = R \cos \phi\)
  • \(\cos \lambda = p_{z}/d \Rightarrow p_{z} = d \cos \lambda = R\cos \phi \cos \lambda\)
  • \(\sin \lambda = p_{x}/d \Rightarrow p_{x} = d \sin \lambda = R \cos \phi \sin \lambda\)
  • \(p_{y} = R \sin \phi\)
• Note:
  • This is very similar to the conversion of spherical coordinates to Cartesian coordinates (but it uses the latitude not the colatitude)

svaexamples/cartography/projections.c (Fragment: lonlat2xyz)

/**
 * Convert from spherical longitude/latitude to 
 * 3D rectangular coordinates (with the origin at the 
 * center of the sphere)
 *
 * @param lon   The longitude in degrees (IN)
 * @param lat   The latitude  in degrees (IN)
 * @param x     The x-coordinate (OUT)
 * @param y     The y-coordinate (OUT)
 * @param x     The z-coordinate (OUT)
 */
void lonlat2xyz(double lon, double lat, double &x, double &y, double &z)
{
   double     coslat, coslon, sinlat, sinlon;   

   coslat = cos(lat*RADIANS_PER_DEGREE);
   coslon = cos(lon*RADIANS_PER_DEGREE);
   sinlat = sin(lat*RADIANS_PER_DEGREE);
   sinlon = sin(lon*RADIANS_PER_DEGREE);
   
   x = (coslat * sinlon);
   y = (sinlat);
   z = (coslat * coslon);
}
        

• What We Need To Do:
  • Represent geographic features
  • Read country boundaries from a file
  • Create longitude/latitude lines
  • Draw geographic features
• 3D Visualization:
  • We will use an orthographic visualization (i.e., we will essentially ignore \(z\) coordinates)

svaexamples/cartography/feature.c

#ifndef FEATURE_C
#define FEATURE_C


/**
 * A generic geographic feature (e.g., a polygon, piecewise linear curve)
 * 
 * @author  Prof. David Bernstein, James Madison University
 * @version 1.0
 */

struct feature
{
    int               id;      // Unique identifier 
    char              *type;   // POLYGON, LINE_STRIP, ...  
    int               size;    // Number of vertices

    double            *lon;    // Longitudes
    double            *lat;    // Latitudes

    double            *x;      // Cartesian x-coordinate
    double            *y;      // Cartesian y-coordinate
    double            *z;      // Cartesian z-coordinate

    double            *xProj;  // Projected x-coordinate (for maps)
    double            *yProj;  // Prohected y-coordinate (for maps)
};

#endif
        

• File Format:
  • One record description of each file (number of polygons, number of countries)
  • One record description of each polygon (polygon number, country number,number of vertices,country name)
  • Each polygon has one record for each vertex (longitude, latitude)
• An Important Observation:
  • Each country might consist of several polygons

svaexamples/cartography/initialization.c (Fragment: readPolygons)

/**
 * Create the array of features, read the countries from 
 * the file, create the polygons, and insert them 
 * into the array of features
 *
 * @return   The index of the next empty feature
 */
int readPolygons()
{
   char        name[80];   
   double      lat, lon;   
   feature     *f;
   FILE        *in;
   int         countryNumber, k, numberOfCountries;
   int         numberOfPolygons, polygonNumber, size;
   


   in = fopen("world.txt", "r");
   
   fscanf(in, "%d,%d", &numberOfPolygons, &numberOfCountries);
   
   numberOfFeatures = numberOfPolygons + GRID_SIZE;
   
   features = new feature[numberOfFeatures];
   
   for (k=0; k<numberOfPolygons; k++)
   {
      fscanf(in, "%d,%d,%d,%s", &polygonNumber, &countryNumber, &size, name);
      
      f       = new feature;      
      initializeFeature(f, size);
      f->id   = countryNumber;
      f->type = "POLYGON";   
      f->size = size;      
      for (int i=0; i<size; i++)
      {
         fscanf(in, "%lf,%lf", &lon, &lat);
         setCoordinate(f, i, lon, lat);
      }
      features[k] = *f;
   }
   fclose(in);  


   return k;
}
        

svaexamples/cartography/initialization.c (Fragment: createGrid)

/**
 * Create the lines that comprise the longitude/latitude 
 * grid and insert them into the array of features
 *
 * @param firstIndex  The index to use for the first line
 */
void createGrid(int firstIndex)
{
   feature     *f;
   int         k, size;


   k = firstIndex;

   // Create the longitude lines
   for (int lon=-180; lon<=180; lon+=10) // 37 lines
   {
      size = 19;      
      f = new feature;      
      initializeFeature(f, size);
      f->id   = 214;      
      f->type = "LINE_STRIP";
      f->size = size;
      
      int i = 0;      
      for (int lat=-90; lat<=90; lat+=10)
      {
         setCoordinate(f, i, (double)lon, (double)lat);
         i++;         
      }
      features[k] = *f;      
      k++;      
   }


   // Create the latitude lines
   for (int lat=-90; lat<=90; lat+=10) // 19 lines
   {
      size = 37;      
      f = new feature;      
      initializeFeature(f, size);
      f->id   = 214;      
      f->type = "LINE_STRIP";
      f->size = size;      

      int i = 0;      
      for (int lon=-180; lon<=180; lon+=10)
      {
         setCoordinate(f, i, (double)lon, (double)lat);
         i++;         
      }          
      features[k] = *f;      
      k++;      
   }
}
        

svaexamples/cartography/globe.c (Fragment: display)

/**
 * Display the Earth
 */
void onDisplay()
{
   feature                             *f;
   GLenum                              type;   
   GLfloat                             v[3];


   for (int i=0; i<numberOfFeatures; i++)
   {
      // Get the Feature
      f = &features[i];

      // Determine the Feature's type
      if      (strcmp(f->type, "POLYGON") == 0)    type = GL_POLYGON;
      else if (strcmp(f->type, "LINE_STRIP") == 0) type = GL_LINE_STRIP;
      else                                         type = GL_POINTS;

      // Begin the graphics primitive
      glBegin(type);
      {
         glColor3fv(colors[f->id]);

         for (int j=0; j<f->size; j++)
         {
            // Use the 3D rectangular coordinates
            v[0] = f->x[j];
            v[1] = f->y[j];
            v[2] = f->z[j];

            glVertex3fv(v);
         }
      }
      glEnd();
   }
}
        

svaexamples/cartography/globe.c (Fragment: main)

/**
 * The entry point of the application.
 *
 * @param argc  The number of command line arguments
 * @param argv  The array of command line arguments
 * @return      A status code
 */
int main(int argc, char **argv)
{
   float min[] = {-3.0, -3.0, -3.0};
   float max[] = { 3.0,  3.0,  3.0};
   int   k;   

   
   // Create the array of colors
   initializeColors();   

   // Create the Feature set
   k = readPolygons();
   createGrid(k);   
   
   // Create and initialize the SVA window
   svaCreate(argc, argv, min, max);
   svaShowAxes(false);
   
   // Start the SVA visualization
   svaStartVisualization();
}
        

• Objective:
  • Project points on the surface of the Earth onto a map
• Some Intuition:
  • Shine a light onto or through a transparent Earth ` and capture the shadows cast by the opaque features
  • The parameters are: the shape of the screen (called the projection surface), the position of the projection surface, and the location of the light source

• The Most Common:
  • Conformal (i.e., angles are preserved)
  • Equal Area (i.e., areas are in constant proportion)
  • Equidistant (i.e., distances are in constant proportion)
• An Important Mathematical Result:
  • A single projection can not be both conformal and equal area

• Parameters:
  • \(\lambda_{0}\) is the standard longitude (i.e., the horizontal center of the projection)
• Projection:
  • \(x = R (\lambda - \lambda_{0}) \)
  • \(y = R \sin(\phi)\)
• Inverse:
  • \(\lambda = \lambda_{0} + \frac{x}{R}\)
  • \(\phi = \sin^{-1}(y / R)\)

svaexamples/cartography/projections.c (Fragment: cylindricalEqualArea)

/**
 * Project from spherical longitude/latitude to 
 * 2D rectangular coordinates (with the origin at the 
 * equator and Greenwich meridian)
 *
 * @param lon   The longitude in degrees (IN)
 * @param lat   The latitude  in degrees (IN)
 * @param x     The x-coordinate (OUT)
 * @param y     The y-coordinate (OUT)
 */
void cylindricalEqualArea(double lon, double lat, double &x, double &y)
{
   x = lon*RADIANS_PER_DEGREE; 
   y = sin(lat*RADIANS_PER_DEGREE);   
}
        

• What We Need To Do:
  • Everything for visualizing a globe
  • Perform the projection
• An Important Observation:
  • Some projections will "tear" polygons

svaexamples/cartography/map.c (Fragment: display)

/**
 * Display the Earth
 */
void onDisplay()
{
   feature                             *f;
   GLenum                              type;   
   GLfloat                             v[3];

   for (int i=0; i<numberOfFeatures; i++)
   {
      // Get the Feature
      f = &features[i];

      // Determine the Feature's type
      if      (strcmp(f->type, "POLYGON") == 0)    type = GL_POLYGON;
      else if (strcmp(f->type, "LINE_STRIP") == 0) type = GL_LINE_STRIP;
      else                                         type = GL_POINTS;

      // Begin the graphics primitive
      glBegin(type);
      {
         glColor3fv(colors[f->id]);

         for (int j=0; j<f->size; j++)
         {
            // Use the projected coordinates
            v[0] = f->xProj[j];
            v[1] = f->yProj[j];
            v[2] = 0.0;

            glVertex3fv(v);
         }
      }
      glEnd();
   }
}