USGS Earthquake Data from Lisp

Kevin Layer has a fun little article in Franz' Tech Corner about processing U.S. Geological Survey earthquake data with ACL.

cl-user(8): (get-quake-info
	     :period :week
	     :larger-than nil
	     :within 1.0
	     :reference-location 
	     (place-to-location
	      "555 12th St, Oakland, CA, 94607, US"))
;; Downloading data...done.
January 09, 2007 22:33:42 GMT: Hollister, California: magnitude=3.1
January 09, 2007 22:16:21 GMT: Los Altos, California: magnitude=2.0
January 09, 2007 21:39:52 GMT: Hollister, California: magnitude=1.3
January 09, 2007 16:56:54 GMT: Aromas, California: magnitude=1.4
January 09, 2007 13:20:02 GMT: Cobb, California: magnitude=1.4
January 09, 2007 11:52:11 GMT: Hayward, California: magnitude=1.2
January 09, 2007 03:01:03 GMT: El Cerrito, California: magnitude=1.9
January 08, 2007 11:30:55 GMT: San Jose, California: magnitude=2.0
January 07, 2007 14:22:56 GMT: Walnut Creek, California: magnitude=1.2
January 07, 2007 09:20:38 GMT: Kenwood, California: magnitude=3.1
January 06, 2007 12:58:01 GMT: Yountville, California: magnitude=1.8
January 06, 2007 07:04:01 GMT: Oakland, California: magnitude=2.4
January 05, 2007 22:55:53 GMT: Walnut Creek, California: magnitude=1.3
January 05, 2007 20:18:51 GMT: Livermore, California: magnitude=2.4
January 05, 2007 17:54:17 GMT: Ladera, California: magnitude=1.1
January 05, 2007 11:09:43 GMT: Cobb, California: magnitude=2.4
January 05, 2007 10:40:26 GMT: Concord, California: magnitude=1.3
January 04, 2007 20:38:03 GMT: Los Altos, California: magnitude=1.7
January 04, 2007 20:21:41 GMT: Aromas, California: magnitude=2.6
January 04, 2007 13:42:36 GMT: Palo Alto, California: magnitude=1.2
January 04, 2007 03:42:10 GMT: Oakland, California: magnitude=1.9
nil

Unfortunately ACL 8.0 Express' heap limit prevented the code from compiling for me. It does some compile-time processing of ZIP code data, and while Kevin's comments mention that the code uses 1/10 as much memory as his first version it's still too much for the free version of ACL.

So for your convenience, this patch modifies the code to do the ZIP code processing at load-time: usgs-acl-express.patch

Once I got it working, I wanted to modify the output a bit to look more like the USGS output, which gives the offset direction and distance of the epicenter from the nearest city:

 MAG    DATE    LOCAL-TIME  LAT     LON    DEPTH    LOCATION
        y/m/d     h:m:s     deg     deg     km

 1.5  2007/01/19 16:54:03 34.187N 117.496W  6.7    9 km ( 6 mi) WSW of Devore, CA
 1.7  2007/01/18 13:37:46 33.823N 117.517W  0.0    7 km ( 4 mi) SE  of Corona, CA
 1.9  2007/01/17 16:08:31 34.320N 118.537W  9.5    6 km ( 4 mi) NNW of Granada Hills, CA
 1.3  2007/01/17 15:42:23 33.978N 117.160W 17.5    9 km ( 5 mi) NE  of Moreno Valley, CA
 1.5  2007/01/17 14:48:04 34.053N 117.260W 17.5    1 km ( 1 mi) WNW of Loma Linda, CA

Kevin's code already figures out the city closest to the quake, so I just pass that to his place-to-location function to get the latitude and longitude of the city from Google (there was one glitch in that the ZIP code data he uses to figure out the nearest city includes some entries that are not actually ZIP codes and don't correspond to any city; I modified the code to not return those). I don't know if Google's city-level geocoding returns the coordinates of the center of mass or what, but testing with Google Maps give some evidence that it returns some basically central point. From that I use distance-between to get distance in miles.

To get the direction of the epicenter from the reference city, I first calculate the bearing from the city to the quake using a formula from Ed Williams' Aviation Formulary:

