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