Standard Deviational Ellipse differs from QGIS, sfdep

[ljwolf]

Hi,

I've gotten a bugreport from Nick Tate by email:

description

I spent some time Friday/Monday looking at the sfdep::st_ellipse() and st_dev_ellipse() functions. I think the rotation angle needs to be multiplied by -1 for the results to make sense. BUT the parameters of the ellipse are not quite the same as the equivalent functions in PySAL pointpats: in the single example I ran, it was the major axis that was quite different whereas the other parameters were the same. And both seem to be based on Ned Levine's CrimeStat algorithm (I think).

I tried an independent (I think) Python implementation in QGIS which agreed with the R function parameters. Of course, it may be that the pointpats algorithm is slightly different...

I have asked for an MWE. At the same time, I've taken a quick look at the code.

sfdep flips the angle to the positive direction.

However, we use numpy.arctan(), which is only defined in [-90º,90º]. This means it will treat 135º as if it is -45º. This is important, since some algorithms assume 0 < theta < 2pi while others assume |theta| < pi/2, and you have to identify this when the authors define their trigonometry functions. I remember this rough edge when testing esda.shape before publishing it in my PhD.

# these will be the same since arctan assumes the angle is in [-90º,90º]
#                     y    x
theta1 = numpy.arctan(-1 / 1)
theta2 = numpy.arctan(1 / -1)
# these will differ, since arctan2 chooses the quadrant correctly based on the signs of x,y
#                      y   x
theta3 = numpy.arctan2(-1, 1)
theta4 = numpy.arctan2(1, -1)

print(f"""
The angle formed by 1/-1 should differ from the angle formed by 1/-1
arctan: {numpy.rad2deg(theta1), numpy.rad2deg(theta3)}
arctan2: {numpy.rad2deg(theta2), numpy.rad2deg(theta4)}
""")

yields

# The angle formed by 1/-1 should differ from the angle formed by 1/-1
# arctan: (-45.0, -45.0)
# arctan2: (-45.0, 135.0)

I think this is likely the source of the issue given what Nick says & sfdep's implementation. numpy.cos(numpy.arctan(y/x)) is always positive, while numpy.cos(numpy.arctan2(y,x))'s sign will be the same as x.