(defun bearing-to (location1 location2)
  (flet ((radians (deg) (* deg (/ pi 180.0)))
	 (degrees (rad) (* rad (/ 180.0 pi))))
    (let ((lat1 (radians (location-latitude location1)))
	  (lon1 (radians (location-longitude location1)))
	  (lat2 (radians (location-latitude location2)))
	  (lon2 (radians (location-longitude location2))))
      (degrees (mod (atan (* (sin (- lon2 lon1)) (cos lat2))
			  (- (* (cos lat1) (sin lat2))
			     (* (sin lat1) (cos lat2) (cos (- lon2 lon1)))))
		    (* 2 pi))))))

Then I find the nearest cardinal direction to that bearing:

(defparameter *directions*
  '(("N"  . 0)
    ("NE" . 45)
    ("E"  . 90)
    ("SE" . 135)
    ("S"  . 180)
    ("SW" . 225)
    ("W"  . 270)
    ("NW" . 315)))

(defun bearing-to-direction (bearing)
  (let ((best nil)
	(best-diff 0))
    (dolist (dir *directions*)
      (let ((direction (car dir))
	    (heading (cdr dir)))
	(when (or (null best)
		  (< (abs (- bearing heading)) best-diff))
	  (setf best direction)
	  (setf best-diff (abs (- bearing heading))))))
    best))

The final result:

CL-USER> (get-quake-info :period :week
			 :larger-than nil
			 :within 1.0
			 :reference-location (place-to-location
					      "los angeles, ca, us"))
;; Downloading data...done.
January 18, 2007 21:37:46 GMT: 4.6 miles SE of Corona, California: magnitude=1.7
January 18, 2007 00:25:25 GMT: 3.5 miles NW of Boron, California: magnitude=1.5
January 18, 2007 00:08:31 GMT: 3.2 miles N of Porter Ranch, California: magnitude=1.9
January 17, 2007 22:48:04 GMT: 0.3 miles N of Loma Linda, California: magnitude=1.5
January 17, 2007 17:36:43 GMT: 5.9 miles E of Gorman, California: magnitude=1.3
January 17, 2007 00:38:13 GMT: 3.0 miles NW of Boron, California: magnitude=2.2
January 16, 2007 19:57:58 GMT: 6.7 miles N of Santa Paula, California: magnitude=2.1
January 16, 2007 00:47:12 GMT: 3.1 miles NW of Boron, California: magnitude=1.7
January 15, 2007 01:08:47 GMT: 3.2 miles NW of Boron, California: magnitude=1.7
January 14, 2007 19:49:13 GMT: 12.4 miles E of Anaheim, California: magnitude=1.4
January 14, 2007 14:50:31 GMT: 17.1 miles NW of Frazier Park, California: magnitude=1.4
January 14, 2007 09:06:07 GMT: 6.2 miles N of Crestline, California: magnitude=1.7
January 13, 2007 21:01:04 GMT: 6.8 miles SW of Huntington Beach, California: magnitude=2.1
January 13, 2007 20:57:40 GMT: 1.3 miles NW of Newport Beach, California: magnitude=1.9
January 13, 2007 18:43:57 GMT: 6.8 miles N of Piru, California: magnitude=1.7
January 13, 2007 16:31:55 GMT: 4.8 miles E of Yorba Linda, California: magnitude=1.6
January 13, 2007 06:39:57 GMT: 1.0 miles NW of Loma Linda, California: magnitude=2.1
NIL

It looks basically right, but when you compare to even just the first few USGS results there are some differences.

• Devore, CA is about 60 miles east of LA, which probably puts it outside the 1.0 degrees bound, longitude-wise.
• The magnitude 1.7 quake near Corona seems to match up.
• My code has a 1.5 quake near Boron that doesn't seem to match anything on the USGS website.
• The USGS has the magnitude 1.9 quake on 1/18/2007 4 miles north of Granada Hills while my code has it as 3.2 miles north of Porter Ranch; Luckily Granada Hills and Porter Ranch seem to be right next to one another.
• The 1.5 quake near Loma Linda matches, though the USGS calls it 1 mile WNW while I have it as 0.3 miles north.

I am satisfied by this level of consistency.

This patch adds the additional code and includes the ACL 8.0 Express patch: usgs-enhanced.patch.

Posted by jjwiseman at January 19, 2007 06:32 PM