Calculate areas in which a person can reach a certain point within a certain time limit

area.php

<?php

mb_internal_encoding("utf8");

$askedG=0;

//takes as argiment a csv file with format "id,lat,lng"
$files = file($argv[1]);
foreach($files as $line) {
  echo "Working $line\n";
  $line = explode(",",trim($line));  
  $geoj = getArea(doubleval($line[1]),doubleval($line[2]));
  file_put_contents("../data/geojson/".$line[0].".geojson",$geoj);
	$f = fopen ( "../data/geojson/.".$line[0].".geojson.gz", 'w' );
	fwrite ( $f,  gzencode ( $geoj , 9 ));
	fclose ( $f );

  sleep(10);
  if ($askedG>2000) break;
}

echo "\nAsked all $askedG.\n";


function getArea($lat,$lng) {
  global $askedG;
  $source = array($lat,$lng);

  $c = array("type"=>"FeatureCollection",
    "features"=> array());

  $points=array();
  for ($j=0;$j<32;$j+=1) {
    $start = $j%4==0 ? 10 : ($j%2==0 ? 50 : 110 );
    for ($i=$start;$i<=210;$i+=20)
      $points[] = getOffset($source,$i,$j);
  }

  $asked=0;
  $blocked = array();
  $points2 = array();

  for ($k=10;$k<=210;$k+=20) {
    $dest = array();
    for ($i=0;$i<count($points);$i++)
      if ($points[$i][2]==$k && !in_array($points[$i][3],$blocked)) 
        $dest[]=$points[$i];

    $size = count($dest);
    if ($size==0)
      continue;
    if ($asked+$size>=100) {
      echo "Sleeping... (asked $asked) \n";
      sleep(11);
      $asked=0;
    }

    $asked += $size;
    $askedG += $size;
    $destT="";
    for ($i=0;$i<count($dest);$i++)
      $destT[]=$dest[$i][0].",".$dest[$i][1];
    $destT = implode("|",$destT);

    $res = file_get_contents("https://maps.googleapis.com/maps/api/distancematrix/json?origins=$lat,$lng&destinations=$destT&mode=driving&language=en&sensor=false&key=[here be google api key]");
    sleep(1);

    if (!$res) die("error with google");
    $res = json_decode($res);
    if (!$res->status || $res->status!="OK") die("error with google");
    for ($i=0;$i<count($res->rows[0]->elements);$i++) {
      if ($i>=count($dest)) break;
      if ($dest[$i][4]=$res->rows[0]->elements[$i]->status!="OK") continue;
      $dest[$i][4]=ceil($res->rows[0]->elements[$i]->duration->value/60);
      $points2[]=$dest[$i];
      if ($dest[$i][4]>=1.5*60)
        $blocked[]=$dest[$i][3];
    }
  }

  $ways=array();
  $paths=array(array(),array(),array());

  for ($i=0;$i<count($points2);$i++) {
   
 /*   $c["features"][]=array("type" => "Feature",
      "properties" => array("time"=>$points2[$i][4],"dist"=>$points2[$i][2],"bear"=>$points2[$i][3]),
      "geometry" => array(
        "type" => "Point",
        "coordinates" => array($points2[$i][0],$points2[$i][1]))
    );
*/

    if(!$ways[$points2[$i][3]])
      $ways[$points2[$i][3]]=array();
    $ways[$points2[$i][3]][]=$points2[$i];
  }


//  return cleanGeoJson(json_encode($c));

  for ($j=0;$j<32;$j+=1) {
    $way=$ways[$j];
    if (count($way)==0) 
      continue;

    $p=false;
    $p0=false;
    $p1=false;
    $p2=false;
    usort($way,"waysort");
    for ($i=0;$i<count($way);$i++) {
      if ($way[$i][4]>=90 && !$p2) {
        if ($p==false)
          $endd = $way[$i][2]/$way[$i][4]*90;
        else
          $endd = $p[2]+($way[$i][2]-$p[2])/($way[$i][4]-$p[4])*(90-$p[4]);
        $p2=getOffset($source,$endd,$j);
      }
      if ($way[$i][4]>=60 && !$p1) {
        if ($p==false)
          $endd = $way[$i][2]/$way[$i][4]*60;
        else
          $endd = $p[2]+($way[$i][2]-$p[2])/($way[$i][4]-$p[4])*(60-$p[4]);
        $p1=getOffset($source,$endd,$j);
      }
      if ($way[$i][4]>=30 && !$p0)  {
        if ($p==false)
          $endd = $way[$i][2]/$way[$i][4]*30;
        else
          $endd = $p[2]+($way[$i][2]-$p[2])/($way[$i][4]-$p[4])*(30-$p[4]);
        $p0=getOffset($source,$endd,$j);
      }
      $p=$way[$i];
    }

    if ($p0)
      $paths[0][]=$p0;
    else if ($p[4]>30) {
      $endd = $p[2]/$p[4]*30;
      $paths[0][]=getOffset($source,$endd,$j);
    } else
      $paths[0][]=$p;

    if ($p1)
      $paths[1][]=$p1;
    else if ($p[4]>60) {
      $endd = $p[2]/$p[4]*60;
      $paths[1][]=getOffset($source,$endd,$j);
    } else
      $paths[1][]=$p;

    if ($p2)
      $paths[2][]=$p2;
    else if ($p[4]>90) {
      $endd = $p[2]/$p[4]*90;
      $paths[2][]=getOffset($source,$endd,$j);
    } else
      $paths[2][]=$p;
  }

  $poli=array();
  for ($i=0;$i<count($paths[0]);$i++)
    $poli[]=array($paths[0][$i][0],$paths[0][$i][1]);
  $poli[]=$poli[0];

  if ($poli)
  $c["features"][]=array("type" => "Feature",
    "properties"=> array("time"=>30),
    "geometry" => array(
      "type" => "Polygon", 
      "coordinates" => array($poli))
  );

  $poli=array();
  for ($i=0;$i<count($paths[1]);$i++)
    $poli[]=array($paths[1][$i][0],$paths[1][$i][1]);
  $poli[]=$poli[0];

  if ($poli)
  $c["features"][]=array("type" => "Feature",
    "properties"=> array("time"=>60),
    "geometry" => array(
      "type" => "Polygon", 
      "coordinates" => array($poli))
  );

  $poli=array();
  for ($i=0;$i<count($paths[2]);$i++)
    $poli[]=array($paths[2][$i][0],$paths[2][$i][1]);
  $poli[]=$poli[0];

  if ($poli)
  $c["features"][]=array("type" => "Feature",
    "properties"=> array("time"=>90),
    "geometry" => array(
      "type" => "Polygon", 
      "coordinates" => array($poli))
  );

  return cleanGeoJson(json_encode($c));
}

function waysort($a,$b) {
  if ($a[4] == $b[4])
        return $a[2] <= $b[2] ? -1 : 1;
    return $a[4] < $b[4] ? -1 : 1;
}

function getOffset($source, $km, $way) {
  $r = 6372.797;

  $deg = deg2rad($way/32*360);

	$lat1 = deg2rad($source[0]);
	$lng1 = deg2rad($source[1]);
  $lat1s = sin($lat1);
  $lat1c = cos($lat1);

  $dr = $km/$r;
  $drs = sin($dr);
  $drc = cos($dr);

  $lat2 = asin($lat1s*$drc+$lat1c*$drs*cos($deg));
  $lng2 = $lng1+atan2(sin($deg)*$drs*$lat1c, $drc-$lat1s*sin($lat2));

  $lat2 = rad2deg($lat2);
  $lng2 = rad2deg($lng2);

	return array($lat2, $lng2, $km, $way);
}

function cleanGeoJson($data) {
	return mb_ereg_replace('(-?\d{1,3}\.\d{0,5})\d*\,(-?\d{1,3}\.\d{0,5})\d*','\2,\1',$data);
}
?>

And yeah, it's not pretty, but for a night's work it's ok.