-1

I have an online app I am developing to search for historical wheat collection sites matched to climate http://www.wheat-gateway.org.uk/climate_search_1.php?ord=4&cns=108,114&ctrl_r=1//522/891&ctrl=1

I have the collection sites (some 173,000 and the data of the wheats collected from them) and their climate records and I have a climate record for a set of longitude/latitude points at 15 minute intervals for all land.

As you can see currently I am displaying points from the climate data set as heatmap in Google Maps but I want to convert these into lines, points and polygons as geoJSON file for display in Google maps viz https://developers.google.com/maps/documentation/javascript/datalayer partly for speed (I hope) and partly because for this purpose heatmaps distort when zoomed out and don't work visually zoomed right down either.

I have found https://assemblysys.com/php-point-in-polygon-algorithm/ which looks good to strip out points in search results array once I have found polygons - but first of all how to find them (and lines and points)? Any suggestions?

yours Andy Forbes

Andrew Forbes
  • 81
  • 1
  • 4
  • 1
    If you try to solve this in pure PHP you're going to have a Bad Time™. You're going to want to put all that data into a database with spatial indexes and GIS capabilities. I would suggest postgres and [postGIS](https://postgis.net/). – Sammitch Feb 25 '19 at 18:48
  • 1
    Hi and welcome to StackOverflow. Please post code samples rather than links - you need a minimal, complete, and verifiable example – Yvonne Aburrow Feb 25 '19 at 19:08
  • an error on your page - `Uncaught TypeError: map.getStreetView is not a function at get_res (climate_search_1.php?ord=4&cns=108,114&ctrl_r=1//522/891&ctrl=1:1088)` – Professor Abronsius Feb 26 '19 at 08:48
  • Yep, shows error "map.getStreetView is not a function" but I remove no longer any operating Streetview, if I move get different errors even if they also are not fatal. – Andrew Forbes Feb 26 '19 at 10:21
  • Thanks for your advice Sammitch, you are maybe right but I plunged on anyway and seem to be getting results in php. http://wheat-gateway.org.uk/json_clean.php which is about half the set used in example. Code too long to include here (how else do I post?) but basically I create an array of all possible points in the area, peel them off one by one N>S, W>E, dump if no actual point and examine where there is a point to see if it has neighbours, follow along path of neighbours (anti-clockwise) till get back to start point, check result and remove from search set, continue to next remaining. – Andrew Forbes Feb 26 '19 at 10:39

1 Answers1

0

so basically I cycle though the array of points that can be used as heatmap in Google maps from W>E and N>S and each new one found check if lone point or has neighbours. If has neighbours then trace around to see if part of polygon or line. When vertices back to start point are found run a check on remainder of points to see if they are inside the vertices, remove any found and the vertices' points and put vertices points and lone points into arrays that will output geoJSON file. Move to next point left in points until none left.

The costly part in terms of time is checking if points are inside a polygon. I found https://assemblysys.com/php-point-in-polygon-algorithm/ to be most accurate non-library set of functions to do this - most vertices found. I tried a variety of others including https://geophp.net/ which I found over twice as slow as "php-point-in-polygon-algorithm" funciton though I can't install GEOS on my server as recommended. However I found library phpgeos https://github.com/mjaschen/phpgeo did a much quicker job.

So here is result with both heatmap and points(as diamonds),lines and polygons - you can see on this area heatmap is much quicker.

http://www.wheat-gateway.org.uk/climate_search_both.php?ord=4&cns=110,113,116,124,131,135,150,153,157,161&ctrl_r=1//66/535,2//-5/31&ctrl=1,2

but on a smaller area the extra time is acceptable http://www.wheat-gateway.org.uk/climate_search_both.php?ord=4&cns=108,114&ctrl_r=1//291/1135,2//-5/12&ctrl=1,2

so my solution is to find out how many points in set and go with heatmap over 3000 points. I have managed to tune radius of points in heatmap depending on zoom level so the presentation doesn't give typical distortion when zoomed out eg http://www.wheat-gateway.org.uk/climate_search.php?ord=4&cns=all&ctrl_r=1//1061/1361,2//23.043/27.043&ctrl=1,2&targ=107723

Probably the most tricky bit is finding if there are holes in any polygon found, see the "rings" function here.

// requires library php-geos

$lines = array();
$pts = array();
$polys = array();
$holes = array();
//array to put results in
$cnt_p = 1;


$div = 0.25;
//gap between points in grid
$motion = 8;

$points=array (
  '-9.50:38.75' => '-9.50 38.75',
  '-9.25:38.75' => '-9.25 38.75',
  '-9.25:39.00' => '-9.25 39.00',
  '-9.25:39.25' => '-9.25 39.25',
  '-9.00:38.50' => '-9.00 38.50',
  '-9.00:38.75' => '-9.00 38.75',
  '-9.00:39.00' => '-9.00 39.00',
  '-9.00:39.25' => '-9.00 39.25',
  '-9.00:39.50' => '-9.00 39.50',
  '-9.00:39.75' => '-9.00 39.75',
  '-8.75:38.00' => '-8.75 38.00',
  '-8.75:38.25' => '-8.75 38.25',
  '-8.75:38.50' => '-8.75 38.50',
  '-8.75:38.75' => '-8.75 38.75',
  '-8.75:39.00' => '-8.75 39.00',
  '-8.75:39.25' => '-8.75 39.25',
  '-8.75:39.50' => '-8.75 39.50',
  '-8.75:39.75' => '-8.75 39.75',
  '-8.75:40.00' => '-8.75 40.00',
  '-8.75:40.25' => '-8.75 40.25',
  '-8.75:40.50' => '-8.75 40.50',
  '-8.50:38.25' => '-8.50 38.25',
  '-8.50:38.50' => '-8.50 38.50',
  '-8.50:38.75' => '-8.50 38.75',
  '-8.50:39.00' => '-8.50 39.00',
  '-8.50:39.25' => '-8.50 39.25',
  '-8.50:39.50' => '-8.50 39.50',
  '-8.50:39.75' => '-8.50 39.75',
  '-8.25:38.50' => '-8.25 38.50',
  '-8.25:38.75' => '-8.25 38.75',
  '-8.25:39.00' => '-8.25 39.00',
  '-8.25:39.25' => '-8.25 39.25',
  '-8.25:39.50' => '-8.25 39.50',
  '-8.25:39.75' => '-8.25 39.75',
  '-8.00:37.25' => '-8.00 37.25',
  '-8.00:38.50' => '-8.00 38.50',
  '-8.00:38.75' => '-8.00 38.75',
  '-8.00:39.00' => '-8.00 39.00',
  '-8.00:39.25' => '-8.00 39.25',
  '-8.00:39.50' => '-8.00 39.50',
  '-7.75:38.75' => '-7.75 38.75',
  '-7.75:39.00' => '-7.75 39.00',
  '-7.75:39.25' => '-7.75 39.25',
  '-7.75:39.50' => '-7.75 39.50',
  '-7.75:39.75' => '-7.75 39.75',
  '-7.50:38.75' => '-7.50 38.75',
  '-7.50:39.00' => '-7.50 39.00',
  '-7.50:39.25' => '-7.50 39.25',
  '-7.50:39.50' => '-7.50 39.50',
  '-7.50:39.75' => '-7.50 39.75',
  '-7.50:40.00' => '-7.50 40.00',
  '-7.25:39.00' => '-7.25 39.00',
  '-7.25:39.25' => '-7.25 39.25',
  '-7.25:39.50' => '-7.25 39.50',
  '-7.25:39.75' => '-7.25 39.75',
  '-7.25:40.00' => '-7.25 40.00',
  '-7.00:38.00' => '-7.00 38.00',
  '-7.00:38.50' => '-7.00 38.50',
  '-7.00:39.50' => '-7.00 39.50',
  '-7.00:39.75' => '-7.00 39.75',
  '-7.00:40.00' => '-7.00 40.00',
  '-7.00:40.25' => '-7.00 40.25',
  '-7.00:40.50' => '-7.00 40.50',
  '-7.00:40.75' => '-7.00 40.75',
)

function scope_it() {
    global $scope, $points, $re, $pts, $lines, $polys, $holes;
    if ($points) {
        $re = key($points);
        do_res($re);
        //start search picking N>S and W>E
    }
}

function do_res($re) {
    // start enquiry
    global $points, $results;
    $results = array();
    array_push($results, $points[$re]);
    $targ = $re;
    clock($targ);
    //set off search
}


function clock($targ) {
    // sniff along from neigbour to neigbour to find lines and polygons
    global $poss, $points, $re, $starter, $motion, $pts, $results;
    get_poss($targ);



//get box of points to search for neighbours
    $got = NULL;
    for ($i = 0; $i <= $motion - 1; ++$i) {
        $pos = $i + $starter;
        $find = isset($points[$poss[$pos]]);
        if ($find !== FALSE) {
            $got = $poss[$pos];
            array_push($results, $points[$poss[$pos]]);
            $i = $motion + 1;
            //got a neighbour, stop search 
        } elseif ($pos == $motion) {
            $starter = -$i;
        }
    }



    if ($got === NULL) {
        array_push($pts, $points[$re]);
        $sn2 = explode(":", $re);
        // no neighbours, it must be a point
        unset($points[$re]);
        $starter = 1;
        scope_it();
    } elseif ($got == $re) {
        $starter = 1;
        examine($results);
        // go to check if result is a polygon or a line
    } else {
        // end of object not reached, carry on sniffing along line till back to start
        $starter = $pos + 5;
        if ($starter > 8) {
            $starter = $starter - 8;
        }
        //change start point in next search in order to avoid coming back to previous point in the line (unless all the way around
        clock($got);
    }
}

function get_poss($targ) {
    global $poss, $div, $points;

    $poss = array();
    $sn = explode(":", $targ);

    $sn = explode(":", $targ);
    $poss[1] = number_format($sn[0] - $div, 2) . ":" . $sn[1];
    $poss[2] = number_format($sn[0] - $div, 2) . ":" . number_format($sn[1] + $div, 2);
    $poss[3] = $sn[0] . ":" . number_format($sn[1] + $div, 2);
    $poss[4] = number_format($sn[0] + $div, 2) . ":" . number_format($sn[1] + $div, 2);
    $poss[5] = number_format($sn[0] + $div, 2) . ":" . $sn[1];
    $poss[6] = number_format($sn[0] + $div, 2) . ":" . number_format($sn[1] - $div, 2);
    $poss[7] = $sn[0] . ":" . number_format($sn[1] - $div, 2);
    $poss[8] = number_format($sn[0] - $div, 2) . ":" . number_format($sn[1] - $div, 2);
    //compile array of possible neighbours to a point
}

function get_dia($quo) {
    global $div;
    $poss = array();
    $sn = explode(" ", $quo);
    $div2 = $div / 2;

    $poss[4] = "[" . $sn[0] . ", " . number_format($sn[1] - $div2, 3) . "]";
    $poss[1] = "[" . number_format($sn[0] - $div2, 3) . ", " . $sn[1] . "]";
    $poss[2] = "[" . $sn[0] . ", " . number_format($sn[1] + $div2, 3) . "]";
    $poss[3] = "[" . number_format($sn[0] + $div2, 3) . ", " . $sn[1] . "]";
    $poss[0] = "[" . $sn[0] . ", " . number_format($sn[1] - $div2, 3) . "]";
    //get cordinate of diamond around points to convert them to small polygon/diamond

    return implode(",", $poss);
}

function examine($polygon) {
    // run through polygon check
    global $points, $lines, $polys, $cnt_p, $n, $s, $w, $e, $in_points, $pts, $geofence;
    $poly = 0;
    $in_points = array();
    $n = -180;
    $s = 180;
    $w = 90;
    $e = -90;

    $geofence = new Polygon();
    foreach ($polygon as $k => $v) {
        $sn = explode(" ", $v);
        $geofence->addPoint(new Coordinate($sn[1], $sn[0]));
        $key = $sn[0] . ":" . $sn[1];
        $in_points[$key] = $v;
        unset($points[$key]);
        unset($search[$sn[0]][$sn[1]]);
    }

    $area = $geofence->getArea();

    if ($area > 0) {
        foreach ($points as $key => $point) {

            $pt = explode(" ", $point);
            $Point = new Coordinate($pt[1], $pt[0]);
            if (($geofence->contains($Point)) == 1) {
                $in_points[$key] = $points[$key];
                // put points into array to be searched for holes (rings)
                unset($points[$key]);
                //find area square if new polygon
                if ($pt[0] > $n)
                    $n = $pt[0];
                if ($pt[0] < $s)
                    $s = $pt[0];
                if ($pt[1] < $w)
                    $w = $pt[1];
                if ($pt[1] > $e)
                    $e = $pt[1];
            }
        }

        if (count($in_points) > 9) {
            rings($polygon, $in_points, $cnt_p);
        }
        $polys[$cnt_p] = $polygon;
        ++$cnt_p;
    } else {
        $break = array_unique($polygon);
        foreach ($break as $k => $v) {
            array_push($pts, $v);
        }
    }
    scope_it();
}

function rings($polygon, $in_points, $cnt_p) {
    global $n, $s, $w, $e, $div, $ring_s, $holes, $geofence;
// create array same area sqaure as polygon
    $ring_s = array();
    for ($a = $w; $a <= $e; $a = $a + $div) {
        for ($b = $n; $b >= $s; $b = $b - $div) {
            $val = number_format($b, 2) . " " . number_format($a, 2);
            $key = number_format($b, 2) . ":" . number_format($a, 2);
            $ring_s[$key] = $val;
        }
    }

// remove points in new array outside of polygon
    foreach ($ring_s as $key => $point) {
        $pt = explode(":", $key);
        $Point = new Coordinate($pt[1], $pt[0]);
        if (($geofence->contains($Point)) == 0) {
            unset($ring_s[$key]);
        }
    }

    foreach ($polygon as $k => $v) {
        $sn = explode(" ", $v);
        $key = $sn[0] . ":" . $sn[1];
        $ring_s[$key] = $v;
    }

    foreach ($in_points as $k => $v) {
        $find = isset($ring_s[$k]);
        //deduct from new array ponits found in those to be sorted, leaves any holes in polygon
        if ($find !== FALSE) {
            unset($ring_s[$k]);
        }
    }
    $cn_r = count($ring_s);
    if ($cn_r > 0) {
        $holes[$cnt_p] = array();
        scope_r($cnt_p);
    }
}

function scope_r($cnt_p) {
    global $ring_s, $in_points, $poss, $motion, $re, $results, $starter, $previous;
    // cycle through remains of the new array W>E and N>S looking for ones that are neigbouring vertices of holes in points to de sorted. When one is found jump onto that vertices and find the other points of internal polygon/edge.
    if ($ring_s) {
        $key = key($ring_s);
        $val = reset($ring_s);
        get_poss($key);
        $got = NULL;

        for ($i = 1; $i <= $motion; ++$i) {
            $find = isset($in_points[$poss[$i]]);
            if ($find !== FALSE) {
                $got = $poss[$i];
                $pos = $i;
                $i = $motion + 1;
            }
        }
        if ($got !== NULL) {
            $results = array();
            $re = $got;
            array_push($results, $in_points[$re]);
            $starter = $pos + 4;
            if ($starter > 8) {
                $starter = $starter - 8;
            }
            clock_r($cnt_p, $got);
            // jump to search set of points
        } else {
            unset($ring_s[$key]);
            scope_r($cnt_p);
        }
    }
}

function clock_r($cnt_p, $targ) {
    global $poss, $starter, $in_points, $motion, $re, $results, $holes, $ring_s;
    get_poss($targ);

//get box of points to search for neighbours
    $got = NULL;
    for ($i = 0; $i <= $motion - 1; ++$i) {
        // search through neighbours for internal hole edges
        $pos = $i + $starter;
        $find = isset($in_points[$poss[$pos]]);
        if ($find !== FALSE) {
            $got = $poss[$pos];
            array_push($results, $in_points[$got]);
            $i = $motion + 1;
        } elseif ($pos == $motion) {
            $starter = -$i;
        }
    }

    if ($got == $re) {
        $starter = 1;
        $poly = $results;
        array_push($holes[$cnt_p], array_reverse($poly));
        $geo2 = new Polygon();
        foreach ($poly as $k => $v) {
            $sn = explode(" ", $v);
            $geo2->addPoint(new Coordinate($sn[1], $sn[0]));
        }

        foreach ($ring_s as $key => $point) {
            $pt = explode(":", $key);
            $Point = new Coordinate($pt[1], $pt[0]);
            if (($geo2->contains($Point)) == 1) {
                unset($ring_s[$key]);
            }
        }
//when whole of internal polygon found ($re is start/end point), add it to array of holes per main polygon, reversing array order as these should be clockwise and main polygon anti-clockwise. Delete points internal to hole searching array and find any other holes left to record in that array
        scope_r($cnt_p);
    } else {
        $starter = $pos + 5;
        if ($starter > 8) {
            $starter = $starter - 8;
        }
        clock_r($cnt_p, $got);
    }
}

$starter = 1;
// start position for searching through neighbours
scope_it();


the array shown here outputs http://wheat-gateway.org.uk/json_map_final.php?ord=4&cns=109&ctrl_r=1//591/923,2//-5/31&ctrl=1,2 which in turn creates http://www.wheat-gateway.org.uk/climate_search.php?ord=4&cns=109&ctrl_r=1//591/923,2//-5/31&ctrl=1,2

Andrew Forbes
  • 81
  • 1
  • 4