Robin Kraft

How to make a satellite fire map

2017-12-07T00:00:00+00:00

07 December 2017 - Oakland

How to make a satellite fire map

I learned a lot making my fire map, from raster processing to managing a viral webpage. And the nice folks at Mapbox have given me some space on their blog to write about how I made and updated it. In a nutshell, I made heavy use of virtual raster files and automated everything process with a simple bash script.

You can read all the gory details here!

Get your satellite fire maps here!

2017-10-18T00:00:00+00:00

18 October 2017 - Oakland

Wrapping up: Get your satellite fire map here!

The last week has been intense with all the interest in satellite mapping of the Northern California firestorm. And I’ve updated the imagery in the map a few times, prioritizing recency, sometimes at the expense of visual clarity (curse you Smoke!). I’ve tried to strike a decent balance.

Red = veg

Pease remember that red in the false color images shows vegetation, NOT fires. A lot of people have misinterpreted this and gotten needlessly worried. Not surprising given all the stress we’ve been under. but I don’t want to add to it!

Go here to browse the map.

Interactive satellite map of NorCal fires

2017-10-15T00:00:00+00:00

15 October 2017 - Oakland

Interactive satellite map of NorCal fires

Go here to browse the map.

Spatial joins with Hadoop and Cascalog

2012-05-27T00:00:00+00:00

27 May 2012 - San Francisco

Spatial joins with Hadoop and Cascalog

Last time we looked at a simple spatial join of a single latlon with a polygon geometry of the United States. (intersects? poly pt) works like a charm, but it takes a while to complete for the detailed US polygon geometry. Checking millions of points would take weeks at .27 seconds per point.

But! We have Hadoop! This embarassingly parallel problem would be much more tractable if expressed as map-reduce jobs. We’ll let Cascalog handle mapping and reducing, I just have to write a few queries.

The query

Full code is available on Github. Here’s the main query:

(defn join-n-count
  [poly-path points-path intersect-op]
  (let [[iso poly-geom] (import-poly poly-path)
        pts-tap (points-tap points-path)]
    (<- [?count ?iso]
        (pts-tap ?lat ?lon)
        (add-fields iso :> ?iso)
        (intersect-op [poly-geom] ?lat ?lon)
        (c/count ?count))))

join-n-count has two interesting features.

First, and perhaps most importantly, it loads the polygon geometry before the body of the query. import-poly is another Cascalog query that prepares the polygon geometry and extracts the iso code, returning both as Clojure data structures using the ??- operator. Because we execute the query and load the data structure inside the let statement, we can pass it as an argument to the parameterized (see section two) intersect-op operation.

This mainly means that each machine gets its own copy of the data structure, serialized to the job conf. Limits on the size of the job conf means this won’t work with large data structures. An alternative method would use a dreaded cross-join to intersect each point with the polygon, but that runs through a single reducer.

Second, the intersect-op argument to join-n-count expects a parameterized custom operation for intersection. This makes it easy to simply count fires falling inside the US polygon envelope (or bounding box), then do a full intersection. The query itself doesn’t change, but boy the run times sure do! Envelope intersects are much faster, given the simple (rectangular) polygon geometry.

(deffilterop [intersects-op [poly]] [lat lon]
  "Parameterized intersect"
  (let [pt (latlon->point lat lon)]
    (intersects? poly pt)))

(deffilterop [intersects-env-op [poly]] [lat lon]
  "Parameterized intersect only checks envelope intersect"
  (let [pt (latlon->point lat lon)]
    (intersects? (envelope poly) pt)))

The main query has three parts. First, it loads the polygon boundary, then loads all the fires points in the map stage, performing the intersect operation on the reduce side. Each reducer emits a tuple of the number of fires it saw that fell inside the polygon (or its envelope). The final step combines these separate results to give a final total of fires inside the US polygon or its envelope.

Test environment

I’m using Elastic MapReduce from Amazon Web Services, using m2.4xlarge spot instances with 68.4gb of RAM and 8 virtual cores. This works out to 30 mappers and 24 reducers per machine. I’m working with a 1-machine “cluster” and a 10-machine cluster. Set up is simple:

