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

Monday, 12 September 2016

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

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

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