Download or view MollweideProjection.frink in plain text format
/** This program demonstrates the use of the country polygons in
Country.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]
}
Download or view 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 19966 days, 15 hours, 17 minutes ago.