Thursday, December 22, 2016

Visualize your ideas using Rasem

A major part of a researchers' work is to write papers and articles that describe their work and make posters and presentations to better communicate their ideas. We all believe that "A picture is worth a thousand words" and we are always looking for better ways to visualize our ideas. In this blog article, I present Rasem, a library that I built as I started my PhD and used it in many of my papers and presentation to build nice visualizations like the ones shown below.

The figures shown above are generated from hundreds or thousands of records and are formatted based on some logic related to the algorithm being illustrated in the relevant papers. Furthermore, although not shown in this blog article, these figures are vector images which would be very nice to include in articles that will be printed. These two features, i.e., thousands of elements and user-defined logic, make it impractical to draw these figures manually using any image editor such as Inkscape. This blog article describes a simple tool that makes it very easy to generate all these figures.

Overview and Usage

Rasem uses Ruby code to describe images and generates standard SVG images that can be further edited or converted to EPS or PDF by popular applications including the open source Inkscape and Adobe Illustrator. Even if you're not a ruby developer, you can easily use Rasem as it's very simple to use.

Installation steps

  1. Install the Ruby compiler at
  2. Install Rasem gem by running the command
    gem install rasem

A simple example

Create a text file named 'test.rasem' with the following content
self.width = 100
self.height = 100
circle 20, 20, 5
circle 50, 50, 5
line 20, 20, 50, 50
The contents of the file are self explanatory. Now to generate an image out of this figure, all you need to do is to type the following command.
rasem test.rasem
or simply type
It will generate an SVG file named 'test.svg' which you can simply display in your Internet browser and is shown below.

You can find more examples at
The best feature of Rasem is that the file shown above is a fully functioning Ruby code. Which means you can include additional Ruby constructs and integrate it seamlessly with the visualized image. This is further shown in the following example which visualizes a simple graph.
nodes = [[10,10], [20,20], [10,20]]
radius = 3
with_style :fill=>"white", :stroke=>"black" do
  for node in nodes
    circle node.first, node.last, radius

line nodes[0].first, nodes[0].last, nodes[1].first, nodes[1].last,
The generated image is shown below.

How it works

Rasem relies on the dynamic features of Ruby which allows it to compile and run user-defined code on-the-fly. The line and circle commands shown above are nothing but function calls which are part of the API of Rasem. There are more functions in Rasem including rectangle, polygon, and polyline. These functions generate XML elements that conform to the SVG standard to draw the corresponding elements. That's how you can easily integrate your Ruby code into your figure to add user-defined logic which can affect both the generated figure elements and their style. You can also write a code that opens an external file, parses it, and generates the figure from this data. Rasem is open source and you're welcome to use it or contribute to it.

A Voronoi Diagram example

In this part, we show a more sophisticated example that we used to draw a Voronoi Diagram figure similar to the one shown below. We describe some tricks that can be used to generate similar figures.
The input file was generated using the Voronoi diagram operation in SpatialHadoop and it looks like the following snippet.
189.05078396875297,319.3714811786753    POLYGON ((194.43089260205164 325.5334210386081, 170.08384007573167 321.1403988020476, 198.24373337246874 315.0171135767282, 194.43089260205164 325.5334210386081))
3.678985871399587,429.815627730571    LINESTRING (-46.06193792012205 322.23283446689334, -46.06193792012205 322.23283446689334)
Each line has X and Y coordinates for the site followed by a WKT representation of a POLYGON or a LINESTRING that represents a closed or an open Voronoi region, respectively.
We start by providing the code that reads and parses those lines.
def read_regions(filename)
  lines = File.readlines(filename)
  regions = do |l|
    if l =~ /(.*),(.*)\t(.*)/
      {:x => $1.to_f, :y => $2.to_f, :wkt => $3}

To draw each entry, we provide the following simple code that draws the Voronoi region as either a polygon or a polyline.

def draw_region(svg, region)
  vertices = region[:wkt].scan(/[\-\d\.]+\s+[\-\d\.]+/).map{|coords| coords.split(" ").map(&:to_f)}
  if region[:wkt].index 'POLYGON'
    svg.polygon vertices
  elsif region[:wkt].index 'LINESTRING'
    svg.polyline vertices

Finally, the following function takes a file and draws all the regions in it. It also draws all sites as small circles.
def draw_vd(in_file, out_file, mbr)
  outs =, "w")
  mbr_width = mbr[2] - mbr[0]
  mbr_height = mbr[3] - mbr[1]
  regions = read_regions(in_file){:width=>mbr_width, :height => mbr_height}, outs) do
    rectangle mbr[0], mbr[1], mbr_width, mbr_height, :stroke=>"black", :fill=>:none
    # Sub-group 1: Points
    for region in regions
      circle region[:x], region[:y], 1
      draw_region svg, region