git clone git@github.com:robinkraft/spatialog.git
cd spatialog
lein uberjar

zip -d spatialog-0.1.0-SNAPSHOT-standalone.jar org/apache/xerces/\*
zip -d spatialog-0.1.0-SNAPSHOT-standalone.jar META-INF/services/\*

screen -Lm hadoop jar spatialog-0.1.0-SNAPSHOT-standalone.jar clojure.main

From the REPL:

(use 'spatialog.easyjoin)
(in-ns 'spatialog.easyjoin)
;; sample query
(??- (parameterized-join 
      "s3n://formaresults/test/gadm/usa.csv" 
      "s3n://formaresults/test/firesnoheader/firesnoheader" 
      intersects-env-op))

The lines above deleting files from the jar file are my workaround to an issue with JTS’s dependency on an old version of Apache Xerces. This causes Cascalog jobs to fail for reasons I have not yet uncovered, but all is not lost! If you’re working locally, just delete xercesImpl-2.4.0.jar file from the lib folder in your project. Some day I hope to understand why Xerces causes problems. Stacktrace here.

Results

Of 46 million fires detected from 2000-2012 using the MODIS sensors on Terra and Aqua, 8,331,135 fell inside the bounding envelope of the US border polygon. Of those, 950,055 were detected actually detected within the border of the United States.

Performance

I ran three queries on each cluster on two separate fires datasets, one with just 2700 lines, the other with the full 46 million. Take the exact times with a grain of salt, rather focus on the orders of magnitude. The first run, the envelope intersect on the small dataset and small cluster, established a baseline of 8 min. Running the exact intersect on the polygon for these 2700 fires took 9 min.

The 1-machine cluster C1 took 11 min. for the envelope intersect for the 46 million fires - not meaningfully different from the 2700 fires above. The 10-machine cluster C10 took 13 min. I chalk the longer time for the larger cluster up to the small files problem. Hadoop is designed to power through massive files, not tiny ones, so performance suffers.

The interesting times are for the exact intersect for all fires. C1 never completed because I got tired of waiting. C10 ended up finishing in 5h15, compared to about 13 min. for the envelope intersect. We can infer that C1 likely would have finished in upwards of 50 hours. In total, C10 used 17.5 CPU days for this job, about 1/3 less than I had expected given the performance of my iMac. The total cost of EC2 time was $45, with another $25 for the EMR service.

Discussion

All in all, I’m pleased that this finally worked, although the 17.5 CPU days it took seems excessive. But since I had no other viable options I can’t complain. The complexity of the polygon (40+ thousand vertices) was the main contributor to the 5h15 it took to complete on the larger cluster. An otherwise identical query intersecting polygon envelopes took only 13 minutes, or 25x less time. Were it not for the small size of the files that add to Hadoop overhead, I would expect the envelope query to approach the 100x less time it takes to do the envelope intersection at the REPL. For larger files, Hadoop’s IO model really gets fast, but for small files it never gets a chance to ramp up.

This Hadoop job only does a small part of what I’d like to do: count up fires for all countries. But the run time for what I have done makes it seem infeasible to do the join for all countries without moving to highly customized workflows and custom code. Then again, I didn’t try building a spatial index, something that could make this more feasible.

As it stands, the success of this exploration of spatial joins rests on the artificial simplicity of the question. But I am convinced that between Hadoop distributed caches, simplified polygons, and some kind of indexing, it will be possible to count up fires by country without spending too much or waiting months for the job to finish.

Fires data and spatial joins with Clojure and JTS

2012-05-26T00:00:00+00:00

26 May 2012 - San Francisco

Fires data and spatial joins with Clojure and JTS

I’ve got a soft spot for spatial joins. This is partly because they’re so darn useful, but also because a spatial join with raster data (e.g. the ArcGIS Sample tool) was one of the first GIS operation I coded up by hand, using Python and Numpy. Ah, those were the days …

But today, we’re going to take the easy route, and use libraries to do a spatial join with Clojure. Let’s find out how many fires were detected in the US since late 2000.

Data

I’ve got some 46 million points representing fires detected by the MODIS sensors flying on NASA’s Terra and Aqua satellites over the last 11 years. You can get more info and all the data from the UMD FIRMS web site. We’ll be processing a 2.6gb textfile of fire points, and using polygon data from the Global Administrative Areas database (GADM).

I happen to have the GADM data sitting on a CartoDB instance, so let’s just grab some well-known-text in CSV format using the CartoDB SQL API:

http://your-account.cartodb.com/api/v2/sql/?q=SELECT iso,ST_asEWKT(the_geom)
AS geom FROM gadm0 WHERE iso='USA'&format=csv

That was easy! You can grab a zipped copy here if you want to follow along at home. I’ve stripped out the iso,geom header line for convenience.

The file looks like this:

USA,"SRID=4326;MULTIPOLYGON(((179.585615158081    
52.0173091888428,179.652040481567 52.0247783660889 ...)))"

