Castor project's earthquakes map with D3js

The Castor project is a submarine natural gas storage facility located in front of the eastern Iberian Peninsula coast. The idea is to store the Algerian gas inside an old oilfield cavity. At least this is what I understood (sorry, geologists).
Somehow, when the facility started working, a series of earthquakes have started to occur. At the beginning, the platform owners said it wasn't related to their activity, but now everybody agrees that it is, and the activity has stopped, but not the earthquakes.

I didn't find a nice map about the earthquakes epicenters, so I thought that D3js would be a good option to do it.

The star represents the aproximate platform location, and the circles, the epicenters. It's easy to see why they are related to the platform activity.

The animated map is at my bl.ocks.org page. The explanations are in Catalan, but basically say the same as here.

Getting the data

The data about the significant earthquakes around Catalonia can be found at the Catalan Geologic Institute web site, but the format for the reports is not very convenient to get the data, so I made this python script to get it:

# -*- coding: utf-8 -*-
#http://jramosgarcia.wordpress.com/2013/10/01/que-es-el-proyecto-castor/
import urllib
import json
import dateutil.parser
import datetime

def get_data(url):
    filehandle = urllib.urlopen(url)
    html = filehandle.read()
    filehandle.close()

    lines = html.splitlines()

    for i in range(len(lines)):
    
        if lines[i].find('Latitud') > 0:
            lat = lines[i].strip().split(" ")[1].replace("º","")
        if lines[i].find('Longitud') > 0:
            lon = lines[i].strip().split(" ")[1].replace("º","")   
        if lines[i].find('mol del dia') > 0:
            date = lines[i + 1].strip().replace(" >/div<","")
        if lines[i].find('Hora origen') > 0:
            hour = lines[i].strip().split(" ")[4]                   
        if lines[i].find('Magnitud') > 0:
            magnitude = lines[i+1].strip().split(" ")[0]

    date_array  = date.split("/")
    hour_array = hour.split(":")

    date_time = datetime.datetime(int(date_array[2]),int(date_array[1]),int(date_array[0]),int(hour_array[0]), int(hour_array[1]))

    data = {'lat':lat, 'lon':lon, 'datetime': date_time.isoformat(), 'magnitude': magnitude}
    return json.dumps(data)

if __name__ == "__main__":
    url_list = [
    'http://www.igc.cat/web/gcontent/ca/sismologia/sismescomact/comhistcat/20130910175510/comact.html',
    'http://www.igc.cat/web/gcontent/ca/sismologia/sismescomact/comhistcat/20130913095842/comact.html',
    'http://www.igc.cat/web/gcontent/ca/sismologia/sismescomact/comhistcat/20130918142943/comact.html',
    'http://www.igc.cat/web/gcontent/ca/sismologia/sismescomact/comhistcat/20130920104607/comact.html',
    'http://www.igc.cat/web/gcontent/ca/sismologia/sismescomact/comhistcat/20130924091301/comact.html',
    'http://www.igc.cat/web/gcontent/ca/sismologia/sismescomact/comhistcat/20130925125029/comact.html',
    'http://www.igc.cat/web/gcontent/ca/sismologia/sismescomact/comhistcat/20130929084140/comact.html',
    'http://www.igc.cat/web/gcontent/ca/sismologia/sismescomact/comhistcat/20130929192416/comact.html',
    'http://www.igc.cat/web/gcontent/ca/sismologia/sismescomact/comhistcat/20130930005900/comact.html',
    'http://www.igc.cat/web/gcontent/ca/sismologia/sismescomact/comhistcat/20130930051316/comact.html',
    'http://www.igc.cat/web/gcontent/ca/sismologia/sismescomact/comhistcat/20131001045206/comact.html',
    'http://www.igc.cat/web/gcontent/ca/sismologia/sismescomact/comhistcat/20131001055709/comact.html',
    'http://www.igc.cat/web/gcontent/ca/sismologia/sismescomact/comhistcat/20131002121626/comact.html',
    'http://www.igc.cat/web/gcontent/ca/sismologia/sismescomact/comhistcat/20131002232928/comact.html',
    'http://www.igc.cat/web/gcontent/ca/sismologia/sismescomact/comhistcat/20131003012732/comact.html',
    'http://www.igc.cat/web/gcontent/ca/sismologia/sismescomact/comhistcat/20131003031301/comact.html',
    'http://www.igc.cat/web/gcontent/ca/sismologia/sismescomact/comhistcat/20131004113222/comact.html',
    'http://www.igc.cat/web/gcontent/ca/sismologia/sismescomact/comhistcat/20131004120323/comact.html'
    ]

    f = open("data.json","w")
    f.write("[")
    json_data = ""
    for url in url_list:
     json_data = json_data + get_data( url ) + ", "
       
    f.write(json_data[0:-2])
    f.write("]")
    f.close()
   

• The web pages are in the list, since the interesting reports have to be chosen one by one. It would be nice to have a better way to do it.
• Then, the get_data function just parses the text in the way that all the reports are parsed properly.The data is stored in a json file to make easier it's use from D3js.

Using D3js to visualize the data

I used this example by Mike Bostock to create the background  map. The tiles origin has been changed because the example tiles are not available at this zoom level, and to have more points of interest to situate the earthquake locations.

