search

drewfustin / isocronut

For a given geospatial location, calculate an isochrone (same time) contour around it.
97 stars 44 forks source link

Doubt on isochrones for 2 hours #5

Closed wolframteetz closed 8 years ago

wolframteetz commented 8 years ago

Something's weird about the isochrones... the 90 and 120 minute isochrones are almost identical to the west although there is a highway in that direction. I strongly doubt that these are correct. Need to investigate the cause.

(html file with 60,90,120min isochrones and 48 directional vectors)

Unfortunately, my daily limit's gone again --

    <!DOCTYPE html>
    <html>
    <head>
    <meta name="viewport" content="initial-scale=1.0, user-scalable=no">
    <meta charset="utf-8">
    <title>Isochrone</title>
    <style>
      html, body, #map-canvas {
        height: 100%;
        margin: 0px;
        padding: 0px
      }
    </style>

    <script src="https://maps.googleapis.com/maps/api/js?v=3.exp&signed_in=true"></script>

    <script>
    function initialize() {
      var mapOptions = {
        zoom: 8,
        center: new google.maps.LatLng(48.1089197,11.451958),
        mapTypeId: google.maps.MapTypeId.TERRAIN
      };

      var map = new google.maps.Map(document.getElementById('map-canvas'), mapOptions);

      var marker = new google.maps.Marker({
          position: new google.maps.LatLng(48.1089197,11.451958),
          map: map
      });

      var isochrone;

      var isochroneCoords = [
    new google.maps.LatLng(48.69084,11.45411), 
new google.maps.LatLng(48.6821922,11.611889), 
new google.maps.LatLng(48.6821922,11.611889), 
new google.maps.LatLng(48.64269,11.78697), 
new google.maps.LatLng(48.4982456,11.7956948), 
new google.maps.LatLng(48.564644,11.9831995), 
new google.maps.LatLng(48.5040218,12.0576007), 
new google.maps.LatLng(48.4573868,12.1469998), 
new google.maps.LatLng(48.28024,12.07553), 
new google.maps.LatLng(48.3974082,12.7495727), 
new google.maps.LatLng(48.2167837,12.1330527), 
new google.maps.LatLng(48.1665,12.12216), 
new google.maps.LatLng(48.1101058,12.1113211), 
new google.maps.LatLng(48.0780006,11.9669667), 
new google.maps.LatLng(48.0053526,12.0136777), 
new google.maps.LatLng(47.9346692,12.014565), 
new google.maps.LatLng(47.81522,12.19909), 
new google.maps.LatLng(47.8640182,12.0093171), 
new google.maps.LatLng(47.7172118,11.9563338), 
new google.maps.LatLng(47.7505202,11.85938), 
new google.maps.LatLng(47.74573,11.76186), 
new google.maps.LatLng(47.7168545,11.719957), 
new google.maps.LatLng(47.75063,11.5962), 
new google.maps.LatLng(47.72323,11.52728), 
new google.maps.LatLng(47.7058469,11.409518), 
new google.maps.LatLng(47.62628,11.35382), 
new google.maps.LatLng(47.6341305,11.233276), 
new google.maps.LatLng(47.5975289,11.1884435), 
new google.maps.LatLng(47.6645158,11.0703333), 
new google.maps.LatLng(47.7863276,11.0604693), 
new google.maps.LatLng(47.8739757,11.0244206), 
new google.maps.LatLng(47.9118491,11.0575774), 
new google.maps.LatLng(47.9139271,10.8490461), 
new google.maps.LatLng(47.8855514,10.6467416), 
new google.maps.LatLng(48.0051291,10.5949358), 
new google.maps.LatLng(48.1005133,10.561436), 
new google.maps.LatLng(48.17436,10.68676), 
new google.maps.LatLng(48.2272487,10.7888558), 
new google.maps.LatLng(48.3265549,10.6511387), 
new google.maps.LatLng(48.4659587,10.5135426), 
new google.maps.LatLng(48.5297439,10.6963314), 
new google.maps.LatLng(48.5617904,10.6791623), 
new google.maps.LatLng(48.56507,10.91982), 
new google.maps.LatLng(48.5628425,11.0627245), 
new google.maps.LatLng(48.57973,11.14955), 
new google.maps.LatLng(48.5626401,11.2630624), 
new google.maps.LatLng(48.5561561,11.3652545), 

      ];

      isochrone = new google.maps.Polygon({
        paths: isochroneCoords,
        strokeColor: '#000',
        strokeOpacity: 0.5,
        strokeWeight: 1,
        fillColor: '#0000FF',
        fillOpacity: 0.1
      });

      isochrone.setMap(map);

          var isochrone2;

      var isochroneCoords2 = [
    new google.maps.LatLng(49.01559,11.63351), 
new google.maps.LatLng(48.9352388,11.7886699), 
new google.maps.LatLng(49.0347825,12.0631613), 
new google.maps.LatLng(48.9854354,12.1959296), 
new google.maps.LatLng(48.77459,12.23209), 
new google.maps.LatLng(48.7879195,12.4686735), 
new google.maps.LatLng(48.80049,12.83738), 
new google.maps.LatLng(48.5510303,12.556838), 
new google.maps.LatLng(48.3889649,12.4925531), 
new google.maps.LatLng(48.31527,12.64295), 
new google.maps.LatLng(48.1847363,12.4313311), 
new google.maps.LatLng(48.1005422,12.4657629), 
new google.maps.LatLng(48.0077528,12.4679199), 
new google.maps.LatLng(47.87637,12.69482), 
new google.maps.LatLng(47.7629182,12.6453054), 
new google.maps.LatLng(47.7070091,12.4768688), 
new google.maps.LatLng(47.6197486,12.3883036), 
new google.maps.LatLng(47.5761214,12.2262449), 
new google.maps.LatLng(47.5683328,12.062927), 
new google.maps.LatLng(47.6892024,11.7704106), 
new google.maps.LatLng(47.6446218,11.7433702), 
new google.maps.LatLng(47.3445665,11.6890377), 
new google.maps.LatLng(47.3445665,11.6890377), 
new google.maps.LatLng(47.6825636,11.5765039), 
new google.maps.LatLng(47.6825636,11.5765039), 
new google.maps.LatLng(47.4433478,11.1744672), 
new google.maps.LatLng(47.4617723,11.0053989), 
new google.maps.LatLng(47.4035774,10.8806966), 
new google.maps.LatLng(47.6316531,10.8201512), 
new google.maps.LatLng(47.61059,10.71506), 
new google.maps.LatLng(47.6366,10.54275), 
new google.maps.LatLng(47.6644577,10.3069488), 
new google.maps.LatLng(47.7724165,10.2581077), 
new google.maps.LatLng(47.789715,9.8857952), 
new google.maps.LatLng(47.9675891,9.948535), 
new google.maps.LatLng(48.069478,9.9515912), 
new google.maps.LatLng(48.223185,10.1784321), 
new google.maps.LatLng(48.3497101,10.0558781), 
new google.maps.LatLng(48.5405915,9.9226157), 
new google.maps.LatLng(48.6223985,10.2391536), 
new google.maps.LatLng(48.6641069,10.3632798), 
new google.maps.LatLng(48.7518434,10.4708748), 
new google.maps.LatLng(48.8600826,10.5008528), 
new google.maps.LatLng(48.93889,10.71853), 
new google.maps.LatLng(48.9532292,10.9188556), 
new google.maps.LatLng(48.8405607,11.218338), 
new google.maps.LatLng(49.2316261,11.2074204), 
new google.maps.LatLng(49.09484,11.45046), 

      ];

      isochrone2 = new google.maps.Polygon({
        paths: isochroneCoords2,
        strokeColor: '#000',
        strokeOpacity: 0.5,
        strokeWeight: 1,
        fillColor: '#000',
        fillOpacity: 0.25
      });

      isochrone2.setMap(map);

  var isochrone3;

      var isochroneCoords3 = [
    new google.maps.LatLng(49.26565,11.45254), 
new google.maps.LatLng(49.25586,11.68265), 
new google.maps.LatLng(49.216316,11.8367523), 
new google.maps.LatLng(49.186964,12.1129601), 
new google.maps.LatLng(49.10618,12.33497), 
new google.maps.LatLng(49.02192,12.52515), 
new google.maps.LatLng(48.9199723,12.6965073), 
new google.maps.LatLng(48.8169289,12.8017088), 
new google.maps.LatLng(48.6804397,12.9679205), 
new google.maps.LatLng(48.5328439,13.0575411), 
new google.maps.LatLng(48.3720899,13.1193429), 
new google.maps.LatLng(48.2452541,13.1805614), 
new google.maps.LatLng(48.097362,13.1838851), 
new google.maps.LatLng(47.9448083,13.1724472), 
new google.maps.LatLng(47.7967697,13.1164647), 
new google.maps.LatLng(47.65502,13.03983), 
new google.maps.LatLng(47.5251713,12.9483501), 
new google.maps.LatLng(47.4035357,12.8251358), 
new google.maps.LatLng(47.2857166,12.6596103), 
new google.maps.LatLng(47.1893922,12.4765612), 
new google.maps.LatLng(47.1427419,12.2762211), 
new google.maps.LatLng(47.037923,12.100363), 
new google.maps.LatLng(46.990808,11.9637515), 
new google.maps.LatLng(46.9702609,11.594396), 
new google.maps.LatLng(47.0028056,11.5042219), 
new google.maps.LatLng(46.8684818,11.3173944), 
new google.maps.LatLng(46.9917225,11.0124301), 
new google.maps.LatLng(47.0211632,10.7792317), 
new google.maps.LatLng(47.1202927,10.6461183), 
new google.maps.LatLng(47.1202927,10.6461183), 
new google.maps.LatLng(47.4098581,10.2797448), 
new google.maps.LatLng(47.4656903,10.1044633), 
new google.maps.LatLng(47.5211163,9.9694213), 
new google.maps.LatLng(47.6475674,9.8293321), 
new google.maps.LatLng(47.776013,9.7674612), 
new google.maps.LatLng(47.9226998,9.7532823), 
new google.maps.LatLng(48.0951473,9.7901525), 
new google.maps.LatLng(48.2494496,9.7289785), 
new google.maps.LatLng(48.39828,9.76904), 
new google.maps.LatLng(48.5430757,9.7913518), 
new google.maps.LatLng(48.6799436,9.9370948), 
new google.maps.LatLng(48.825928,10.1004936), 
new google.maps.LatLng(48.9176446,10.2259365), 
new google.maps.LatLng(49.0195622,10.3680885), 
new google.maps.LatLng(49.0872864,10.5405129), 
new google.maps.LatLng(49.1751,10.77501), 
new google.maps.LatLng(49.2322124,10.9798272), 
new google.maps.LatLng(49.2533173,11.2388746), 

      ];

      isochrone3 = new google.maps.Polygon({
        paths: isochroneCoords3,
        strokeColor: '#000',
        strokeOpacity: 0.5,
        strokeWeight: 1,
        fillColor: '#000',
        fillOpacity: 0.25
      });

      isochrone3.setMap(map);
    }

    google.maps.event.addDomListener(window, 'load', initialize);
    </script>

    </head>
    <body>
    <div id="map-canvas"></div>
    </body>
    </html>
