1

I use this python code to calculate a greatcircle: https://code.google.com/p/pyroms/source/browse/trunk/pyroms/greatcircle.py?r=39 After initializing the greatcircle I can calculate points on the line.

I would like to do the same with an arbitrary smallcircle. I would like to define the smallcircle using its center location (in latlon) and one location on the line (also latlon), so not using a diameter. Say as an input I would use a certain point spacing in meters, how would I then proceed to calculate the other locations on the smallcircle line?

Thanks!

1 Answer 1

1

Python:

# From: http://www.movable-type.co.uk/scripts/latlong.html

def distance(lon1, lat1, lon2, lat2):
    # http://code.activestate.com/recipes/576779-calculating-distance-between-two-geographic-points/
    # http://en.wikipedia.org/wiki/Haversine_formula
    dlat = lat2 - lat1
    dlon = lon2 - lon1
    q = sin(dlat/2)**2 + (cos(lat1) * cos(lat2) * (sin(dlon/2)**2))
    return 2 * atan2(sqrt(q), sqrt(1-q))

from math import *

def SmallCircle(radius,clon,clat,llon,llat,meters):
    # clon = circle central longitude (deg)
    # clat = circle central latitude (deg)
    # llon = longitude of point on circle line (deg)
    # llat = latitude of point on circle line (deg)
    # meters = point spacing (meters). This is not the exact spacing, it will be used to calculate nrpoints.
    # radius = radius of sphere

    dlat = radians(llat-clat)
    dlon = radians(llon-clon)
    clat = radians(clat)
    clon = radians(clon)
    llat = radians(llat)
    llon = radians(llon)

    a = sin(dlat/2) * sin(dlat/2) + sin(dlon/2) * sin(dlon/2) * cos(clat) * cos(llat)
    c = 2 * atan2(sqrt(a), sqrt(1-a))
    d = radius * c

    # bearing
    y = sin(dlon) * cos(llat)
    x = cos(clat)*sin(llat) - sin(clat)*cos(llat)*cos(dlon)
    brng = atan2(y, x)

    # circle
    angle = distance(clon,clat,llon,llat)
    circleradius = sin(angle) * radius
    circlelength = 2*pi*circleradius

    # loop circle
    list = []
    list.append([degrees(llon),degrees(llat)])
    nrpoints = int(round(circlelength / meters))
    for i in range(1, nrpoints):
        brng = brng + (float(i) / float(nrpoints))*2*pi
        lat = asin( sin(clat)*cos(d/radius) + cos(clat)*sin(d/radius)*cos(brng) )
        lon = clon + atan2(sin(brng)*sin(d/radius)*cos(clat), cos(d/radius)-sin(clat)*sin(lat))
        list.append([degrees(lon),degrees(lat)])

    return list

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.