diff --git a/src/pg/sql/16_getis.sql b/src/pg/sql/16_getis.sql new file mode 100644 index 0000000..37f288b --- /dev/null +++ b/src/pg/sql/16_getis.sql @@ -0,0 +1,16 @@ +-- Getis-Ord's G +-- Hotspot/Coldspot Analysis tool +CREATE OR REPLACE FUNCTION + CDB_GetisOrdsG( + subquery TEXT, + column_name TEXT, + w_type TEXT, + num_ngbrs INT, + permutations INT, + geom_col TEXT, + id_col TEXT) +RETURNS TABLE (z_val NUMERIC, p_val NUMERIC, rowid BIGINT) +AS $$ + from crankshaft.clustering import getis + return getis_ord(subquery, column_name, w_type, num_ngbrs, permutations, geom_col, id_col) +$$ LANGUAGE plpythonu; diff --git a/src/py/crankshaft/crankshaft/clustering/__init__.py b/src/py/crankshaft/crankshaft/clustering/__init__.py index ed34fe0..d9682fa 100644 --- a/src/py/crankshaft/crankshaft/clustering/__init__.py +++ b/src/py/crankshaft/crankshaft/clustering/__init__.py @@ -1,3 +1,4 @@ """Import all functions from for clustering""" from moran import * from kmeans import * +from getis import * diff --git a/src/py/crankshaft/crankshaft/clustering/getis.py b/src/py/crankshaft/crankshaft/clustering/getis.py new file mode 100644 index 0000000..29a2c50 --- /dev/null +++ b/src/py/crankshaft/crankshaft/clustering/getis.py @@ -0,0 +1,51 @@ +""" +Moran's I geostatistics (global clustering & outliers presence) +""" + +# TODO: Fill in local neighbors which have null/NoneType values with the +# average of the their neighborhood + +import pysal as ps +import plpy +from collections import OrderedDict + +# crankshaft module +import crankshaft.pysal_utils as pu + +# High level interface --------------------------------------- + +def getis_ord(subquery, attr, + w_type, num_ngbrs, permutations, geom_col, id_col): + """ + Getis-Ord's G + Implementation building neighbors with a PostGIS database and Getis-Ord's G + hotspot/coldspot analysis with PySAL. + Andy Eschbacher + """ + + # geometries with attributes that are null are ignored + # resulting in a collection of not as near neighbors if kNN is chosen + + qvals = OrderedDict([("id_col", id_col), + ("attr1", attr), + ("geom_col", geom_col), + ("subquery", subquery), + ("num_ngbrs", num_ngbrs)]) + + query = pu.construct_neighbor_query(w_type, qvals) + + try: + result = plpy.execute(query) + # if there are no neighbors, exit + if len(result) == 0: + return pu.empty_zipped_array(3) + except plpy.SPIError, e: + plpy.error('Query failed: %s' % e) + + attr_vals = pu.get_attributes(result) + weight = pu.get_weight(result, w_type, num_ngbrs) + + # calculate LISA values + getis = ps.esda.getisord(attr_vals, weight, star=True) + + return zip(getis.z_sim, getis.p_sim, weight.id_order)