MollweideProjection.frink

View or download MollweideProjection.frink in plain text format


/** This program demonstrates the use of the country polygons in
    Countries.frink to draw a map of the world.  It can be readily altered to
    draw the map in your favorite projection.  This demonstrates the
    Mollweide projection:


    https://mathworld.wolfram.com/MollweideProjection.html
    https://en.wikipedia.org/wiki/Mollweide_projection
*/


use Country.frink

g = new graphics
g.stroke[0.001]
g.font["SansSerif", "bold", 1 degree]

// Iterate through all countries.
for [code, country] = Country.getCountryList[]
{
   cc = new color[randomFloat[0,1], randomFloat[0,1], randomFloat[0,1], .8]
   first = true
   for poly = country.borders  // Iterate through polygons in a country.
   {
      p = new filledPolygon   // This polygon is the filled country
      po = new polygon        // This is the outline of the country
      for [long, lat] = poly  // Iterate through points in polygon
      {
         [x,y] = latLongToXYMollweide[lat degree, long degree]
         p.addPoint[x, -y]
         po.addPoint[x, -y]
      }

      // Draw filled countries
      g.color[cc]
      g.add[p]

      // The polygons in Country.frink are sorted so the first polygon is the
      // largest. Just label the largest.
      if first
      {
         [cx, cy] = p.getCentroid[]
         g.color[0,0,0]
         g.text[code, cx, cy]
      }

      // Draw country outlines
      g.color[0.2, 0.2, 0.2, 0.8]
      g.add[po]

      first = false
   }
}

// Optional:  Draw the ellipse that makes the boundary
// You could do this first but with fillEllipse... to draw water
//g.color[0,0,0,.3]
//g.drawEllipseSides[-2 sqrt[2], -sqrt[2], 2 sqrt[2], sqrt[2]]

g.show[]
g.write["Mollweide.svg", 1000, undef]
g.write["Mollweide.png", 1000, undef]
g.write["Mollweide.html", 1000, 500]


/** Convert a latitude, longitude, and optional center longitude (long0)
    into x,y coordinates.

    The x coordinate ranges from -2 sqrt[2] to 2 sqrt[2]
    The y cooridnate ranges from -sqrt[2] to sqrt[2]]

    returns
     [x ,y]
*/

latLongToXYMollweide[lat, long, long0 = 0 degrees West] :=
{
   // The Mollweide projection uses an "auxiliary angle" theta where theta is
   // given by
   //   2 theta + sin(2 theta) = pi sin[lat]
   //
   // Where theta is found iteratively by iterating:
   // theta = theta-(2 theta + sin[2 theta] - pi sin[lat])/(2 + 2 cos[2 theta])
   //
   // or, by introducing a different supplementary angle theta1 where
   // theta = 1/2 theta1
   // The equations can be iterated as
   // theta1 = theta1 - (theta1 + sin[theta1] - pi sin[lat]) / (1 + cos[theta1])
   // with an initial guess
   // theta1 = 2 arcsin[2 lat / pi]
   theta1 = 2 arcsin[2 lat / pi]
   slat = sin[lat]

   do
   {
      ct = cos[theta1]
      if ct == -1
         break
      delta = - (theta1 + sin[theta1] - pi slat) / (1 + ct)
      theta1 = theta1 + delta
   } while abs[delta] > 1e-7  // error of 0.63 m on the earth
   
   theta = 1/2 theta1

   x = 2 sqrt[2]/pi (long - long0) cos[theta]
   y = sqrt[2] sin[theta]

   return [x,y]
}


View or download MollweideProjection.frink in plain text format


This is a program written in the programming language Frink.
For more information, view the Frink Documentation or see More Sample Frink Programs.

Alan Eliasen was born 18864 days, 11 hours, 24 minutes ago.