fit_shadow.py

import numpy as np
import shapefile
import sys
from scipy.spatial import cKDTree
import scipy.optimize


def ellipse(x, theta):
    phi = x[4]
    r = 1 / np.sqrt((np.cos(theta)) ** 2 + (np.sin(theta)) ** 2)
    xx = r * np.cos(theta)
    yy = r * np.sin(theta)
    data = np.array([xx, yy])
    S = np.array([[x[2], 0], [0, x[3]]])
    R = np.array([[np.cos(phi), -np.sin(phi)], [np.sin(phi), np.cos(phi)]])
    T = np.dot(S, data)
    data = np.dot(R, T)
    data[0] += x[0]
    data[1] += x[1]
    return data.T


def objfn(x, args):
    theta, xy = args[0], args[1]
    Xd = ellipse(x, theta)
    tree = cKDTree(Xd)
    r, ii = tree.query(list(xy))
    return r


def fit(xs, theta):
    xc = xs.mean(0)
    x0 = np.array([xc[0], xc[1], 2, 1, 0])  # intial conditions for fit, [lon, lat, major, minor, rotation]
    x, success = scipy.optimize.leastsq(objfn, x0, args=[theta, xs], ftol=1e-9, xtol=1e-9, maxfev=10000)
    return x


# Download the low res shape file from: https://svs.gsfc.nasa.gov/4518
path = 'path_data/eclipse2017/'
shp = open(path + 'w_umbra17.shp')
dbf = open(path + 'w_umbra17.dbf')
sf = shapefile.Reader(shp=shp, dbf=dbf)

theta = np.linspace(0, 2 * np.pi, 360)  # point to generate for the theoretical ellipse

# For pylab plotting in pylab. run %pylab in the ipython.
if len(sys.argv) > 1:
    shape_record = sf.shapeRecords()[int(sys.argv[1])]
    xs = np.array([[xy[0], xy[1]] for xy in shape_record.shape.points])  # shadow outline xy coordinates
    x = fit(xs, theta)
    Xd = ellipse(x, theta)

    print('Coords: %.2f %.2f' % (x[0], x[1]))
    print('SemiAxis: %.2f %.2f' % (x[2], x[3]))
    print('Rotation: %.2f' % (x[4]))

    clf()
    plot(xs[:, 0], xs[:, 1], '-g')
    plot(Xd[:, 0], Xd[:, 1], '-b')
    draw()

else:
    shadows = {}
    for i in range(20):
        shape_record = sf.shapeRecords()[i]
        xs = np.array([[xy[0], xy[1]] for xy in shape_record.shape.points])  # shadow outline xy coordinates
        print(shape_record.record)
        x = fit(xs, theta)
        shadows[shape_record.record[1]] = [s for s in x]
        print('Coords: %.2f %.2f' % (x[0], x[1]))
        print('SemiAxis: %.2f %.2f' % (x[2], x[3]))
        print('Rotation: %.2f' % (x[4]))