How to make a circle geometry in KML or WKT

There is no circle geometry in KML.

In remote sensing, some algorithms e.g. MODIS wild fire detection, etc. produce sparse point data. But, I find it more meaningful to visualise these point data on a map, Google Earth layer for instance, as circles whose radii are based on the satellite imagery resolution.

As suggested in the KML FAQ, I approximate the circle geometry using regular n-sided polygon. The Python code below is originates from las3rjock’s posting on StackOverflow. I have adapted it slightly to accept input in metres, and to produce KML fragment and WKT. It is good enough when radius is small compared to the radius of the Earth.

Alternatively, get KML circle which can also do stars, or use the online generator.

from math import asin,cos,pi,sin
import re
class GeomHelper:
    rEarth_km = 6371.01 # Earth's average radius in km
    epsilon = 0.000001 # threshold for floating-point equality
    @staticmethod
    def deg2rad(angle):
        return angle*pi/180
    @staticmethod
    def rad2deg(angle):
        return angle*180/pi
    @staticmethod
    def pointRadialDistance(lat1_deg, lon1_deg, bearing_deg, distance_m):
        """
        Return final coordinates (lat2,lon2) [in degrees] given initial coordinates
        (lat1,lon1) [in degrees] and a bearing [in degrees] and distance [in m]
        """
        rlat1 = GeomHelper.deg2rad(lat1_deg)
        rlon1 = GeomHelper.deg2rad(lon1_deg)
        rbearing = GeomHelper.deg2rad(bearing_deg)
        rdistance = distance_m / (GeomHelper.rEarth_km*1000) # normalize linear distance to radian angle
        rlat = asin( sin(rlat1) * cos(rdistance) + cos(rlat1) * sin(rdistance) * cos(rbearing) )
        if cos(rlat) == 0 or abs(cos(rlat)) < GeomHelper.epsilon: # Endpoint a pole
            rlon=rlon1
        else:
            rlon = ( (rlon1 - asin( sin(rbearing)* sin(rdistance) / cos(rlat) ) + pi ) % (2*pi) ) - pi
        lat = GeomHelper.rad2deg(rlat)
        lon = GeomHelper.rad2deg(rlon)
        return (lat, lon)

    @staticmethod
    def make_polygon(lon_deg, lat_deg, radius_m, sides=8, format='wkt'):
        theta = 360/sides
        if format is 'wkt':
            separator = ' '
        else:
            separator = ','
        result = []
        for i in range(sides):
            lat, lon = GeomHelper.pointRadialDistance(lat_deg, lon_deg, i*theta, radius_m)
            result.append("%f%s%f" % (lon, separator, lat))
        # Close the polygon
        result.append(result[0])
        if format is 'wkt':
            return 'POLYGON((%s))' % ','.join(result)
        else:
            # KML
            return '%s' % ' '.join(result)
    @staticmethod
    def make_polygon_from_wkt(point, radius_m, **kwargs):
        """ 'point' is a string like POINT(30,-25.45) """
        m = re.search('^point\((\S+) ([^)]+)\)$', point, re.I)
        if m:
            return GeomHelper.make_polygon(float(m.group(1)), float(m.group(2)), radius_m, **kwargs)
        else:
            raise ValueError("Invalid WKT POINT: %s" % point)
    @staticmethod
    def make_polygon_from_kml(point, radius_m, **kwargs):
        """ 'point' is a string like <POINT><coordinates>27.113,-25.389</coordinates></POINT> """
        m = re.search('<point><coordinates>([^,]+),([^<]+)', point, re.I)
        if m:
            return GeomHelper.make_polygon(float(m.group(1)), float(m.group(2)), radius_m, **kwargs)
        else:
            raise ValueError("Invalid KML POINT: %s" % point)

Advertisement