/** 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] }