Reproducing Spatial Joins using Hadoop and EC2

Finally I got some time to continue researching the use of Hadoop for the processing of the GBIF dataset with respect to the protected areas of the world.

The problem
There is a growing number of point data, and an increasing need to perform cross referencing of the GBIF data with external sources, in this case the protected areas of the world.

The last 3 years of working on reasonably large databases (>100G) has told me that this join was not going to be an easy one, despite the obvious partitioning strategy that can be invoked here (spatial partitioning).
Therefore Javi and I set about this task by

1. Using PostGIS and the traditional approach of joining 2 tables of 150M point records with 120,000 polygon records - finding will be posted shortly
   
2. Using Hadoop and the Amazon EC2 to process the join

Using Hadoop
This is ongoing but I thought I'd post the findings so far.

Strategy 1: 

Partition by cell
Start by splitting the 150M records from one file into a file per 1 degree cell.
Loop over each polygon:

• Determine 1 degree cells covered polygon
  
• Use the input files as a merge input to hadoop and pass in the Polygon for the Map to do the contains()
  

Result  

Hadoop will not handle so many file splits no matter what configuration I try so far. Mailing list responses suggest this is not the ideal use of Hadoop and it prefers single large files, rather than many smaller ones.

Amendment

Using 10 degree cells, it does work but takes a LONG time... so long I killed it. This was expected since it is really the equivalent of only a DB join with an index on an INTEGER teg_deg column.

Strategy 2: 

1-Degree in-memory index of polygons
So this approach goes through the join the opposite way. The polygons are loaded into memory and a

HashMap<Integer, Collection<Polygon>>

for each Map operation to use.  

The logic is effectively give me polygons intersecting the 1 degree square for this point, and then test them properly. This worked ok, but requires a LOT of memory for the polygons to be held.

Result

About 6 hours to process 60,000 polygons.

Strategy 3: 

RTree index of polygons
Very similar to strategy 2, this approach indexes all the polygons by RTree. Thus for each point, the index provides only candidate polygons whose bounding box contains the point, and the number of candidate polygons to test is far smaller than that of strategy 2.

Result

45 minutes for a 130M x 30,000 cross reference but requires the larger EC2 instance sizes. The RTree index means it will need processed in batches.
Hadoop requires tuning, just like a database does for large joins. Here is the final tuning for strategy 3:

• mapred.child.java.opts=-Xmx1G // for the R-Tree index
• mapred.output.compress=false // for easy reading of output
• mapred.tasktracker.map.tasks.maximum=4 // we have 7G per node, so 4G for 4 Map tasks is ok
• mapred.job.reuse.jvm.num.tasks=-1 // Reuse the Map as it has the index built (takes 90 seconds for startup)
• dfs.block.size=134217728 // use larger blocks - in retropsect this was a probable mistake and will be removed for next test

More to follow...

Thanks to the Hadoop-core Users mailing list for their insights, and Andrea of OpenGEO for pointing me at the JTS RTree implementation

Convex Hull over GBIF "points"

I found yesterday a post where Convex Hull was explained and the source code made available for Javascript and PHP. "The convex hull may be easily visualized by imagining an elastic band stretched open to encompass the given points; when released, it will assume the shape of the required convex hull." (from wikipedia). I talked about it months ago with Tim and I thought that would be fun to adapt it to AS3 and use it in BiodiversityAtlas. Months ago we talked that one way to create "polygons" from gbif straight points was to do a Convex Hull over them.

But when finished I realized it might dont make much sense... the convex hull over GBIF points always display enormous polygons over land and water that dont look very nice. Look at the next picture.

This is Puma concolor convex hull over GBIF data. Well this is exactly what a convex Hull is so I am not sure what I was expecting.

There is other types of convex hull like the Orthogonal convex hull that might make more sense in some scenarios, but in general this does not look good to me.

By the way, the convex hull algorithm i am using, quickhull, is amazingly fast.

If you want to try for yourself I have uploaded the test app for you to try. If anybody is interested in the source code let me know and I will post it.

Launch application.

Mapreduce for species in a Protected Area

The World Commission on Protected Areas (WCPA) publish shapefiles for protected areas, so I thought I'd run these against the GBIF index and get some species lists.
Since I was running on the full 135 million record GBIF index, and the IUCN Categories I-VI National shapefiles containing 60,753 polygons, I went back to Hadoop running in single node mode (i.e. the simplest of simple) as I am not sure the Mapreduce 'Lite' I wrote would cut it.

Input file: 135 Million records of 12 DwC fields (13G)
Stack: Hadoop for Mapreduce, Geotools for Polygon

I went for the brut force approach; the Map held a List<Polygon> and I looped each time over the polygons for each GBIF point record (60,000 x 135M = 8,100,000,000,000 comparisons ;o)

Of course this did not perform... but I wanted to post some benchmarks:

Reference: A Map that does nothing, and the GBIF index input: 500secs
Per Polygon: 1500secs to produce the species by polygon
(Note: This is running Mapreduce in a non clustered, single server environment)

I tried this using the polygon bounding boxes only, to remove Geotools from the equation in each Map operation and the results were the same.

It is clear that to do this kind of processing, the data needs to be sliced up to reduce the number of combinations.
I am pondering using the same tiling algorithm that the map guys seem to have settled upon (e.g. GE superoverlay style). By doing this, a preselect based on the intersect of the protected area bounding box with the tiled data would result in massively reduced processing. I am thinking 7 zoom levels resulting in 32768 distinct tiles so therefore working only on 2.8x2.8 degree cells. Currently GBIF models to 1x1 degree and 0.1x0.1 degree cell but this does not easily port to Mapreduce for partitioning purposes of the input data.
(As a side effect, the tiles could be processed in parallel for specific mapping views ready for immediate serving through Geoserver...)