draw_vd "vd", "vd.svg", [0, 0, 500, 500] 

Now, let's make it more interesting. We would like to give different colors, as shown above, for safe and unsafe regions as described in the SpatialHadoop Voronoi diagram blog post. We define the following function that decides whether a site is safe or not.

def is_safe(region, mbr)
  site = [region[:x], region[:y]]
  vertices = region[:wkt].scan(/[\-\d\.]+\s+[\-\d\.]+/).map{|coords| coords.split(" ").map(&:to_f)}
  for vertex in vertices
    # Distance from the vertex to the site
    d_site = distance_to(site, vertex)
    # Minimum distance from the vertex to the MBR boundary
    d_mbr_min = [vertex[0] - mbr[0],
             mbr[2] - vertex[0],
             vertex[1] - mbr[1],
             mbr[3] - vertex[1]].min
    return false if d_mbr_min < d_site
  return true

def distance_to(p1, p2)
  dx = p1[0] - p2[0]
  dy = p1[1] - p2[1]
  return Math::sqrt(dx * dx + dy * dy)

We can now easily make a test before drawing each region to decide the color we use to draw each region. However, as a part of this tutorial we will do it in a different way. Instead of applying the styles in the Rasem code, we will apply them in Inkscape. We will use Rasem only to separate safe and unsafe regions. We will use the "group" feature in Rasem and SVG to put each set of regions together as shown in the following code.

def draw_vd(in_file, out_file, mbr)
  outs =, "w")
  mbr_width = mbr[2] - mbr[0]
  mbr_height = mbr[3] - mbr[1]

  regions = read_regions(in_file){:width=>mbr_width, :height => mbr_height}, outs) do
    rectangle mbr[0], mbr[1], mbr_width, mbr_height, :stroke=>"black", :fill=>:none
    # Group 1: Safe regions
    group do
      # Sub-group 1: Points
      group do
        for region in regions
          if is_safe(region, mbr)
            circle region[:x], region[:y], 1
      # Sub-group 2: Regions
      group(:fill => :none) do
        for region in regions
          if is_safe(region, mbr)
            draw_region svg, region

    # Group 2: non-safe regions
    group do
      # Sub-group 1: Points
      group do
        for region in regions
          if !is_safe(region, mbr)
            circle region[:x], region[:y], 1
      # Sub-group 2: Regions
      group(:fill => :none) do
        for region in regions
          if !is_safe(region, mbr)
            draw_region svg, region

The code is much longer than before but it is very simple. We create two big groups for safe and unsafe regions. Then, we create two subgroups within each one for sites and regions. This makes a total of four subgroups, namely, safe sites, safe regions, unsafe sites, and unsafe regions. This makes it easier to select apply a style to each individual group in Inkscape. This technique is preferable than applying the styles in the code because you can easily choose the most appealing colors in a GUI and get an instant feedback on your choice rather than modifying the code and regenerating the image.

Exporting standard SVG images is a very powerful feature in Rasem as it allows users to further edit the generated images in their favorite image editor such as Inkscape or Adobe Illustrator. This way, you can easily integrate different figures together, add some text, or some other annotations.

As usual the source code and sample files are provided on the following link.

Further readings

Rasem project page on Github
Voronoi diagram and Delaunay triangulation in SpatialHadoop
Example files

Thursday, March 31, 2016

Around the world in one hour! (revisit)

In this blog post, we revisit an earlier blog post about extracting data from OpenStreetMap Planet.osm file. We still use the same extraction script in Pigeon but we make it modular and easier to reuse. We make use of the macro definitions in Pig to extract common code into a separate file. In the following part, we first describe the OSMX.pig file which contains the reusable macros. After that, we describe how to use it in your own Pig script.


The osmx.pig file contains all the common code that is used to extract points, ways, or relations from an OSM file. It contains the following functions.


This macro extracts all the nodes from an OSM file. It returns a dataset that contains tuples of the following format.



This macro returns all ways in the file. Each way is returned as a series of line segments which connect two consecutive nodes on the way. It returns a dataset with tuples of the following format.
segment_idlongA generated unique ID for each segment
id1longThe ID of the starting node
latitude1doubleLatitude of the starting node
longitude1doubleLongitude of the starting node
id2longThe ID of the ending node
latitude2doubleLatitude of the ending node
longitude2doubleLongitude of the ending node
way_idlongThe ID of the way that contains this segment
tagsmap[(chararray)]All the tags of the way


This macro returns all ways in the file. However, unlike LoadOSMWaysWithSegments, it returns one tuple for each segment which contains the entire geometry of the way. Each tuple is formatted as follows.

way_idlongThe ID of the way as it appears in the OSM file
first_node_idlongThe ID of the first node in this way
last_node_idlongThe ID of the last node in this way
geombytearrayThe geometry of the way
tagsmap[(chararray)]The tags of the way as they appear in the OSM file