A line from the fires data represents one 1km^2 pixel, and looks like this:

-18.735,144.626,324.8,1.1,1.0,05/13/2012,0030,T,78,5.0,297.2,24.9

A full description of the attributes is available on the FIRMS site. We’re only interested in the first two fields, latitude and longitude - the rest includes brightness, confidence, date, time, etc.

Setup

Add [cljts "0.1.0"] to project.clj to get cljts and JTS. Check out this post for some simple use cases. Cljts pulls in the raw JTS library if you want to use it, but also gives you some convenient helper functions. Thank you sunng! Be sure to check out the jts and cljts docs.

Basic processing

Cljts has a nice read-wkt-str function that we’ll be using here. Let’s set up the environment first:

(use '[cljts.io :only (read-wkt-str)])
(use '[cljts.relation :only (intersects?)])

Then we load the polygon, remembering to strip out some of the metadata provided by CartoDB:

(defn load-poly
  [path]
  (let [poly-line (slurp path)
        [iso, poly-field] (clojure.string/split poly-line #"," 2)
        wkt (subs poly-field 11)]
    (read-wkt-str wkt)))

(def poly (load-poly "/Users/robin/data/usa_poly.csv"))
(envelope poly)
=> #<Polygon POLYGON ((-179.15055847168 18.9098606109619,
    -179.15055847168 71.3906841278076, 179.773408889771
    71.3906841278076, 179.773408889771 18.9098606109619,
    -179.15055847168 18.9098606109619))>
(count (coordinates poly))
=> 40411

Now we’ve got a JTS multipolygon object. Let’s create a fire point based on a coordinate:

(def pt (point (c 144.626 -18.735))) ;; note the xy order
=> #<Point POINT (144.626 -18.735)>

Now we can check intersection:

(intersects? poly pt)
=> false

That point happens to be in northeast Australia. So let’s try another random latlon, from Yellowstone National Park:

(intersects? poly (point (c -110.560913 44.801327)))
=> true

Now we need a spatial join function to handle the polygon attributes.

(defn spat-join 
  [poly-map point]
  (if (intersects? (:geom poly-map) (point)
      {:geom point :attributes (:attributes poly-map)})))

(def usa-poly-map {:geom poly :attributes {:iso "USA"}})
(def yellowstone-pt (point (c -110.560913 44.801327)))
=> #<Point POINT (-110.560913 44.801327)>

(spat-join usa-poly-map yellowstone-pt)
=> {:geom #<Point POINT (-110.560913 44.801327)> :attributes {:iso "USA"}}

(spat-join usa-poly-map (point (c 1 1)))
=> nil

It works! Now we just need to loop through 46 million fires points, storing the count by iso code, and then we would have a spatial join.

JTS performance

So, looping through fire points line by line, intersecting each one with the polygon geometry in memory, would technically work. But it would take a while. On my 2010 iMac, intersect? returns false in as few as .0025 milliseconds. I suspect that intersect? first checks whether a point is within the bounding envelope of the polygon of interest - that operation screams. If the point doesn’t fall within the bounding envelope, intersects? immediately returns false. If it does, I ended up waiting 268 milliseconds on average. Granted, the border of the US and its territories is a complex polygon, with 40+ thousand vertices, but I would love to see it faster.

I don’t know enough (read anything) about computational geometry to know how hard it would be to do better than this. JTS uses uses a simple polygon intersection algorithm to check for point location. Tim Robertson, who helped inspire this exploration, has a nice implementation of a slab decomposition algorithm that is incredibly fast. But as I said at the beginning, for now I want to stick with libraries, hoping to retain some semblance of a general solution.

Conclusions

This exercise is going to be a problem. I happen to know (foreshadowing the next blog post) that of 46 million fires, 8.33 million fall within the USA’s bounding envelope, which apparently stretches across most of the northern hemisphere, from Alaska to perhaps Guam. From the GADM site:

So the question is, ignoring super-fast non-intersecting operations, how long should the exact intersection take?

(/ (/ (/ (* 8300000 .268) 60) 60) 24.)
=> 25.8 days

That’s crazy! For completeness, note that ArcGIS chokes with a memory error after about 5 hours of work on a spatial join.

If only this ran in parallel … (stay tuned)

I’m satisfied with the ease of using JTS in Clojure thanks to cljts, but this level of performance just won’t do for the task at hand. Let’s not forget that I’ve simplified the question. My real question - how many fires were detected in each country? - would take ages to complete given the size of the bounding envelopes of the US and Russia alone. There has to be a better way.

TL;DR

Cljts is a great way to use JTS Topology Suite in Clojure, and has some convenient helper functions for IO and spatial relationships. But a naive loop through 46 million fires points to count up fires intersecting the US would take nearly 26 days using JTS’ intersection algorithm. A custom solution could be much faster, at the loss of generality. Look for a future post on trying this with Hadoop and Cascalog.

JTS Topology Suite and Clojure

2012-04-15T00:00:00+00:00

15 April 2012 - San Francisco

JTS Topology Suite and Clojure

This is the first in a series of posts about using the JTS Topology Suite (JTS) in Clojure, particularly in conjunction with Cascalog and Hadoop. To start I’ll briefly look at the setup needed for using JTS with Clojure. Further along I’ll explore replicating GIS operations using Cascalog/Hadoop.

JTS Topology Suite is “an API for modelling and manipulating 2-dimensional linear geometry”. That really just means you can use it to do typical GIS operations, like intersecting polygons or determining whether a line touches a polygon, among many other things.

Written in Java, JTS should work directly with Clojure, but there are already at least three partial wrappers for JTS that provide greater convenience.

jts-clj provides a Clojure-ish API for creating JTS geometry instances.
cljts’s geometry API feels like less idiomatic Clojure at first glance, but it does wrap some of JTS’ spatial analysis and spatial relationship functions.
geoscript-clj aims more broadly at working with geospatial data in Clojure, using GeoTools and JTS.

I haven’t worked with these enough to pick a favorite, but I started with cljts because of the spatial functions and because it seems simpler than geoscript-clj. I’ll probably revisit this at some point.

To get started with cljts, add [cljts "0.1.0"] to your project.clj file. Then you can create a simple coordinate geometry:

user> (use 'cljts.geom)
user> (c 10 20)
#<Coordinate (10.0, 20.0, NaN)>

And there you have it. Create a polygon by stringing together coordinates in a vector, being careful to close the polygon by making the first and last coordinates the same. Note that polygon expects two rings - the external border and a hole - but you can create a polygon without a hole by replacing the hole linear-ring with nil:

user> (polygon (linear-ring [(c 20 40) (c 20 46) (c 34 56) (c 20 40)]) nil)
#<Polygon POLYGON ((20 40, 20 46, 34 56, 20 40))>

cljts provides an interface to JTS’ DE9-IM spatial relationship methods. You can check intersection using intersects?.

user> (def poly (polygon (linear-ring [(c 20 40) (c 20 46) (c 34 56) (c 20 40)]) nil))
user> (def pt (point (c 21 45)))
user> (intersects? pt poly)
true

Other functions for determining spatial relationships are described in the cljts.relation documentation.