1

I'm attempting to generate both a uniform distribution and a normal distribution of points that fall within an ellipse shape, made from a centre point, semi major and minor and rotation in degrees. I have code to do this (translated from MATLAB) and it plots as expected, but I was just looking for some assurance that it's doing what it is supposed to by others more knowledgeable than me. Would there be any need to account for projection issues or plotting this geographically that I'm not accounting for here?

normal uniform

from pyproj import transform, Proj
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns

def trans(lat, lon):
    return(transform(Proj(init='epsg:4326'), Proj(init='epsg:3857'), lon, lat))

def ellp(lat, lon, smaj, smin, rot):
    return(trans(lon, lat), smaj, smin, rot)

def dists(mean, sma, smi, a, d, norm=True):
    #https://www.quora.com/Is-it-possible-to-generate-random-points-within-a-given-rotated-ellipse-without-using-rejection-sampling
    theta = np.radians(a)

    if norm == True:
        rp = np.random.rand(1, d)
    else:
        rp = np.sqrt(np.random.rand(1, d))

    phi = np.random.rand(1,d)*2*np.pi

    x_c = rp*np.cos(phi)
    y_c = rp*np.sin(phi)

    x_e = x_c * sma
    y_e = y_c * smi

    x = np.array(x_e*np.cos(theta) - y_e*np.sin(theta)) + mean[1]
    y = np.array(x_e*np.sin(theta) + y_e*np.cos(theta)) + mean[0]

    return(y,x)

j = ellp(30, 40, 2000, 500 , 40)
f = dists(j[0], j[1], j[2], j[3], 10000, norm=False)
g = dists(j[0], j[1], j[2], j[3], 10000, norm=True)

plt.figure(figsize=(10,10))
ax = sns.scatterplot(x=f[0][0], y=f[1][0])
plt.axis('equal')
plt.show()

plt.figure(figsize=(10,10))
ax = sns.scatterplot(x=g[0][0], y=g[1][0])
plt.axis('equal')
plt.show()
  • 3
    The fact that you've used Web Mercator in your ellipse plotting algorithm likely dooms your effort, since it is seriously flawed for distance calculations -- If it looks right in 3857, it's probably wrong. – Vince Apr 14 at 13:23
  • Ah good to know, perhaps then a UTM conversion instead so that it can be more accurately translated would make a difference? – saph_top Apr 14 at 13:26
  • 1
    UTM is only valid in narrow longitude bands. An appropriate Equal-Area projection might be necessary. – Vince Apr 14 at 14:02

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.