This is the code:

<script>
/**
Based on Mike Bostock's http://bl.ocks.org/mbostock/4150951
*/
var width = 960,
    height = 500;

var projection = d3.geo.mercator()
    .center([0.5491, 40.4942])
    .scale(20000);

var path = d3.geo.path()
    .projection(projection);

var tile = d3.geo.tile()
    .scale(projection.scale() * 2 * Math.PI)
    .translate(projection([0, 0]))
    .zoomDelta((window.devicePixelRatio || 1) - .5);

var svg = d3.select("body").append("svg")
    .attr("width", width)
    .attr("height", height);

  var tiles = tile();

  var defs = svg.append("defs");

var magnitudeScale = d3.scale.linear().domain([2,5]).range([5,30]);
d3.json("data.json", function(error, locations) {

  svg.append("g")
      
    .selectAll("image")
      .data(tiles)
    .enter().append("image")
      .attr("xlink:href", function(d) { return "http://" + ["a", "b", "c", "d"][Math.random() * 4 | 0] + ".tiles.mapbox.com/v3/examples.map-vyofok3q/" + d[2] + "/" + d[0] + "/" + d[1] + ".png"; })
      .attr("width", Math.round(tiles.scale) + 1)
      .attr("height", Math.round(tiles.scale) + 1)
      .attr("x", function(d) { return Math.round((d[0] + tiles.translate[0]) * tiles.scale); })
      .attr("y", function(d) { return Math.round((d[1] + tiles.translate[1]) * tiles.scale); });

  svg.append("g")
      .append('path')
      .attr("d","m 0,0 -8.47858,-5.22254 -8.31623,5.47756 2.34696,-9.67752 -7.77927,-6.21653 9.92909,-0.75852 3.50829,-9.31953 3.78972,9.20873 9.94748,0.45679 -7.58687,6.44982 z")
      .attr("stroke","black")
      .attr("stroke-width",2)
      .style("fill", d3.rgb(90, 90, 90))
      .attr("transform", "translate("+projection([0.66879, 40.33503])[0]+","+projection([0.66879, 40.33503])[1]+")");
  
  var locationsG = svg.append("g");
  

  function addLocation(loc){
    locationsG.append('circle')
      .attr('class','location')
      .attr("r", 5)
      .attr("cx", projection([loc.lon, loc.lat])[0])
      .attr("cy", projection([loc.lon, loc.lat])[1])
      .style("fill", d3.rgb(255, 0, 0).darker(2))
      .style("opacity", 0.8)
      .transition()
      .duration(1000)
      .attr("r", magnitudeScale(loc.magnitude))
      .transition()
      .delay(2000)
      .duration(2000)
      .style("opacity",0.3);

    locationsG
      .append("text")
      .text(loc.magnitude)
      .attr("x", projection([loc.lon, loc.lat])[0] - 10)
      .attr("y", projection([loc.lon, loc.lat])[1] + 5)
      .attr("font-family", "sans-serif")
      .attr("font-size", "12px")
      .attr("fill", "black")
      .style("opacity",0)
      .transition()
      .duration(1000)
      .style("opacity",1)
      .transition()
      .delay(2000)
      .duration(2000)
      .style("opacity",0);
  }

  //addLocation({"lat": "40.43", "magnitude": "2.7", "lon": "0.7", "datetime": "2013-10-09T06:43:16"});
  
  var intDate = new Date("2013-09-10T00:00:00Z");
  var maxDate = new Date("2013-10-04T00:00:00Z");
  var usedLocations = new Array();

  var dateTitle = svg
    .append("text")
    .attr("id", "dataTitle")
    .text(intDate.toLocaleDateString())
    .attr("x", 70)
    .attr("y", 20)
    .attr("font-family", "sans-serif")
    .attr("font-size", "20px")
    .attr("fill", "black");

  var interval = setInterval(function() {

    dateTitle.text(intDate.toLocaleDateString());

    for (i = 0; imaxDate){
        clearInterval(interval);
      }

    }
    
    intDate.setDate(intDate.getDate() + 1);
  }, 1000);
  
});

•  Lines 30 and 48 are the most important ones to create the tile map. In this case, no zoom or pan is set, so the code is quite simple.
• The line 57 creates the star indicating the platform location. I made the star using Inkscape and captured the code using the tool to see the generated svg. To situate the symbol on the map, a transform attribute is used, translating the symbol to the pixel calculated projecting the location coordinates.
• To add the earthquake epicenters:
• The function addLocation (line 68) adds the circle and the text indicating the earthquake magnitude. To locate them, the projection function is used again. Two transitions are used to make the circle grow and then get brighter. The same for making the text disappear.
• An interval is created (line 118) to change the date (one day per second). For every iteration, the date is evaluated and the earthquakes occured during this day plotted using the addLocation function. I didn't find a "more D3js" solution. The earthquakes already plotted are stored in an array to avoid plotting them more than once.
  The date on the map is also changed for every iteration.

Some related posts and pages:

• D3.geo.tile
• D3 World maps 
• Let's make a map
• Three little circles
• Que es el proyecto Castor? (in Spanish). Very good explanation about the Castor project, with some illustrations for the non-Spanish speakers.
• Quakes rattle Spain, offshore gas storage project blamed