wolframteetz commented 8 years ago

Looking at the code carefully

rmax = [1.25 * duration] * number_of_angles  # rmax based on 75 mph speed

Looking at the map, I'd say you imposed a 75 KMH speed limit Which is quite unrealistic for Autobahn - so in Germany, you will get perfect circles for longer distances ;)

Cannot check if altering this to 290kmh (for my motorbike ;) (no really, 150 kmh, so 2.5 is fine) works, will do tomorrow once my quota is back

drewfustin commented 8 years ago

Let me know if switching it up to a higher limit works. I was well aware that this was completely not an ideal way to handle this.

Also, I'm going to post below a different take on this whole project I worked on for a different company. It uses Mapbox, is entirely JS, allows you to click on a point on a map and have it draw isochrones dynamically. Mapbox just released their directions API to the public (I had a private beta earlier). In doing so, or because of something else, my code now throws Ajax errors and I haven't fixed it yet (I'm no JS guy). So, it doesn't work. But maybe it's useful. It's certainly quicker and more robust if fixed.

drewfustin commented 8 years ago

The basics of this algorithm: for a central point on a map (clicked by the user), throw a bunch of points around it and test the duration to each of these points. Use isolines turf.js to draw contours using these durations as values to be contoured over. Quick, easy. But now with Ajax errors.

    <!DOCTYPE html>
    <html>
    <head>
    <meta name="viewport" content="initial-scale=1.0, user-scalable=no">
    <meta charset="utf-8">
    <title>Isochrone</title>

    <script src="https://api.tiles.mapbox.com/mapbox.js/v2.2.2/mapbox.js"></script>
    <script src="turf.min.js" charset="utf-8"></script>

    <link href='https://api.tiles.mapbox.com/mapbox.js/v2.2.2/mapbox.css' rel='stylesheet' />
    <style>
      body { margin:0; padding:0; }
      #map { position:absolute; top:0; bottom:0; width:100%; }
    </style>

    </head>
    <body>

    <div id="map"></div>

    <script>
    var accessToken = 'GetYourAPITokenFromMapboxAndPutItHere';

    L.mapbox.accessToken = accessToken;

    var map = L.mapbox.map('map', 'mapbox.streets')
        .setView([31.2, -99.3], 6);

    map.on('click', function(e) {
        create_isochrones(map, [e['latlng']['lat'], e['latlng']['lng']]);
    });

    function create_isochrones(map, origin) {

        map.setView(origin, 8);

        var well = {
            type: 'Feature',
            geometry: {
                type: 'Point',
                coordinates: [
                  origin[0],
                  origin[1]
                ]
            },
            properties: {}
        };

        L.mapbox.featureLayer(turf.flip(well)).addTo(map);

        var m = [50, 150, 300];
        var s = [50, 100, 100];

        for (var i in m) {
            if (i > 0) {
                var circle = turf.buffer(well, m[i], 'miles');
                var extent = turf.extent(circle);
                var additions = turf.random('points', s[i], {bbox: extent});
                points['features'].push.apply(points['features'], additions['features']);
            } else {
                var circle = turf.buffer(well, m[i], 'miles');
                var extent = turf.extent(circle);
                var points = turf.random('points', s[i], {bbox: extent});
            };
        };

        var deferreds = [];
        var destinations = {
            type: 'FeatureCollection',
            features: []
        };

        for (i = 0; i < points.features.length; i++) {
            var url = 'http://api.tiles.mapbox.com/v4/directions/mapbox.driving/' + well.geometry.coordinates[1].toString() + ',' + well.geometry.coordinates[0].toString() + ';' + points.features[i].geometry.coordinates[1].toString() + ',' + points.features[i].geometry.coordinates[0].toString() + '.json?access_token=' + L.mapbox.accessToken;
            deferreds.push(
                $.ajax({
                  url: url
                })
                .success(function(data) {
                    var min_dur = 9999999;
                    for (r = 0; r < data['routes'].length; r++) {
                        min_dur = Math.min(min_dur, data['routes'][r]['duration'] / 60);
                    };
                    if (data['destination']) {
                        var d = {
                            type: 'Feature',
                            geometry: {
                                type: 'Point',
                                coordinates: [
                                    data['destination']['geometry']['coordinates'][1],
                                    data['destination']['geometry']['coordinates'][0]
                                ]
                            },
                            properties: {
                                'z': min_dur
                            }
                        };
                    destinations['features'].push(d);
                    };
                })
                .fail(function() {
                  alert("Ajax failed to fetch data")
                })
            );
        }

        $.when.apply($, deferreds).done(
            function() {
                var breaks = [30, 60, 90, 120];
                var colors = {};
                colors[breaks[0]] = 'rgb(255, 0, 0)';
                colors[breaks[1]] = 'rgb(255, 255, 0)';
                colors[breaks[2]] = 'rgb(0, 255, 0)';
                colors[breaks[3]] = 'rgb(0, 0, 255)';
                var isolined = turf.isolines(destinations, 'z', 200, breaks);

                for (i = isolined.features.length - 1; i >= 0; i--) {
                    var polygon_options = {
                        color: '#000',
                        opacity: 0.25,
                        weight: 1,
                        fillColor: colors[isolined.features[i].properties.z],
                        fillOpacity: 0.25
                    };
                    L.polygon(isolined.features[i].geometry.coordinates, polygon_options).addTo(map);
                }
            }
        );

    };

    </script>

    </body>
    </html>
wolframteetz commented 8 years ago

Somewhere I still have a nice Radial Basis Function network fitting written in Octave, need to transfer that to JS and include it to make really nice interpolations ;)

Before I saw your blog I had actually thought about doing exactly what you have just posted - only that I didn't have early access to the Mapbox api. I'm not a JS expert either, but I surely know more than I do about python.. too busy right now but that looks like that way to go. I will definitely take a look at it later.

wolframteetz commented 8 years ago

Fixed & Pull requested the example

wolframteetz commented 8 years ago

Confirmed, if you move to any country with metric system, better replace 75 "mph" with 150 "kph" ;)

drewfustin commented 8 years ago

Cool. That's good to know. You're amazing. Appreciate the help!