countryHeatMap.frink

View or download countryHeatMap.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

/** This is a simple data struture that holds the data for a country that we
    read from a file. */

class CountryData
{
   var iso2
   var value
   var name

   new[iso, displayName, val] :=
   {
      iso2 = iso
      name = displayName
      value = val
   }
}

filename = "countries.txt"

// This is a dictionary where the key is the 2-letter ISO code and the
// value is a CountryData object.
countryDict = new dict

// Read the data from the file and make a dictionary out of it.
largestValue = 0
for line = lines[filenameToURL[filename]]
{
   // Split line on semicolons
   [iso2, name, val] = split[%r/;/, line]
   val = eval[val]   // Turn string into a number
   if isUnit[val]
   {
      countryDict@iso2 = new CountryData[iso2, name, val]

      // Keep track of largest value for scaling
      if val > largestValue
         largestValue = val
   }
}


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

// Iterate through all countries' map data.
for [code, country] = Country.getCountryList[]
{
   cd = countryDict@code   // Fetch the CountryData from the dictionary
   if cd != undef
   {
      // Autoscale to largestValue
      red   = 1
      green = 1-(cd.value / largestValue)
      blue  = green
      cc = new color[red, green, blue]
   }
   
   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["heatmap.svg", 1000, undef]
g.write["heatmap.png", 1000, undef]
g.write["heatmap.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 countryHeatMap.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 18664 days, 19 hours, 45 minutes ago.