1

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
4

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

enter image description here

enter image description here

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

enter image description here

enter image description here

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))

enter image description here

enter image description here

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)

enter image description here

Now, what is the best ?

enter image description here

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

  • Very informative answer. The ellipse_wang.py was especially useful. I am currently trying to create my own tools and learn what goes on behind the scenes when I press the buttons. This gives me plenty to research and I have been able to get the ellipse into a polygon feature class. – Clubdebambos Apr 26 '17 at 14:42

Your Answer

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

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