;;;; Mapping the continent divide. ;;;; Written by Brian Hayes, 2008-2009, based on ;;;; earlier work done in 2000. ;;;; Author contact: brian@bit-player.org ;;;; I'm making this freely avalable for any use. Acknowledgement ;;;; of authorship would be appreciated. ;;;; For background on where this idea came from, and references ;;;; to related work by others, see "Dividing the Continent," ;;;; by Brian Hayes, American Scientist Vol. 88, No. 6, pp. 481-485, ;;;; http://www.americanscientist.org/issues/pub/dividing-the-continent ;;;; See also http://bit-player.org/2009/long-division ;;;; The digital elevation maps on which these programs work ;;;; are available at http://www.ngdc.noaa.gov/mgg/topo/globe.html ;;;; The code was developed in Clozure Common Lisp ;;;; Version 1.2-r11300M (DarwinX8664) [see http://trac.clozure.com/openmcl] ;;;; I see no reason it shouldn't run in any other CL, but please ;;;; send feedback if I'm wrong about that. ;;; UTILITIES ;;; Some weird stuff from my init.lisp file that ;;; I can't live without. ;; Syntax for the Scheme-style named let. ;; This is supposed to transform ;; ;; (defun factn (n) ;; (nlet loop ((n n) (prod 1)) ;; (if (<= n 1) ;; prod ;; (loop (1- n) (* n prod))))) ;; ;; into ;; ;; (defun factn (n) ;; (labels ((loop (n prod) ;; (if (<= n 2) ;; prod ;; (loop (1- n) (* n prod))))) ;; (prod n 1))) (defmacro nlet (name argbinds &rest body-exprs) `(labels ((,name ,(mapcar #'car argbinds) ,@body-exprs)) (,name ,@(mapcar #'cadr argbinds)))) ;; The old-fashioned while loop (defmacro while (test &rest body) `(do () ((not ,test)) ,@body)) ;; And the Pascal-style writeln -- why not? (defun writeln (destination &rest args) (nlet loop ((args args)) (cond ((null args) (terpri destination)) ((null (cdr args)) (princ (car args) destination) (terpri destination)) (t (princ (car args) destination) (princ " " destination) (loop (cdr args))))) (values)) ;; Works just like the primitive 'list', except that it ;; ignores any arguments that have the value nil. Thus ;; (list-if 1 2 3 nil 5 nil) returns (1 2 3 5) (defun list-if (&rest args) (remove-if #'not args)) ;; Encoding bytes as hexhex values for ASCII Postscript files. ;; Given a byte (that is, an integer in [0, 255]), return a ;; two-character string with its hexadecimal representation. (defun byte->hexhex (byte) (let ((vec '#("0" "1" "2" "3" "4" "5" "6" "7" "8" "9" "A" "B" "C" "D" "E" "F"))) (multiple-value-bind (hi lo) (floor byte 16) (concatenate 'string (elt vec hi) (elt vec lo))))) ;;; Here is the single big procedure that does the analysis ;;; of the map to identify the continental divide. ;;; Some preliminary explanation. We start with a big array ;;; of pixel elevations, which are nonnegative integers. (See notes ;;; below on preliminary data preparation.) For each pixel ;;; we need a tag variable that indicates the pixel's current state, ;;; and so we set up a second array of identical size and ;;; shape to hold the tags. The only other major data structure ;;; is the pixel-queue, which we can view as an abstract data ;;; type with enqueue and dequeue procedures and a queue-empty-p ;;; predicate. ;;; The basic routine goes like this. Initially, all pixels ;;; have the tag "dry", then we choose some seed or starter ;;; pixels that are assigned tags corresponding to the oceans ;;; or other basins. Each seed is associated with a color. In ;;; an initialization step we identify the neighbors of each ;;; seed pixel and put them in the pixel-queue. Then, until the ;;; queue is empty, we repeat the following steps: ;;; ;;; 1. Withdraw a pixel p from the queue. ;;; 2. If p is tagged anything other than "dry", discard it. ;;; (Because we must have seen it before.) ;;; 3. List all the colored neighbors of p. (The list ;;; should not be empty, because pixels get onto ;;; the queue by virtue of being neighbors of ;;; colored pixels.) ;;; 4. If all the colored neighbors are the same color, ;;; assign that color to p also. Then find p's ;;; neighbors, and add them to the queue. ;;; 5. Evidently p has neighbors of different colors, ;;; so it is on the divide. Tag it red. ;;; 6. Repeat. ;;; Guide to arguments: ;;; ;;; mapvec : vector of elevations as nonnegative integers ;;; cols : width of the map in pixels ;;; rows : height of the map in pixels ;;; (assert (= (* cols rows) (length mapvec))) ;;; print-sched : list of elevations where the program is to ;;; stop and print the current state of the map ;;; The initial and final states need not be ;;; listed, they'll be printed anyway. ;;; starters : list of pixels from which the filling of ;;; the oceans will begin. Each element of this ;;; list can be either an integer in the range ;;; 0 <= p < (length mapvec) or a list of such ;;; integers. In the latter case, all the pixels ;;; in the sublist are assigned the same color. ;;; For example, '(0 100 (200 300)) creates a ;;; map with four starter pixels but only three ;;; distinct colors -- pixels 200 and 300 share ;;; the same color. (The sublist facility is ;;; included to solve the Gulf of California ;;; problem: Because the map doesn't extend all ;;; the way to the bottom of Baja California, ;;; the program sees the Gulf of California as ;;; a body of water separate from the Pacific ;;; Ocean.) Internally, colors are integers ;;; assigned sequentially from 1. How these tags ;;; are intepreted visually is determined by the ;;; clut function. ;;; clut-fn : functional argument giving a color lookup table. ;;; This is a function of two arguments, a pixel ;;; elevation and a pixel tag. It should return ;;; three values -- hex-encoded levels of red, ;;; green and blue. Sample functions below. Note ;;; that you need to pass this as (function clut-fn-name) ;;; or #'clut-fn-name. ;;; filename-stem : map files are named 'filename-stem-nnn.eps' where ;;; nnn is an integer giving the water level in the ;;; map. Each nnn corresponds to one of the elevation ;;; levels in print-sched (or the next higher level ;;; if there is no pixel at the specified level). (defun global-warming (mapvec ; the one-dimensional array of elev data cols rows ; x and y dimensions of the map print-sched ; elevations at which we print a snapshot of the map starters ; seed pixels that get initial colors clut-fn ; how to color the map filename-stem) ; how to name map.eps files (let* ((pixel-count (length mapvec)) (max-level (reduce #'max mapvec)) ; highest pixel in the landscape (pending-pixels (make-array (1+ max-level) :element-type 'list :initial-element nil)) ; the pixel queue (least-elevation 0) ; lowest occupied level ; in the queue (tagvec (make-array pixel-count :element-type 'fixnum :initial-element 0)) ; a tag for each pixel (red-tag -1) ; for pixels on a divide (dry-tag 0) ; for pixels not yet underwater (start-of-last-row (* cols (1- rows))) ; used in south, below (modlastcol (1- cols))) ; used in east, below (macrolet ; convenient syntax for accessing ((get-tag (p) ; the arrays of elevations and tags `(aref tagvec ,p)) ; we don't need a 'set-elev' macro (set-tag (p tag-value) ; because we're not going to alter `(setf (aref tagvec ,p) ,tag-value)) ; the landscape (get-elev (p) `(aref mapvec ,p))) (labels ; local procedures ((north (p) (if (< p cols) nil (- p cols))) ; the apparatus needed to (south (p) (if (>= p start-of-last-row) nil (+ p cols))) ; treat our one-dimensional (west (p) (if (zerop (mod p cols)) nil (1- p))) ; vector of pixels as if it (east (p) (if (= modlastcol (mod p cols)) nil (1+ p))) ; were a two-dimensional array (list-neighbors (p) (list-if (north p) (south p) (west p) (east p))) ; 'list-if' drops nil neighbors ;; the next three procedures define the interface to the pixel-queue ;; it's not really a queue but an array of stacks, with an auxiliary ;; variable 'least-elevation' that keeps track of the lowest nonempty ;; element of the array. The guarantee is that when we pull a pixel ;; from the queue (using dq-pixel) it will be one of those at the lowest ;; nonvacant elevation. We don't care which of the pixels at that elevation ;; is returned, just that the sequence of dequeued pixels is always in ;; elevation order. (nq-pixel (p) (let ((elev (get-elev p))) ; look up p's elevation (push p (aref pending-pixels elev)) ; push it onto the corresponding stack (when (< elev least-elevation) ; if this is the lowest, update (setf least-elevation elev)))) ; the pointer to least-elevation (dq-pixel () (if (pq-empty-p) ; empty queue returns nil nil (prog1 (pop (aref pending-pixels least-elevation)) ; return top item on lowest nonempty stack (when (endp (aref pending-pixels least-elevation)) ; if that stack is now empty, we need (nlet loop ((j (1+ least-elevation))) ; to readjust 'least-elevation' by (cond ((> j max-level) ; searching for next nonempty stack. (setf least-elevation 0)) ; if none found, set 'least-elevation' ((not (null (aref pending-pixels j))) ; to zero (setf least-elevation j)) (t (loop (1+ j))))))))) (pq-empty-p () ; queue is empty if 'least-elevation' (and (zerop least-elevation) ; points to 0, and the stack at (null (aref pending-pixels least-elevation)))) ; elevation 0 is empty ;; Here's where we do all the work. This is called on each ;; pixel drawn from the pixel-queue, and assigns it a color ;; based on the colors of the adjacent pixels. (flood-a-pixel (p) (when (= (get-tag p) dry-tag) ; if not, we've already dealt with this pixel so skip it (let* ((neighbors (list-neighbors p)) ; get all the neighbors and (neighbor-tags ; their tags, then keep only (mapcar #'(lambda (x) (aref tagvec x)) neighbors)) ; those >= 1 -- the ocean colors (wet-neighbor-tags ; We don't care about neighbors (remove-if #'(lambda (n) (< n 1)) neighbor-tags))) ; 0 (dry) or -1 (on the divide) (when (null wet-neighbor-tags) ; Shouldn't happen. How did a pixel (cerror "Uh oh. All neighbors of ~D are dry.~%" p)) ; without wet neighbors get enqueued? (if (not (apply #'= wet-neighbor-tags)) ; neighbors or differet colors? (set-tag p red-tag) ; Yes: we're on the divide (progn ; No: we're in the sea (set-tag p (car wet-neighbor-tags)) ; match neighbors' color (nlet loop ((neighbors neighbors) (neighbor-tags neighbor-tags)) ; loop to add all (unless (endp neighbors) ; the dry neighbors (when (= (car neighbor-tags) dry-tag) ; of the newly flooded (nq-pixel (car neighbors))) ; pixel to the queue (loop (cdr neighbors) (cdr neighbor-tags))))))))) ;; This is called for each of the starter pixels (install-a-starter (p tag-value) (format t "Installing starter ~D with tag ~D~%" p tag-value) (set-tag p tag-value) ; the initial drop of color (let ((neighbors (list-neighbors p))) ; add neighbors of colored (mapcar #'nq-pixel neighbors)))) ; pixel to the pixel-queue ;; end of local macros and functions ;; The initialization loop that assigns colors (integers 1, 2, 3...) to ;; the starter pixels and calls 'install-a-starter' on them. (nlet init-loop ((basin-tag 1) (starters starters)) (unless (endp starters) (if (listp (car starters)) ; handle sublists of seed pixels (mapcar #'(lambda (x) (install-a-starter x basin-tag)) (car starters)) (install-a-starter (car starters) basin-tag)) (init-loop (1+ basin-tag) (cdr starters)))) (format t "Flooding...~%") ;; print map in the initial state (format t "~%Printing initial map at elevation ~D.~%" least-elevation) (let ((pathname (concatenate 'string filename-stem "-initial.eps"))) (print-tagged-map mapvec tagvec cols rows clut-fn pathname)) ;; The main loop. Pull a pixel from the queue -- it will be one ;; of those at the lowest unflooded elevation -- and hand it off ;; to flood-a-pixel. Then, if we've reached one of the benchmark ;; elevations in print-sched, print the current state of the map. (while (not (pq-empty-p)) (flood-a-pixel (dq-pixel)) ; handle one pixel (when (and (not (endp print-sched)) ; is it time to print map? (>= least-elevation (car print-sched))) ; (format t "~%Printing map at elevation ~D.~%" least-elevation) ; okay, print it (let ((pathname (concatenate 'string filename-stem ; "-" ; (format nil "~D" least-elevation) ; ".eps"))) ; (print-tagged-map mapvec tagvec cols rows clut-fn pathname)) ; (setf print-sched (cdr print-sched)))) ; drop elev from list ;; print map in final state (format t "~%Printing final map at elevation ~D.~%" max-level) (let ((pathname (concatenate 'string filename-stem "-final.eps"))) (print-tagged-map mapvec tagvec cols rows clut-fn pathname)))))) ;;; END OF 'global-warming' ;;; Create an Encapsulated Postscript image of the map. This is called ;;; at various intervals by 'gloabl-warming' The output is in ASCII ;;; Postscript, which should be viewable in Ghostscript, Photoshop, ;;; etc., but you should also be able to inspect it with a text editor. (defun print-tagged-map (mapvec tagvec cols rows clut-fn pathname) (let* ((max-elev (reduce #'max mapvec)) (min-elev (reduce #'min mapvec)) (elev-range (- max-elev min-elev)) (pixel-count (* cols rows))) (labels ((normalize (i) ; convert elevation to range [0.0 .. 1.0] (/ (- i min-elev) elev-range))) (with-open-file (outfile pathname :direction :output :if-does-not-exist :create :if-exists :overwrite) (writeln outfile "%!PS-Adobe-3.0 EPSF-3.0") (writeln outfile "%%Creator: contintental-divide.lisp") (writeln outfile "%%BoundingBox:" 0 0 cols rows) (writeln outfile "%%HiResBoundingBox:" 0 0 cols rows) (writeln outfile "%%EndComments") (writeln outfile "%ImageData:" cols rows 8 3 0 1 2 "\"imageproc\"") (writeln outfile "/pixelstring" (* 3 cols) "string def") (writeln outfile "0 0 translate") (writeln outfile cols rows "scale") (writeln outfile "/imageproc") (writeln outfile " {" cols rows "8") (writeln outfile " [" cols "0 0" (- rows) "0" rows "]") (writeln outfile " { currentfile pixelstring readhexstring pop }") (writeln outfile " false 3 colorimage } def") (writeln outfile "imageproc") (nlet loop ((j 0) (linewidth 0)) (cond ((>= j pixel-count)) ((>= linewidth 72) (terpri outfile) (loop j 0)) (t (multiple-value-bind (r g b) (funcall clut-fn (normalize (aref mapvec j)) (aref tagvec j)) (princ r outfile) (princ g outfile) (princ b outfile)) (loop (+ j 1) (+ linewidth 6))))) (terpri outfile) (writeln outfile "showpage"))))) ;;; Clut functions. Give me an elevation and the corresponding tag, I'll ;;; give you back a color. The output conssts of three strings, each ;;; of which is two characters encoding the hex value of an RGB color ;;; channel. ;;; ;;; The clut functions given here are the ones used for the published ;;; versions of the maps. They set the luminosity of a pixel to the ;;; square root of the pixel's normalized elevation. (A linear mapping ;;; looks too dark.) Then tints are added based on the tag values. (defun clut2-rg-sqrt (elev tag) (let ((lum (byte->hexhex (round (* (sqrt elev) 255)))) (lum+ (byte->hexhex (round (+ 127 (* (sqrt elev) 128)))))) (case tag (0 (values lum lum lum)) (2 (values lum lum lum+)) (1 (values lum lum+ lum)) (-1 (values "FF" "00" "00")) (-2 (values "0F" "0F" "0F")) (t (values "FF" "FF" "00"))))) (defun clut3-rgy-sqrt (elev tag) (let ((lum (byte->hexhex (round (* (sqrt elev) 255)))) (lum+ (byte->hexhex (round (+ 127 (* (sqrt elev) 128)))))) (case tag (0 (values lum lum lum)) (3 (values lum+ lum+ lum)) (2 (values lum lum lum+)) (1 (values lum lum+ lum)) (-1 (values "FF" "00" "00")) (-2 (values "0F" "0F" "0F")) (t (values "FF" "FF" "00"))))) (defun clut4-rgyc-sqrt (elev tag) (let ((lum (byte->hexhex (round (* (sqrt elev) 255)))) (lum+ (byte->hexhex (round (+ 127 (* (sqrt elev) 128)))))) (case tag (0 (values lum lum lum)) (4 (values lum lum+ lum+)) (3 (values lum+ lum+ lum)) (2 (values lum lum lum+)) (1 (values lum lum+ lum)) (-1 (values "FF" "00" "00")) (-2 (values "0F" "0F" "0F")) (t (values "FF" "FF" "00"))))) ;;; TEST DATA #| Here's a toy elevation map, with a zigzag divide down the middle, a peak on one side and a basin on the other: 0 0 0 0 1 1 1 2 1 1 0 0 0 0 0 0 0 0 1 1 1 2 1 1 1 1 0 1 2 1 0 0 1 0 1 3 2 1 1 1 0 0 2 0 1 0 2 4 2 1 2 3 2 1 1 1 0 2 2 2 0 0 1 1 2 2 2 3 2 2 1 0 0 0 0 0 0 0 0 1 1 2 2 3 1 1 0 0 0 0 0 0 0 1 1 1 2 2 2 2 1 1 0 0 0 We encode the data in an elevation array as follows: (setf test-map (make-array (* 15 7) :element-type 'fixnum :initial-contents '(0 0 0 0 1 1 1 2 1 1 0 0 0 0 0 0 0 0 1 1 1 2 1 1 1 1 0 1 2 1 0 0 1 0 1 3 2 1 1 1 0 0 2 0 1 0 2 4 2 1 2 3 2 1 1 1 0 2 2 2 0 0 1 1 2 2 2 3 2 2 1 0 0 0 0 0 0 0 0 1 1 2 2 3 1 1 0 0 0 0 0 0 0 1 1 1 2 2 2 2 1 1 0 0 0))) Now we call 'global-warming' on this data: (global-warming test-map 15 7 '(1 2) '(0 104) #'clut2-rg-sqrt "test-map") The output is four files, named test-map-initial.eps, test-map-1.eps, test-map-2.eps and test-map-final.eps. Here are the contents of the final file: %!PS-Adobe-3.0 EPSF-3.0 %%Creator: contintental-divide.lisp %%BoundingBox: 0 0 15 7 %%HiResBoundingBox: 0 0 15 7 %%EndComments %ImageData: 15 7 8 3 0 1 2 "imageproc" /pixelstring 45 string def 0 0 translate 15 7 scale /imageproc { 15 7 8 [ 15 0 0 -7 0 7 ] { currentfile pixelstring readhexstring pop } false 3 colorimage } def imageproc 007F00007F00007F00007F0080BF8080BF8080BF80FF00008080BF8080BF00007F00007F 00007F00007F00007F007F00007F00007F0080BF8080BF8080BF80FF00008080BF8080BF 8080BF8080BF00007F8080BFB4B4DA8080BF007F00007F0080BF80007F0080BF80FF0000 B4B4DA8080BF8080BF8080BF00007F00007FB4B4DA00007F8080BF007F00B4DAB4FFFFFF B4DAB480BF80B4DAB4FF0000B4B4DA8080BF8080BF8080BF00007FB4B4DAB4B4DAB4B4DA 007F00007F0080BF8080BF80B4DAB4B4DAB4B4DAB4FF0000B4B4DAB4B4DA8080BF00007F 00007F00007F00007F007F00007F00007F00007F0080BF8080BF80FF0000B4B4DADDDDEE 8080BF8080BF00007F00007F00007F00007F007F00007F00007F0080BF8080BF8080BF80 FF0000B4B4DAB4B4DAB4B4DA8080BF8080BF00007F00007F00007F showpage |# #| Further notes on acquiring and preparing data for mapping. The NGDC web site offers elevation data in several formats; I have worked only with the plain-ASCII versions. What you get when you download the data is a long list of integers, guaranteed to be in the range [-32768 .. 32767]. (Since they represent terrain elevations in meters, the real range is considerably smaller.) The value -37768 is reserved for missing data; no pixels were flagged as missing in any of the data sets I downloaded. All of the maps in the American Scientist article and in the bit-player posting cover the same geographic range: 125W to 65W longitude and 50N to 25N longitude. I decided (perhaps unwisely -- see below) to replace all negative elevations with 0. Once you have the data read into a Lisp sequence (either a list or a vector), you can zero-out the negative values as follows: (setf mapvec (substitute-if 0 #'minusp mapvec)) In this and other operations, you might want to be careful to invoke commands in such a way that the REPL doesn't try to print the whole vector, which could take a while for vectors of several million elements. |# ;;; Afterthoughts. How I might have done it differently if I had it ;;; to do over again. ;;; ;;; 1. Zeroing out the negative values in the elevation array ;;; was a mistake. I did it because I didn't really care about ;;; the bathymetry and I didn't want those underwater features ;;; to distract from the coastline in the printed maps. But ;;; the bathymetry could have been suppressed at output time, ;;; leavning more options open. ;;; ;;; 2. Keeping the one-dimensional form of the data was a mistake. ;;; Several parts of the computation would have been clearer if ;;; I'd created a two-dimensional array. And 2D array access ;;; is probably more efficient than calculating the indices ;;; by hand. ;;; ;;; 3. To avoid the messy problems of edge pixels (making sure that ;;; neighbors exist, etc.) it might make sense to embed the ;;; whole array in a slightly larger array. Then pixels beyond ;;; the edge would exist, so that trying to access them would not ;;; raise an array-bounds error. But they would be flagged as ;;; invalid. I tried this at one point and then reverted. I'd ;;; be inclined to try again. ;;; ;;; 4. A binary format for the output images would save lots of ;;; space and time. It would also offer the option of 16-bit-per- ;;; channel encoding. ;;; ;;; 5. Everything in this program is based on the 4-neighborhood ;;; (N, S, E, W). It might be fun to see what happens if we work ;;; instead in the 8-neighborhood (N, S, E, W, NE, NW, SE, SW).