mapStatesSinusoidal.frink

View or download mapStatesSinusoidal.frink in plain text format


/** This program draws a map of U.S. states using a sinusoidal projection
    centered on the longitudinal centroid of the continental U.S.

    It requires the GeoJSON file from:

    https://web.archive.org/web/20130615162524if_/http://eric.clst.org:80/wupl/Stuff/gz_2010_us_040_00_500k.json

    The GeoJSON format is defined in RFC 7946:
    https://tools.ietf.org/html/rfc7946
*/


// This allows us to attach visualvm for profiling
//input=input["Press start", ""]

// This is the center longitude which is approximately the center longitude
// of the contiguous 48 states.
centerLong := -98

s1 = now[]
s = now[]
println["Reading file."]

// The location you've extracted the file to:
f = read["file:/home/eliasen/prog/frink/samples/gz_2010_us_040_00_500k.json"]
e = now[]
println["Finished in " + format[e-s, "s", 1] + "\n"]

s = now[]
b = parseJSON[f]
println["JSON parsed."]
e = now[]
println["Finished in " + format[e-s, "s", 1] + "\n"]

// The "type" key indicates what the top-level container is, hopefully a
// FeatureCollection
filetype = b@"type"
println["This file contains a $filetype"]

g = new graphics
g.font["SansSerif", 2]

if filetype == "FeatureCollection"
   plotFeatureCollection[b, g]
else
   if filetype == "Feature"
      plotFeature[b, g]


println["Writing svg"]
   
s=now[]   
g.write["statesSinusoidal.svg", 1920, undef]
e=now[]
println["Finished in " + format[e-s, "s", 1] + "\n"]
   
println["Writing svgz"]
s=now[]
g.write["statesSinusoidal.svgz", 1920, undef]
e=now[]
println["Finished in " + format[e-s, "s", 1] + "\n"]

println["Writing HTML5"]
s=now[]
g.write["statesSinusoidal.html", 1920, undef]
e=now[]
println["Finished in " + format[e-s, "s", 1] + "\n"]

println["Writing png"]
s=now[]
g.write["statesSinusoidal.png", 1920, undef]
e=now[]
println["Finished in " + format[e-s, "s", 1] + "\n"]

e1 = now[]
println["Total time: " + format[e1-s1, "s", 1] + "\n"]   

g.show[]
g.invertGrays[].show[]   


/** This plots a FeatureCollection, whose "features" key should contain an
    array of Feature objects. */

plotFeatureCollection[fc, g is graphics] :=
{
   for feature = fc@"features"
      plotFeature[feature, g]
}

/** This plots a Feature object.  A Feature has a "properties" key that
    contains a dictionary that describes stuff (in this case, the "tzid" field
    may be the only member,) and a key called "geometry" which contains the
    object (in this case, a Polygon or MultiPolygon.)
*/
 
plotFeature[feature, g is graphics] :=
{
   for [key, value] = feature
   {
      print["$key:\t"]
      if key == "type"   // This is a string, hopefully "Feature"
         print[value]
      
      if key == "properties" // This contains a dictionary of key-value
         // pairs.  We most likely want the "NAME" pair
      {
         print["NAME: " + value@"NAME" + "\t"]
         print["STATE: " + value@"STATE" + "\t"]
      }

      g.color[randomFloat[0,1], randomFloat[0,1], randomFloat[0,1]]
      // A Feature has a key called "geometry" which describes the type
      if key == "geometry"
      {
         type = value@"type"
         print[type]
         if (type == "Polygon")
         {
            coordinates = value@"coordinates"
            g.add[makePolygon[coordinates]]
         } else
         if (type == "MultiPolygon")
         {
            coordinates = value@"coordinates"
            g.add[makeMultiPolygon[coordinates]]
         } else
         println["Unknown type $type"]
      }

      println[]
   }
   println[]
}

/** Make a "polygon", given an object containing a GeoJSON coordinates array.
    In the GeoJSON specification, a "polygon" may actually be an array of
    concentric disconnected polygons with the first one being a surrounding
    polygon and the latter ones being "holes" in this object.  If there are
    no holes, then this is returned as a Polygon object, otherwise as a
    GeneralPath.
*/

makePolygon[coordinates] :=
{
   length = length[coordinates]

   // If length == 1, this can be a simple polygon with no holes
   if length == 1
   {
      ret = new filledPolygon
      for [x,y] = coordinates@0
      {
         // This is a kludge to get the couple of Aleutian islands that are
         // west of 180 degrees west to not make the whole map wrap around
         if x > 0
            x = x - 360
         ret.addPoint[(x-centerLong) cos[y degrees], -y]
      }
      return ret
   }

   // Otherwise, this is a complex GeneralPath with holes
   ret = new filledGeneralPath
   outer = coordinates@0
   for [x, y] = outer   // Draw the outer polygon
   {
      // This is a kludge to get the couple of Aleutian islands that are
      // west of 180 degrees west to not make the whole map wrap around
      if x > 0
         x = x - 360
      ret.addPoint[(x-centerLong) cos[y degrees], -y]
   }

   ret.close[]

   for inner = slice[coordinates, 1, undef]  // Draw all the inner polygons
   {
      for [x, y] = inner
      {
         // This is a kludge to get the couple of Aleutian islands that are
         // west of 180 degrees west to not make the whole map wrap around
         if x > 0
            x = x - 360
         ret.addPoint[(x-centerLong) cos[y degrees], -y]
      }

      ret.close[]
   }

   return ret
}

// This draws a set of polygons and returns it as a graphics.
makeMultiPolygon[coordinates] :=
{
   g = new graphics
   for polygon = coordinates
      g.add[makePolygon[polygon]]

   return g
}


View or download mapStatesSinusoidal.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 18492 days, 22 hours, 34 minutes ago.