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);
}
?>

I wrote this piece of code to generate polygons around a marker on the map showing in what area people would be able to reach that marker within 30, 60 and 90 mins. One result can be seen here:
https://gist.github.com/yurukov/380d9b8a61ea140d99eb

The algo is not as good as it could be. It works as follows: First it calculates a list of evenly distributed radial points around the center. Then it groups them by distance and calls Google Distance API to measure the time one needs to reach that point. Should the time for a point exceed 90, I don't process the rest from the angle/direction. Then it groups them by direction and attempts to place a point for 30, 60 and 90 mins based on the times it takes to reach the precomputed points on that one line. When we connect these final points, we get a polygon.

I tried improving the algorithm by geocoding the addresses the Distance API returns. Thus I hoped to adjust the precomputed points to the actual location where a visitor would reach. This caused more problems however, Especially in areas close to the sea like London and Istanbul.

The polygons can be improved to show splines so that the areas are more curved and easier to read. This should best be done in real time in JS as it would require quite a few points in the file itself. D3.js would be a good choice here.

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