Download or view 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
}
Download or view 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 19963 days, 7 hours, 24 minutes ago.