Back to main page
Back to draw contents
This section shows how to plot cartographic maps with the draw package. The graphic object to be used is geomap; its syntax is:
Some of the examples below can take a while to be rendered.
A very simple map defined by hand. All these plots could be made by means of the graphic object points (with option points_joined set to true) and lists instead of arrays. But problems arise, in terms of processing time, when working with large sets of coordinates to draw real maps.
/* Vertices of boundary #0: {(1,1), (2,5), (4,3)} */
bnd0: make_array(flonum,6)$
( bnd0[0]:1.0, bnd0[1]:1.0, bnd0[2]:2.0,
bnd0[3]:5.0, bnd0[4]:4.0, bnd0[5]:3.0 )$
/* Vertices of boundary #1: {(4,3), (5,4), (6,4), (5,1)} */
bnd1: make_array(flonum,8)$
( bnd1[0]:4.0, bnd1[1]:3.0, bnd1[2]:5.0, bnd1[3]:4.0,
bnd1[4]:6.0, bnd1[5]:4.0, bnd1[6]:5.0, bnd1[7]:1.0)$
/* Vertices of boundary #2: {(5,1), (3,0), (1,1)} */
bnd2: make_array(flonum,6)$
( bnd2[0]:5.0, bnd2[1]:1.0, bnd2[2]:3.0,
bnd2[3]:0.0, bnd2[4]:1.0, bnd2[5]:1.0 )$
/* Vertices of boundary #3: {(1,1), (4,3)} */
bnd3: make_array(flonum,4)$
( bnd3[0]:1.0, bnd3[1]:1.0, bnd3[2]:4.0, bnd3[3]:3.0)$
/* Vertices of boundary #4: {(4,3), (5,1)} */
bnd4: make_array(flonum,4)$
( bnd4[0]:4.0, bnd4[1]:3.0, bnd4[2]:5.0, bnd4[3]:1.0)$
/* Pack all together in boundaries_array */
boundaries_array: make_array(any,5)$
(boundaries_array[0]: bnd0, boundaries_array[1]: bnd1,
boundaries_array[2]: bnd2, boundaries_array[3]: bnd3,
boundaries_array[4]: bnd4 )$
draw2d(geomap([0,1,2,3,4]))$
/* Highlighting borders */
draw2d(geomap([0,1,2,3,4]),
line_width = 3,
color = blue,
geomap([2,3,4]))$
/* Country names */
draw2d(geomap([0,1,2,3,4]),
label(["Country-1",2.5,3],
["Country-2",5,2.8],
["Country-2",3,1]))$
Most of the time one is interested in plotting real maps. Package worldmap contains coast lines and boundaries coordinates in (longitude, latitude) format in region [-180,180]×[-90,90]. When the package is loaded, global variable boundaries_array is set to an array of length 2801, each element being another array containing the coordinates of a boundary. See the following results:
(%i1) /* Load package */
load(worldmap)$
(%i2) /* Now boundaries_array contains
2801 boundaries, from 0 to 2800 */
arrayinfo(boundaries_array);
(%o2) [declared, 1, [2800]]
(%i3) /* Info on boundary 1351:
it contains 310 elements,
or 155 pairs (longitude, latitude) */
( bnd: boundaries_array[1351],
arrayinfo(bnd) );
(%o3) [declared, 1, [309]]
(%i4) /* The latest boundary in boundaries_array
contains three pairs (longitude, latitude) */
boundaries_array[2800];
(%o4) {Array: #(-180.0 -90.0 180.0 -90.0 179.99998 -78.301682)}
In the next examples it is assumed that you executed load(worldmap) first; package draw does not automatically load worldmap. With those maps involving a great number of boundaries the user must be patient, since gnuplot needs to process large files, and that takes a while.
Let's now plot boundaries from 1000 to 2000.
draw2d(geomap(makelist(k,k,1000,2000)));
The complete map. Running clisp, Maxima needed 8.627688 seconds to execute the following command (function time), but you have to wait some more time due to the gnuplot process (file data.gnuplot is about 1.7 MB).
draw2d(geomap(makelist(k,k,0,2800)));
Not only boarders and coastlines are defined in worldmap, but also countries and continents.
c1: gr2d(geomap([Canada,United_States,Mexico,Cuba]))$
c2: gr2d(geomap([Africa]))$
c3: gr2d(geomap([Oceania,China,Japan]))$
c4: gr2d(geomap([France,Portugal,Spain,Morocco,Western_Sahara]))$
draw(terminal = eps_color,
columns = 2,
c1,c2,c3,c4);
There are in worldmap some utilities of interest, especially for defining new geographical objects. Sicily (south Italy) is not defined in worldmap, but we know that it is in the rectangular region defined by vertices 10.4E 41.5N and 20.7E 35.4N (just tell Maxima to plot Italy and move the mouse around). Now call function region_boundaries to know which vertices are within this region and plot them:
(%i10) z: region_boundaries(10.4,41.5,20.7,35.4); (%o10) [1846, 1863, 1864, 1881, 1888, 1894] (%i11) draw2d(geomap(z))$
With this output, you don't know how to define Sicily as a Maxima variable, but you can call function numbered_boundaries, also defined in worldmap, to see the numbers of the polygonal segments:
numbered_boundaries(z)$
We know now that Sicily is formed by boundary number 1881. Let's make use of it.
Sicily: [1881]$
/* cities coordinates and populations from wikipedia */
labs: label(["Agrigento",13.583,37.317 - 0.1],
["Caltanisseta",14.067,37.4833 - 0.1],
["Catania",15.066,37.517 - 0.1],
["Enna",14.276,37.5633 - 0.1],
["Messina",15.55,38.183 - 0.1],
["Palermo",13.367,38.117 - 0.1],
["Ragusa",14.75,36.933 - 0.1],
["Siracusa",15.283,37.083 - 0.1],
["Trapani",12.517,38.017 - 0.1] )$
pop1: points([[13.583,37.317], /* Agrigento */
[14.067,37.4833], /* Caltanisseta */
[14.276,37.5633], /* Enna */
[14.75,36.933], /* Ragusa */
[12.517,38.017]])$ /* Trapani */
pop2: points([[15.55,38.183], /* Messina */
[15.283,37.083]])$ /* Siracusa */
pop3: points([[15.066,37.517]])$ /* Catania */
pop4: points([[13.367,38.117]])$ /* Palermo */
draw2d(yrange = [36.5,39],
xrange = [12.2,15.7],
title = "Sicilian cities populations",
line_width = 3,
geomap(Sicily),
point_type = filled_circle,
color = green,
key = " < 100000", point_size = 2, pop1,
key = "100000 - 300000", point_size = 3, pop2,
key = "300000 - 500000", point_size = 4, pop3,
key = " > 500000", point_size = 5, pop4,
color = blue,
labs )$
Playing together with image and geomap. Variable European_Union is defined in worldmap. You can download the flags and the mermaid from here.
draw2d(
terminal = eps_color,
/* background */
border=false,
fill_color=cyan,
rectangle([-28.80,27.53],[34.79,70.26]),
/* boundaries */
line_width = 3,
color = "dark-violet",
geomap(European_Union),
/* flags */
image(read_xpm("austria.xpm"),13.3,46.9,2,1.2),
image(read_xpm("belgium.xpm"),3.66,49.77,2,1.2),
image(read_xpm("bulgaria.xpm"),24,42,2,1.2),
image(read_xpm("cyprus.xpm"),31,33,2,1.2),
image(read_xpm("czech_republic.xpm"),14.18,49.15,2,1.2),
image(read_xpm("denmark.xpm"),6.86,55.46,2,1.2),
image(read_xpm("estonia.xpm"),27.02,58,2,1.2),
image(read_xpm("finland.xpm"),26.03,63.5,2,1.2),
image(read_xpm("france.xpm"),1.3,46.32,2,1.2),
image(read_xpm("germany.xpm"),8,50.5,2,1.2),
image(read_xpm("greece.xpm"),19.86,36.37,2,1.2),
image(read_xpm("hungary.xpm"),17.65,45.79,2,1.2),
image(read_xpm("ireland.xpm"),-9.69,50.24,2,1.2),
image(read_xpm("italy.xpm"),12.39,40.82,2,1.2),
image(read_xpm("latvia.xpm"),27.71,55.22,2,1.2),
image(read_xpm("lithuania.xpm"),24.89,53.59,2,1.2),
image(read_xpm("luxembourg.xpm"),6.14,48.90,2,1.2),
image(read_xpm("malta.xpm"),13.43,33.69,2,1.2),
image(read_xpm("netherlands.xpm"),4.01,52.88,2,1.2),
image(read_xpm("poland.xpm"),17.86,51.73,2,1.2),
image(read_xpm("portugal.xpm"),-11.56,39.48,2,1.2),
image(read_xpm("romania.xpm"),23.85,45.18,2,1.2),
image(read_xpm("slovakia.xpm"),19.37,48.72,2,1.2),
image(read_xpm("slovenia.xpm"),14.16,44.65,2,1.2),
image(read_xpm("spain.xpm"),-4.79,39.82,2,1.2),
image(read_xpm("sweden.xpm"),14.05,62.63,2,1.2),
image(read_xpm("united_kingdom.xpm"),-2.61,52.11,2,1.2),
/* the mermaid */
image(read_xpm("sirena.xpm"),-26,44,9,9) )$
With function make_poly_country in package worldmap it is possible to build polygonal objects from country names, so that they can be filled with different colours. The return value of make_poly_country is a list of polygons, since some countries need more than just one polygonal because of the islands. A simple example follows.
mymap: append([ color = white], /* borders */
[ fill_color = red], make_poly_country(Bolivia),
[ fill_color = cyan], make_poly_country(Paraguay),
[ fill_color = green], make_poly_country(Colombia),
[ fill_color = blue], make_poly_country(Chile),
[ fill_color = "#23ab0f"], make_poly_country(Brazil),
[ fill_color = goldenrod], make_poly_country(Argentina),
[ fill_color = "midnight-blue"], make_poly_country(Uruguay))$
apply(draw2d, mymap)$
There is also make_poly_continent to build polygonal objects from continent names or lists of countries.
mymap2: append([ grid = true,
color = white /* borders */],
[ fill_color = red], make_poly_continent(North_America),
[ fill_color = green], make_poly_continent(Central_America),
[ fill_color = blue], make_poly_continent(South_America))$
apply(draw2d, mymap2)$
In 3d, maps are projected on the unit sphere (center (0,0,0) and radius 1).
draw3d(line_width=3,
color = blue, geomap(United_States),
color = red, geomap(Spain),
color = "light-goldenrod", geomap(Germany),
color = coral, geomap(Russia),
color = "#58a30f", geomap(Portugal),
color = magenta, geomap(Brazil),
color = orange, geomap(Japan),
color = purple, geomap(Australia) )$
World coastlines and islands projected on the sphere. Variable World_coastlines is defined in package worldmap.
c:%pi/180$
draw3d(surface_hide=true,
color =yellow,
parametric_surface(
cos(phi*c)*cos(theta*c),
cos(phi*c)*sin(theta*c),
sin(phi*c),
theta,-180,180,
phi,-90,90),
color=blue,
line_width=3,
geomap(World_coastlines))$
Australia projected on two different spheres.
draw3d(color = red,
geomap(Australia), /* projected on unit sphere */
color = blue,
line_width = 3,
geomap([Australia],
[spherical_projection,2,2,2,3]))$
Cylindrical and spherical projections.
draw3d(color=red,
geomap([America_coastlines,Eurasia_coastlines],
[spherical_projection,2,2,2,3]),
color=blue,
geomap([America_coastlines,Eurasia_coastlines],
[cylindrical_projection,2,2,2,3,4]));
Another example of cylindrical and spherical projections. The sphere is also plotted.
draw3d(surface_hide=true,
rot_horizontal = 60,
rot_vertical = 131,
color=yellow,
parametric_surface(
cos(phi)*cos(theta),
cos(phi)*sin(theta),
sin(phi),
theta,-%pi,%pi,
phi,-%pi/2,%pi/2),
color=red,
geomap([South_America,Africa,Australia],
[spherical_projection,0,0,0,1]),
color=blue,
geomap([South_America,Africa,Australia],
[cylindrical_projection,0,0,0,1,2]))$
Conic projection.
a:90$
r:1$
[cx,cy,cz]:[0,0,0]$
c: float(%pi/180)$
h:r/sin((a/2) * c)$
R:r/cos((a/2) * c)$
m: sqrt(h^2+R^2)$
draw3d(color=yellow,
/*surface_hide=true,*/
parametric_surface(cx+R*cos(theta)*t,
cy+R*sin(theta)*t,
cz+(1-t)*h,theta,-%pi,%pi,t,0,1),
parametric_surface(cx+R*cos(theta)*t,
cy+R*sin(theta)*t,
cz-(1-t)*h,theta,-%pi,%pi,t,0,1),
color=black,
geomap([America_coastlines,Africa_coastlines],
[conic_projection,cx,cy,cz,r,a]))$
Back to main page
Back to draw contents