Standard Deviational Ellipse with Open Source Python (GDAL/OGR etc)

Question

I am trying to mimic the Standard Deviational Ellipse tool from ArcGIS through Open Source Python avenues. The formula is here.

At the moment I am failing to get the correct X and Y axis lengths for the ellipse, the first step in the formula.

When I run the tool in ArcGIS the values are 4607.796039, 7816.667667 but I get 3948.758414, 5057.00957766.

Here's my code so far...

from osgeo import ogr
from shapely.geometry import MultiLineString
from shapely import wkt
import numpy as np
import sys, math

## set the driver for the data
driver = ogr.GetDriverByName("FileGDB")
## path to the FileGDB
gdb = r"C:\Users\******\Documents\ArcGIS\Default.gdb"
## ope the GDB in write mode (1)
ds = driver.Open(gdb, 1)

input_lyr_name = "Birmingham_Burglaries_2016"

output_fc = input_lyr_name + "_standard_ellipse"

## reference the layer using the layers name
if input_lyr_name in [ds.GetLayerByIndex(lyr_name).GetName() for lyr_name in range(ds.GetLayerCount())]:
    lyr = ds.GetLayerByName(input_lyr_name)
    print "{0} found in {1}".format(input_lyr_name, gdb)

if output_fc in [ds.GetLayerByIndex(lyr_name).GetName() for lyr_name in range(ds.GetLayerCount())]:
    ds.DeleteLayer(output_fc)
    print "Deleting: {0}".format(output_fc)

## for each point in the layer
## get the x and y value
## and place in an array
try:
    first_feat = lyr.GetFeature(1)
    if first_feat.geometry().GetGeometryName() in ["POINT", "MULTIPOINT", "POLYGON", "MULTIPOLYGON"]:
        xy_arr = np.ndarray((len(lyr), 2), dtype=np.float)
        for i, pt in enumerate(lyr):
            ft_geom = pt.geometry()
            xy_arr[i] = (ft_geom.Centroid().GetX(), ft_geom.Centroid().GetY())

    ## for lines we get the midpoint of a line
    elif first_feat.geometry().GetGeometryName() in ["LINESTRING", "MULTILINESTRING"]:
        xy_arr = np.ndarray((len(lyr), 2), dtype=np.float)
        for i, ln in enumerate(lyr):
            line_geom = ln.geometry().ExportToWkt()
            shapely_line = MultiLineString(wkt.loads(line_geom))
            midpoint = shapely_line.interpolate(shapely_line.length/2)
            xy_arr[i] = (midpoint.x, midpoint.y)

except Exception:
    print "Unknown geometry for {}".format(input_lyr_name)
    sys.exit()

## the mean center (average x and average y coordinate)
avg_x, avg_y = np.mean(xy_arr, axis=0)

print "Mean Center: {0}, {1}".format(avg_x, avg_y)

sum_of_sq_diff_x = 0.0
sum_of_sq_diff_y = 0.0

for x, y in xy_arr:
    # (x - xmean)squared 
    diff_x = math.pow(x - avg_x, 2)
    # (y - ymean)squared
    diff_y = math.pow(y - avg_y, 2)
    # sum the differences sqaured from above
    sum_of_sq_diff_x += diff_x
    sum_of_sq_diff_y += diff_y

# x axis length
sum_of_results_x = (sum_of_sq_diff_x/lyr.GetFeatureCount())
standard_distance_x = math.sqrt(sum_of_results_x)
# y axis length
sum_of_results_y = (sum_of_sq_diff_y/lyr.GetFeatureCount())
standard_distance_y = math.sqrt(sum_of_results_y)

print standard_distance_x, standard_distance_y

Answer 1

The Standard Deviational Ellipse from ArcGIS is not the only solution. In the geospatial world, there are two major algorithms (Yuill and CrimeStat III) and many intermediate solutions (QGIS: des Ellipses de Déviation Standard (SDE), un plugin, « Standard Deviational Ellipse », des scripts R (processing) et Python et une approche critique..., in French)

I am failing to get the correct X and Y axis lengths for the ellipse

Using the different solutions, the resulting values are all different (but scientifically correct) and you can use the QGIS Plugin: Standard Deviational Ellipse to compute all.

1) The one proposed by Yuill, R. S.(1971) in "The Standard Deviational Ellipse: An Updated Tool for Spatial Description", Geografiska Annaler 53B(1),28-39)

Equation

Before the 10.x version, ArcGIS uses this Yuill algorithm (ArcGIS 9.3: How Directional Distribution: Standard Deviational Ellipse (Spatial Statistics) works)

2) the one proposed by Ned Levine in CrimeStat III (Ned Levine, 2010, A Spatial Statistics Program for the Analysis of Crime Incident Locations (version 3.3). Ned Levine & Associates, Houston, TX.; National Institute of Justice, Washington, DC)

Equation

In addition to the plugin, you can also use the R aspace package (calc_sde) and the aspace_SDE.rsx QGIS processing script.

3) Between the version 9.x and 10.x, ArcGIS has changed its algorithm (ArcGIS 10.3: How Directional Distribution (Standard Deviational Ellipse) works))

4) a new algorithm proposed by Wang et al. (2015) , Confidence Analysis of Standard Deviational Ellipse and Its Extension into Higher Dimensional Euclidean Space (ellipse_wang.py (in french)

Now, what is the best ?

For those who are only interested in the angles and eccentricity, no problem