This macro returns all objects in the OSM file. Objects can be one of two cases:
  1. First level relations: This contains relations that contain only ways.
  2. Dangled ways: This contains ways that are not part of any relations.
The returned dataset does not contain second level relations such as relations that contain other relations. The format of the returned dataset is as follows.
object_idlongThe ID of either the relation or the way
geombytearrayThe geometry of the object
tagsmap[(chararray)]The tags of either the way or the relation as they appear in the OSM file


The script planet-extractor.pig provides an example that extracts the datasets that are available on the SpatialHadoop datasets page. The header of this file imports the 'osmx.pig' file as well as the required JAR libraries.

REGISTER spatialhadoop-2.4.jar;
REGISTER pigeon-0.2.1.jar;
REGISTER esri-geometry-api-1.2.jar;
REGISTER jts-1.8.jar;
IMPORT 'osmx.pig';

The next two lines extracts all nodes and writes them to a file.

all_nodes = LoadOSMNodes('$input');

STORE all_nodes INTO '$output/all_nodes.bz2';

This makes it much easier than earlier code where the extraction is interleaved with writing the output.
Similarly, the following few lines extract the road network and writes it to the output.

-- Extract road network
road_network = LoadOSMWaysWithSegments('$input');
road_network = FILTER road_network BY edu.umn.cs.spatialHadoop.osm.HasTag(tags,
road_network = FOREACH road_network GENERATE segment_id,
               id1, latitude1, longitude1,
               id2, latitude2, longitude2,
               way_id, edu.umn.cs.spatialHadoop.osm.MapToJson(tags) AS tags;
STORE road_network INTO '$output/road_network.bz2';

Although the code looks a little bit ugly, it only contains four statements. The first one extracts all the ways as segments using the LoadOSMWaysWithSegments macro. The second statement filters the segments that are related to the road network using the tags attribute. The third statement removes unnecessary columns and the fourth statement writes the output.

Similar to the road network, the next few lines extracts and stores the buildings dataset.

all_objects = LoadOSMObjects('$input');
 buildings = 
FILTER all_objects BY edu.umn.cs.spatialHadoop.osm.HasTag(tags,
STORE buildings INTO '$output/buildings.bz2';

The first statement extracts all the objects from the file. The second statement filters the buildings using the tags attribute. Finally, the third statement stores the output.


This work was partially supported by an AWS in Education Grant.

External Resources

Saturday, February 20, 2016

HadoopViz: Extensible Visualization of Big Spatial Data

With huge sizes of spatial data, a common functionality that users are looking for is to visualize this data to see how it looks like. This gives users the power of quickly exploring new datasets with huge sizes. For example, the video below summarizes 1 trillion points that represent the temperature of every 1 km2 on the earth surface on every day from 2009 to 2014 (total of six years).
This video consists of 72 frames, as one per month. These frames are put together in this video. While one can use a single machine to produce these 72 images, it might take up to 60 hours due to the huge size of the input.
In this blog post, we describe how to use HadoopViz, an extensible visualization framework based on SpatialHadoop, to visualize the same dataset in just three hours using a cluster of 10 machines.
Other than single-level images which are typically of low resolution, HadoopViz can also produce multilevel images where users can interactively zoom in and out to explore huge datasets with a lot of details. For example, the image below is a visualization of a 92GB dataset which represents all the objects extracted from OpenStreetMap dataset. You can pan and zoom in this image to view more details about a specific area.


In a nutshell, HadoopViz uses the parallelization power of MapReduce along with the efficiency of SpatialHadoop to partition the data into smaller parts, visualize each part separately into a smaller image, and then put these partial images together to produce the final image. HadoopViz builds on this idea and provides four key features that make it easy to use and very efficient.
  1. HadoopViz piggybacks data smoothing with visualization allowing it to smooth the data on-the-fly as the image is generated.
  2. HadoopViz automatically decides the best way to partition the data allowing it to scale to generate both small and large images efficiently.
  3. HadoopViz can also visualize multilevel images where users can freely pan and zoom into the image to interactively explore the huge dataset..
  4. Instead of customizing the algorithm for a specific use case, e.g., satellite data, HadoopViz provides an extensible implementation that can support a wide range of visualization types.
Below, we first describe how to generate the visualizations show above using HadoopViz, which ships with the recent version of SpatialHadoop. Then, we describe some technical details about the smoothing, partitioning, and extensibility features.

How to ...

... generate the temperature video

  1. You need to download and setup the most recent version of SpatialHadoop which ships with HadoopViz as its visualization package. Check this page for more details about setting up the most recent version of SpatialHadoop on both Hadoop 1.x and Hadoop 2.x.
  2. Download the temperature dataset you would like to visualize. The temperature dataset we used can be obtained from LP DAAC archive on this link.
    You can use this ruby script to download all the data for the six years if you have a good internet connection and enough storage on your machine. Run it using the following command:
    ruby hdf_downloader.rb time:2009.01.01..2014.12.31
  3. Once you have all the data, you can upload it to your HDFS using 'copyFromLocal' command. Let's assume the data is available at hdfs://user/hadoop/temperature
  4. To visualize the 72 frames, run the following SpatialHadoop command
    shadoop multihdfplot hdfs://user/hadoop/temperature combine:31 dataset:LST_Day_1km hdfs://user/hadoop/frames/ time:2009.01.01..2014.12.31
  5. The frames will be available in the output path hdfs://user/hadoop/frames. Download them using 'copyToLocal' command.
  6. Now, upload the frames to YouTube which will put them together into a video similar to the one shown above.

... generate the multilevel image

  1. Follow step 1 above to download and install SpatialHadoop, if you haven't done already.
  2. Download the 'All objects' dataset at the following link
  3. Upload the file to HDFS using the 'copyFromLocal' command. Let's assume it is uploaded to hdfs://user/hadoop/objects/
    NB: You don't have to decompress the file as SpatialHadoop can decompress it on the fly while visualizing. However, if you upload the compressed file, you need to keep the .bz2 extension to tell SpatialHadoop it is compressed.
  4. To generate a multilevel image with 11 levels similar to the one shown above, type the following command
    shadoop gplot hdfs://user/hadoop/objects -pyramid levels:11 hdfs://user/hadoop/multilevel shape:osm
  5. The generated image will be available at hdfs://user/hadoop/multilevel. Download it to your machine using the 'copyToLocal' command.
  6. To view the image in your browser, open the 'index.html' file available in the output directory.


In visualization, smoothing means the fuse of nearby records according to visualization logic to produce a correct result. For example, satellite datasets typically contain holes which are results of clouds that obstructs the view of the satellites. A smoothing function can recover these holes by estimating the missing values using simple interpolation techniques. The two figures below show an example of how the smoothing function can recover missing points.
Original data without smoothing
Data is smoothed using HadoopViz
HadoopViz support on-the-fly smoothing of the data as the visualization is done. This means that you can easily plug in a different smoothing function and regenerate the image without having to carry out the complex smoothing function as a separate step.


HadoopViz supports two ways of partitioning the data which affect the way it merges intermediate partial images. It can use either the default HDFS partitioning or the spatial partitioning that ships with SpatialHadoop.

Default HDFS Partitioning

By default, when you upload a file to HDFS, it is partitioned into equi-sized chunks of 128MB each. Spatial locations of records are not taken into account and nearby records will typically end up in two different partitions. This means that every partition would possibly cover the entire input space and we will end up overlaying intermediate images to produce the final image as shown below.
Overlay intermediate images

Spatial Partitioning

If we use the spatial partitioning that ships with SpatialHadoop, each partition would only contain data from a small limited space and we will end up stitching intermediate images as shown below.
Stitch intermediate images

Which partitioning technique is better?

While both techniques will end up producing the same final answer, the performance might be different. HadoopViz needs to automatically decide which one to use. First of all, if the data needs to be smoothed, then HadoopViz has to choose spatial partitioning as it is the only one that groups nearby records together in one partition before they can be fused.
If HadoopViz doesn't need to apply a smoothing function, then both techniques are applicable. According to the image size, There's an overhead between the partitioning and merging steps. The default HDFS partitioning is faster than spatial partitioning, but the overlay process is more time consuming than stitch due to the huge sizes of intermediate images. HadoopViz decides to go for spatial partitioning if the image size is huge as the cost of the overlay process becomes more and more time consuming.

Multilevel images

A multilevel image consists of a pyramid of fixed-size tiles, typically, each of size 256x256 pixels. The figure below shows an example of a three-level image with 1, 4, and 16 tiles in its three levels, aka, a pyramid of three levels.
A multilevel of three levels
A naive way to generate a multilevel image is to generate each tile independently using the (single-level) techniques shown above. However, this would require executing the single-level algorithm millions of times. Therefore, HadoopViz provides specialized multilevel visualization algorithms for multi-level images that take into consideration the pyramid structure of multi-level images. Similar to single-level visualization, HadoopViz supports two partitioning techniques, namely, default HDFS partitioning and pyramid partitioning.

Default HDFS Partitioning

If we use default HDFS partitioning, each partition might contain records from all over the input space. In this case, each machine plots all these records to all overlapping tiles in all pyramid levels. The generated tiles are considered partial images as multiple partitions might overlap the same tile. Thus, a final merge step will need to overlay all intermediate partial images for the same tile to produce the final image for that tile.

Pyramid Partitioning

The other option for HadoopViz is to first repartition the data so that all records that overlap with one tile go to one partition. Then, these records are visualized to generate the final image for that tile. No merging is needed here as each tile is only generated by one machine.

Which partitioning technique is better?

Again, there is no clear winner here. It all depends on how many tiles are generated. If only a few tiles are generated, then default HDFS partitioning is better as it only needs to merge a few images. However, if a huge number of tiles are generated, pyramid partitioning is better as it avoids altogether the need for merging intermediate tiles.
HadoopViz splits a huge pyramid into two parts, the top and the base of the pyramid. The top of the pyramid contains only a few tiles and is generated by the default HDFS partitioning technique, while the base contains too many tiles and is generated by the pyramid partitioning technique. The tiles are then put together to produce the final image without any extra processing.


While the above techniques can be customized for every visualization type, it would require a huge coding effort to build and maintain all these implementations. Therefore, HadoopViz proposes a visualization abstraction that is used to describe the visualization logic. This abstraction is then plugged into generic implementations of the above algorithms to produce the image efficiently at scale. In short, if you would like to visualize your own data in a new way, all you need to do is write a small class that extends an abstract class, and you're ready to go with both single-level and multilevel visualization techniques.
A new visualization type is defined by extending the base class Plotter. There are mainly five functions that you would like to implement for a new visualization type.

<S extends Shape> Iterable<S> smooth(Iterable<S> r)

This function takes a set of nearby records, fuses them together, and returns a new set of records. This function can be used to apply a user-specified smoothing logic.

Canvas createCanvas(int width, int height, Rectangle mbr)

This function initializes an empty canvas with the given size in width and height. It also associates this canvas with the given MBR in input space. Notice that Canvas can be virtually anything. We provide a simple abstract Canvas class as a skeleton.

void plot(Canvas layer, Shape shape)

The plot function updates the canvas layer by plotting the given shape on it. Users can define their own visualization logic for one shape based on the format of the shape and the canvas layer.

void merge(Canvas finalLayer, Canvas intermediateLayer)

The merge function merges two intermediate canvases. It updates the finalLayer by merging the intermediateLayer into it.

void writeImage(Canvas layer, DataOutputStream out, boolean vflip)

This writeImage function encodes the canvas layer into a standard image that can be displayed to the end user. The image is written to the given DataOutputStream which typically goes to an output file. If the vflip flag is set to true, the image should be vertically flipped before written to the output. The vflip flag is useful when the y-axis of the input is in a different direction than the final image. For example, in PNG images, the y-axis increases from bottom to top while in geographical coordinates, latitude increases from top to bottom.


This work was partially supported by an AWS in Education Grant.

Further References

  1. SpatialHadoop homepage:
  2. Ahmed Eldawy, Mohamed F. Mokbel and Christopher Jonathan "HadoopViz: A MapReduce Framework for Extensible Visualization of Big Spatial Data". In Proceedings of the 32nd IEEE International Conference on Data Engineering, IEEE ICDE 2016, Helsinki, Finland, May 16-20, 2016
  3. Ahmed Eldawy, Mohamed F. Mokbel and Christopher Jonathan "A Demonstration of HadoopViz: An Extensible MapReduce System for Visualizing Big Spatial Data". In Proceedings of the International Conference on Very Large Databases, VLDB 2015, Kohala Coast, HI, 2015

Wednesday, December 2, 2015

Voronoi diagram and Dealunay triangulation construction of Big Spatial Data using SpatialHadoop

Voronoi Diagram and Delaunay Triangulation

A very popular computational geometry problem is the Voronoi Diagram (VD), and its dual Delaunay Triangulation (DT). In both cases, the input is a set of points (sites). In VD, the output is a tessellation of the space into convex polygons, as one per input site, such that each polygon covers all locations that are closest to the corresponding site than any other site. In DT, the output is a triangulation, where each triangle connects three sites, such that the circumcirlce of each triangle does not contain any other sites. These two constructs are dual in a sense that each edge in the DT connects two sites that share a common edge in VD.

Computation Overhead

These two problems, DT and VD, are known to be computationally intensive. There are several O(n log n) algorithms for both but they are mainly single-machine main-memory algorithms and cannot scale to the current sizes of Big Spatial Data. For example, earth surface modeling is usually done via Delaunay triangulation of elevation data. With the increasing resolution of data collected through LiDAR, we might need to construct a Delaunay triangulation for a few trillions of points which cannot be loaded into one machine.

SpatialHadoop and CG_Hadoop

SpatialHadoop provides a computational geometry library, named CG_Hadoop, which includes MapReduce algorithms for both VD and DT. These algorithms can scale to trillions of points by distributing the computation over a Hadoop cluster with up-to thousands of machines. SpatialHadoop is able to construct a DT of 2.7 billion points extracted from OpenStreetMap in less than an hour using only 20 machines.

Single-machine Divide and Conquer (D&C) Algorithm

SpatialHadoop MapReduce algorithm is based on the  main-memory divide and conquer algorithm. This traditional algorithm, originally proposed by Guibas and Stolfi, divides the input sites into two subsets using a straight line, usually a vertical line, computes the DT recursively for each subset, and then merges the two partial DTs. This algorithm requires all data to fit in the main memory which cannot be done with big spatial datasets.

CG_Hadoop MapReduce Algorithm

A straight forward implementation of the D&C algorithm in MapReduce wouldn't work because of the overhead of the merge step. The final merge step has to be done on a single machine and it would require loading the whole final DT in its main memory, which takes us back to the limitation of the single-machine algorithm. Alternatively, CG_Hadoop proposes a novel approach which employs the following two techniques:
  1. Instead of partitioning the space into two, CG_Hadoop partitions the space using one of the SpatialHadoop partitioning methods, e.g., Grid or R-tree, so that it achieves a higher level of parallelism by assigning each partition to one of the machines.
  2. CG_Hadoop early detects VD portions that are final, i.e., will never be affected or used by subsequent merge process. In other words, CG_Hadoop detects the Voronoi regions, or Delaunay triangles, that might be affected by the merge step, and only sends these portions to the merge step.
These two techniques allow the construction process of both VD and DT to be much more scalable as each machine computes a partial VD one partition at a time. Then, each machine early detects parts of the partial VD that will not be affected by the subsequent merge steps, and early writes them to the output. Furthermore, the merge step becomes much faster and scalable as it deals with only a few non-final parts of the VD rather than all of it.
The figure below illustrates the main idea of the pruning technique. The input points are partitioned using a uniform grid into four partitions. The VD of each partition is computed on a separate machine in the map function. Each machine detects final Voronoi regions, marked in green, and writes them to the final output. Only non-final regions, marked in red, are sent to the next merge step.

Local VD step. Each mapper computes a partial VD for one partition and writes final Voronoi Regions (green) to the output. Non-final regions (red) are sent to the next step for merging

How to Use it

The above algorithm ships with the recent version of SpatialHadoop available on github. First, you need to spatially partition the input set of points using the 'index' command. For example
shadoop index input.points shape:point input.grid sindex:grid
Then, you will call the DT operation as follows:
shadoop dt input.grid output.dt shape:point
It launches a single MapReduce job that generates the DT of the input set of points.

More Details

The VD algorithm runs in four steps, namely, partitioning, local VD, vertical merge, and horizontal merge.

1. Partitioning

The partitioning step uses any space partitioning technique supported by SpatialHadoop, such as a uniform grid index.

2. Local VD

In the local VD step, each partition is assigned to one machine which computes the VD for all points in that partition. We use the traditional Guibas and Stolfi D&C algorithm as the size of each partition is small enough to fit in the main memory. We can also use a more advanced algorithm, such as Dwyer's algorithm. We leave this as an exercise for interested readers. Upon the computation of the local VD, each machine detects final Voronoi regions, removes them from the VD, and writes them to the final output. To detect final regions, we use a very simple, yet efficient, rule that is derived from the VD definition. Simply, a Voronoi region is final if it is not possible to have any other point closer to its generator site. The following figure shows two examples of final and non-final regions.
A final Voronoi region. The nearest possible site outside the partition is still farther away than the generator site of that region
A non-final Voronoi region. There could be a site (S2) in another partition that is closer to some (shaded) areas in that Voronoi region.
To find all non-final regions quickly, without exhaustively applying this rule to all regions, we run a breadth-first search (BFS) from the outer regions which expands only to non-final regions. The trick is that all non-final regions form a contiguous region towards the exterior of the VD. In this way, a VD of 1.4 million regions can be pruned by applying the rule to only 7K regions.

3. Vertical Merge

In this step, partial VDs which are on the same column of the partition grid are merged together into vertical strips. This step can be done in parallel as each machine is assigned all partial VDs in one vertical strip. Notice that each machine receives only non-final regions from the assigned VDs are processed by that machine which reduces the memory and computation overhead. After merging them, the pruning rule is applied again to further remove final regions that will not be processed later, as shown in the figure below.
The vertical merge step. Each vertical strip is processed by a reducer which merges all partial VDs in that strip. Final regions in the strip (blue) are directly written to the output while non-final regions (red) are sent to the next step for the final merge step.

4. Horizontal Merge

In this final step, all partial VDs computed for vertical strips are sent to one machine which merges them to compute the final answer. Again, only non-final regions are processed at this step which makes it possible to run on a single machine. As shown in the figure below, all resulting regions are written to disk as there is no further merge steps.
The final horizontal merge step. All resulting Voronoi regions are written to the final output as there are no further merge step.

The Voronoi regions written by the three steps comprise the final VD, as shown in the figure below.
The final results. Green, blue and black regions are the ones written by the local VD, vertical merge, and horizontal merge steps, respectively.

Performance Evaluation

We run the above algorithm on Amazon EC2 using a real dataset extracted from OpenStreetMap. We compare to the standard D&C algorithm running on a machine with 1TB of memory. We have access to such a powerful machine through Minnesota Supercomputing Institute (MSI). CG_Hadoop is 30 times faster than a single machine with 1 billion points. Furthermore, the single machine fails to process the 2.7 billion dataset while CG_Hadoop finishes the computation in around 40 minutes.


The work in Voronoi diagrams was done with Nadezda Weber while she was studying in University of Minnesota. She is currently a graduate student in Oxford University.
The experiments were run on an Amazon EC2 cluster and were supported by AWS in Education Grant award.

External Links

Monday, November 30, 2015

Reducing the memory footprint of the spatial join operator in Hyracks

This is the fourth blog post in a series that describes how to build an efficient spatial join Hyracks operator in AsterixDB. You can refer to the previous posts below:
  1. An Introduction to Hyracks Operators in AsterixDB
  2. Your first Hyracks operator
  3. A Hyracks operator for plane-sweep join

Scope of this post

In the third post, I described how to implement an efficient plane-sweep join algorithm in a Hyracks operator. That implementation simply caches all data frame, or simply records, in memory before running the plane-sweep algorithm. As the input datasets go larger, this algorithm might require a huge memory footprint which is not desirable with the big data that is handled by AsterixDB. In this blog post, I will describe how to improve the previous operator to run with a limited memory capacity.

Main idea

The plane-sweep join algorithm can be considered as the two-dimensional version of the sort-merge join algorithm. After sorting the two lists, you synchronously iterate over them with two pointers and you always advance the pointer that points to the smaller key. This means that you do not really need to look at the two whole datasets for the algorithm to work. You only need to look at a small window of records in both datasets. While the two datasets can be arbitrarily large, this window can remain, in most cases, of a fixed size. The improvement that we're going to describe in this post revolves around this idea of keeping only the records in this small window in memory.


As usual, the first step is to refactor the current code to prepare it for the new logic.


First, we extracted the CachedFrameWriter class to a separate file and improved it to work with a limited memory capacity. We used a CircularQueue which keeps all the queue elements in a fixed size array. We also added a new feature to the CachedFrameWriter that automatically clears all cached frames before the mark point. Thus, whenever the plane-sweep algorithm puts a mark, the CachedFrameWriter translates it as an indication that earlier frames will never be accessed in the plane-sweep join algorithm.
public void mark() {
  // This mark indicates that we do not need to get back beyond this point
  // We can shrink our queue now to accommodate new data frames
  this.markFrame = this.currentFrame = 0;
  this.markRecord = this.currentRecord;
We also added a couple of more functions that will be used by the plane-sweep join operator.
public boolean canGrowInMemory() {
  return !cachedFrames.isFull() && !reachedEndOfStream;
/**Whether there are more frames to be read from input*/
public boolean isComplete() {
  return reachedEndOfStream;

/**Are there more cached frames to iterate over*/
public boolean noMoreImmediatelyAvailableRecords() {
    return this.currentFrame >= this.cachedFrames.size();

Incremental PlaneSweepJoin Helper Class

We modified the helper class that actually performs the plane-sweep join algorithm to be stateful so that it can run incrementally. It keeps track of its state and can continue from where it stopped. This allows the plane-sweep operator to run some iterations of the plane-sweep join which will advance the pointers on the two datasets. Based on this advance, it can evict some data frames from memory to accommodate for more frames. Then, it can continue from where it stopped. This also required changing the plane-sweep method to be an instance method, rather than class method, so that it can access the state stored in the instance. The PlaneSweepJoin class now has the following instance variables.
public class PlaneSweepJoin {
  private IHyracksTaskContext ctx;
  private CachedFrameWriter[] datasets;
  private IFrameWriter outputWriter;
  private ITuplePairComparator rx1sx1;
  private ITuplePairComparator rx1sx2;
  private ITuplePairComparator sx1rx2;
  private IPredicateEvaluator predEvaluator;
  private FrameTupleAppender appender;

   * An enumerated type that stores the current state of the spatial join
   * algorithm so that it can continue from where it stopped.
  public enum SJ_State {
    /** The spatial join algorithm did not start yet */
    /** The current record in dataset 0 is active
      * and needs to be tested with more records in dataset 1 */
    /** The current record in dataset 1 is active
        and needs to be tested with more records in dataset 0 */
    /** No active records */
    /** The spaital join algorithm has finished. No more work to do. */
  /** Current state of the spatial join algorithm */
  private SJ_State sjState;
Till this point, we didn't do any serious changes to the actual plane-sweep join method. It works as before and it just keeps its state in instance variables instead of local variables.

Remodeled PlaneSweepJoinOperatorDescriptor

Now this is the important part of this blog post where we modify the plane-sweep operator to make it more memory friendly. The main modification is the timing when the plane-sweep function is called. Previously, it was called once after the two inputs are cached completely. Currently, it will be called more often to do part of the join and clean up some memory space. So, we modified the CachedFrameWriter#nextFrame method to call plane-sweep join operator as soon as the in-memory buffer is full.
public void nextFrame(ByteBuffer buffer) throws HyracksDataException {
  // Store this buffer in memory for later use
  ByteBuffer copyBuffer = ctx.allocateFrame(buffer.capacity());
  FrameUtils.copyAndFlip(buffer, copyBuffer);
  if (cachedFrames.isFull()) {
    // run the plane-sweep algorithm in case it can free some buffer entries
    try {
      if (owner.getPlaneSweepJoin().getState() == PlaneSweepJoin.SJ_State.SJ_FINISHED) {
    } catch (InterruptedException e) {
   // TODO If after running the plane-sweep, we still cannot find empty entries,
   // we should start spilling records to disk.
   if (cachedFrames.isFull())
     throw new HyracksDataException("Memory full");
Notice that there are two instances of CachedFrameWriter, once for each input dataset. Therefore, the PlaneSweepJoin#planesweepJoin method can be called by the two CachedFrameWriters. The method has to be internally modified to synchronize these calls. In the code shown above, a CachedFrameWriter calls the planesweepJoin method when its memory buffer is full. In this case, this method waits until the other CachedFrameWriter also has its memory buffer full, and then it actually carries out the spatial join algorithm, releases some memory buffers, and returns the call back to the sender. This part is shown in the following code snippet.
public synchronized void planesweepJoin(CachedFrameWriter fullDataset)
        throws HyracksDataException, InterruptedException {
  // No more work to do
  if (sjState == SJ_State.SJ_FINISHED)
  CachedFrameWriter otherDataset =
      fullDataset == datasets[0] ? datasets[1] : datasets[0];
  if (otherDataset.canGrowInMemory()) {
  // ... Perform the actual spatial join logic here
  if ((datasets[0].isComplete() &&
       datasets[0].noMoreImmediatelyAvailableRecords()) ||
      (datasets[1].isComplete() &&
    sjState = SJ_State.SJ_FINISHED;
  if (sjState == SJ_State.SJ_FINISHED) {
    appender.flush(outputWriter, true);
  // To wake up the other thread
In the code above, when the planesweepJoin method is called, the sender tells which CachedFrameWriter it is. This tells the planesweepJoin method whether it needs to block and wait until the other dataset has a memory buffer full, or it can go ahead with the join logic rightaway. Since there are only two datasets, the planesweepJoin method determines the other dataset and decides to wait until the other dataset has a memory full (i.e., cannot grow in-memory anymore). This condition will match when the first CachedFrameWriter is full but the second one is not. When the second CachedFrameWriter is also full, it will not match this condition and will carry out the actual spatial join code. At the very end, after the whole method is finished, the notify() method is called to wake up the other sleeping thread. That sleeping thread will return immediately as it knows that the method which woke it up has already performed the spatial join operator.
In this execution pattern, the spatial join can be thought as a cooperative task where any of the two threads can perform it, whenever it is possible. This is usually better than spawning a third thread that keeps joining records which are inserted by the two other threads. Also, it fits the architecture of AsterixDB better as the underlying framework is the one responsible of creating threads.


To test whether this code runs or not, we execute the spatial join operation against two large datasets and on a limited memory that cannot hold the two full datasets. In this case, it has to execute multiple runs of the spatial join algorithm to get the work done. The test is similar to previous one but it has an additional parameter to limit the memory usage.
PlaneSweepJoinOperatorDescriptor join =
  new PlaneSweepJoinOperatorDescriptor(spec,
    new X1X1ComparatorI(), new X1X2ComparatorI(),
    new X1X2ComparatorI(), outputDesc, 10,
    new SpatialOverlapPredicateI());
The sixth parameter (with value 10) tells the PlaneSweepJoinOperatorDescriptor that it has at most 10 memory buffers to cache the input datasets while the dataset needs at least 20 buffers to be completely cached in memory. Of course in a real environment, we would probablly have much more memory buffer.
In the new test case found in the download section, I also added a new case when the two inputs are not sorted. It uses the existing Hyracks sort operator to sort both inputs before feeding them to the spatial join operation. You can find the source code of this new test case in the downloads section below.


As usual, all the source code of this blog post can be found here.