Search code examples
phpgeospatial

Get circle polygon circumference points (latitude and longitude)


I am passing a center point (latitude, longitude) with a radius to this function to get the circumference points (latitude, longitude) so that I make a geometry from them and store it in MYSQL for me to do ST_Contains function later on.

This is getting me an oval shape instead of a circle when drawing them on Google Maps.

public static function convert($center, $radius, $numberOfSegments = 360)
{
    $n = $numberOfSegments;
    $flatCoordinates = [];
    $coordinates = [];
    for ($i = 0; $i < $n; $i++) {
        $flatCoordinates = array_merge($flatCoordinates, static::offset($center, $radius, 2 * pi() * $i / $n));
    }
    $flatCoordinates[] = $flatCoordinates[0];
    $flatCoordinates[] = $flatCoordinates[1];
    for ($i = 0, $j = 0; $j < count($flatCoordinates); $j += 2) {
        $coordinates[$i++] = array_slice($flatCoordinates, $j, 2);
    }

    return [
        'type' => 'Polygon',
        'coordinates' => [array_reverse($coordinates)]
    ];
}

public static function toRadians($angleInDegrees = null)
{
    return $angleInDegrees * pi() / 180;
}

public static function toDegrees($angleInRadians = null)
{
    return $angleInRadians * 180 / pi();
}

public static function offset($c1, $distance, $bearing)
{
    $lat1 = static::toRadians($c1[1]);
    $lon1 = static::toRadians($c1[0]);
    $dByR = $distance / 6378137; // distance divided by 6378137 (radius of the earth) wgs84
    $lat = asin(
        sin($lat1) * cos($dByR) +
            cos($lat1) * sin($dByR) * cos($bearing)
    );
    $lon = $lon1 + atan2(
        sin($bearing) * sin($dByR) * cos($lat1),
        cos($dByR) - sin($lat1) * sin($lat)
    );
    return [static::toDegrees($lon), static::toDegrees($lat)];
}

I do not really need it to be as accurate in terms of the earth curvature with taking into consideration WGS84, as the radius is quite small (25m-150m).


Solution

  • I'm guessing by the variable names that you got your function from this question but you forgot one line of the code, specifically the line that was causing that person problems.

    In addition, your array code in the first function was quite more verbose than it needed to be. I condensed it considerably here, used built-in PHP functions for the radians/degrees conversion, and put the whole thing into a class:

    <?php
    class PointRadiusCircle implements \Iterator
    {
        /** @see https://en.wikipedia.org/wiki/Earth_ellipsoid */
        public const float IERS_EQ_RADIUS = 6378136.6;
    
        private array $coordinates = [];
        private int $position = 0;
    
        /**
         * @param float $longitude
         * @param float $latitude
         * @param int $radius the circle's radius in meters
         * @param int $points number of points in the polygon
         */
        public function __construct(
            public float $longitude,
            public float $latitude,
            public int $radius,
            public int $points = 360
        )
        {
            for ($i = 0; $i < $points; $i++) {
                $bearing = 2 * M_PI * $i / $points;
                $this->coordinates[] = $this->getOffset($longitude, $latitude, $radius, $bearing);
            }
            $this->coordinates[] = $this->coordinates[0];
        }
    
        /**
         * @return array<float[]> An array of coordinates in longitude, latitude order
         */
        public function toArray(): array
        {
            return $this->coordinates;
        }
    
        /**
         * @param bool $pretty whether to format output for human readability
         * @return string GeoJSON Polygon string
         * @throws JsonException
         */
        public function toGeoJson(bool $pretty = false): string
        {
            return json_encode(
                ["type" => "Polygon", "coordinates" => [$this->coordinates]],
                \JSON_THROW_ON_ERROR | ($pretty ? \JSON_PRETTY_PRINT : 0)
            );
        }
    
        public function rewind(): void {
            $this->position = 0;
        }
    
        public function current(): array {
            return $this->coordinates[$this->position];
        }
    
        public function key(): int {
            return $this->position;
        }
    
        public function next(): void {
            ++$this->position;
        }
    
        public function valid(): bool {
            return isset($this->coordinates[$this->position]);
        }
    
        /**
         * @param float $longitude start longitude
         * @param float $latitude start latitude
         * @param float $distance distance in meters
         * @param float $bearing bearing in radians
         * @return float[] coordinates in longitude, latitude order
         */
        private function getOffset(float $longitude, float $latitude, float $distance, float $bearing): array {
            $lon1 = deg2rad($longitude);
            $lat1 = deg2rad($latitude);
            // convert dist to angular distance in radians
            $dByR = $distance / self::IERS_EQ_RADIUS;
    
            $lat = asin(
                sin($lat1) * cos($dByR) +
                cos($lat1) * sin($dByR) * cos($bearing)
            );
            $lon = $lon1 + atan2(
                    sin($bearing) * sin($dByR) * cos($lat1),
                    cos($dByR) - sin($lat1) * sin($lat)
                );
            $lon = fmod(
                    $lon + 3 * M_PI,
                    2 * M_PI
                ) - M_PI;
    
            return [rad2deg($lon), rad2deg($lat)];
        }
    }
    
    $result = new \PointRadiusCircle(-122, 49, 5000, 20);
    foreach ($result as $v) {
        vprintf("%.8f\t%.8f\n", $v);
    }
    // echo $result->toGeoJson();
    

    Output:

    -122.00000000   49.04491576
    -121.97882561   49.04271549
    -121.95972908   49.03633061
    -121.94458291   49.02638756
    -121.93486969   49.01386141
    -121.93153703   48.99997975
    -121.93490598   48.98610195
    -121.94464163   48.97358593
    -121.95978780   48.96365539
    -121.97886190   48.95728064
    -122.00000000   48.95508424
    -122.02113810   48.95728064
    -122.04021220   48.96365539
    -122.05535837   48.97358593
    -122.06509402   48.98610195
    -122.06846297   48.99997975
    -122.06513031   49.01386141
    -122.05541709   49.02638756
    -122.04027092   49.03633061
    -122.02117439   49.04271549
    -122.00000000   49.04491576
    

    As you can see in this test it draws the circle as expected.