From 6bb4f36df5f18d94306cd508dd65639d894a06e1 Mon Sep 17 00:00:00 2001 From: Andy Eschbacher Date: Wed, 30 Mar 2016 08:10:35 -0400 Subject: [PATCH] extracting util code to new submodule --- .../crankshaft/crankshaft/clustering/moran.py | 214 +++--------------- .../crankshaft/pysal_utils/__init__.py | 1 + .../crankshaft/pysal_utils/pysal_utils.py | 149 ++++++++++++ .../crankshaft/test/test_clustering_moran.py | 21 +- 4 files changed, 196 insertions(+), 189 deletions(-) create mode 100644 src/py/crankshaft/crankshaft/pysal_utils/__init__.py create mode 100644 src/py/crankshaft/crankshaft/pysal_utils/pysal_utils.py diff --git a/src/py/crankshaft/crankshaft/clustering/moran.py b/src/py/crankshaft/crankshaft/clustering/moran.py index 7f402a8..2a043c3 100644 --- a/src/py/crankshaft/crankshaft/clustering/moran.py +++ b/src/py/crankshaft/crankshaft/clustering/moran.py @@ -5,10 +5,12 @@ 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 numpy as np import pysal as ps import plpy +# crankshaft module +import crankshaft.pysal_utils as pu + # High level interface --------------------------------------- def moran(subquery, attr_name, @@ -25,7 +27,7 @@ def moran(subquery, attr_name, "subquery": subquery, "num_ngbrs": num_ngbrs} - query = construct_neighbor_query(w_type, qvals) + query = pu.construct_neighbor_query(w_type, qvals) plpy.notice('** Query: %s' % query) @@ -33,22 +35,23 @@ def moran(subquery, attr_name, result = plpy.execute(query) # if there are no neighbors, exit if len(result) == 0: - return empty_zipped_array(2) + return pu.empty_zipped_array(2) plpy.notice('** Query returned with %d rows' % len(result)) except plpy.SPIError: plpy.error('Error: areas of interest query failed, check input parameters') plpy.notice('** Query failed: "%s"' % query) plpy.notice('** Error: %s' % plpy.SPIError) - return empty_zipped_array(2) + return pu.empty_zipped_array(2) ## collect attributes - attr_vals = get_attributes(result) + attr_vals = pu.get_attributes(result) ## calculate weights - weight = get_weight(result, w_type, num_ngbrs) + weight = pu.get_weight(result, w_type, num_ngbrs) ## calculate moran global - moran_global = ps.esda.moran.Moran(attr_vals, weight, permutations=permutations) + moran_global = ps.esda.moran.Moran(attr_vals, weight, + permutations=permutations) return zip([moran_global.I], [moran_global.EI]) @@ -68,20 +71,20 @@ def moran_local(subquery, attr, "subquery": subquery, "num_ngbrs": num_ngbrs} - query = construct_neighbor_query(w_type, qvals) + query = pu.construct_neighbor_query(w_type, qvals) try: result = plpy.execute(query) # if there are no neighbors, exit if len(result) == 0: - return empty_zipped_array(5) + return pu.empty_zipped_array(5) except plpy.SPIError: plpy.error('Error: areas of interest query failed, check input parameters') plpy.notice('** Query failed: "%s"' % query) - return empty_zipped_array(5) + return pu.empty_zipped_array(5) - attr_vals = get_attributes(result) - weight = get_weight(result, w_type) + attr_vals = pu.get_attributes(result) + weight = pu.get_weight(result, w_type, num_ngbrs) # calculate LISA values lisa = ps.esda.moran.Moran_Local(attr_vals, weight, @@ -90,7 +93,6 @@ def moran_local(subquery, attr, # find quadrants for each geometry quads = quad_position(lisa.q) - plpy.notice('** Finished calculations') return zip(lisa.Is, quads, lisa.p_sim, weight.id_order, lisa.y) def moran_rate(subquery, numerator, denominator, @@ -106,7 +108,7 @@ def moran_rate(subquery, numerator, denominator, "subquery": subquery, "num_ngbrs": num_ngbrs} - query = construct_neighbor_query(w_type, qvals) + query = pu.construct_neighbor_query(w_type, qvals) plpy.notice('** Query: %s' % query) @@ -114,19 +116,19 @@ def moran_rate(subquery, numerator, denominator, result = plpy.execute(query) # if there are no neighbors, exit if len(result) == 0: - return empty_zipped_array(2) + return pu.empty_zipped_array(2) plpy.notice('** Query returned with %d rows' % len(result)) except plpy.SPIError: plpy.error('Error: areas of interest query failed, check input parameters') plpy.notice('** Query failed: "%s"' % query) plpy.notice('** Error: %s' % plpy.SPIError) - return empty_zipped_array(2) + return pu.empty_zipped_array(2) ## collect attributes - numer = get_attributes(result, 1) - denom = get_attributes(result, 2) + numer = pu.get_attributes(result, 1) + denom = pu.get_attributes(result, 2) - weight = get_weight(result, w_type, num_ngbrs) + weight = pu.get_weight(result, w_type, num_ngbrs) ## calculate moran global rate lisa_rate = ps.esda.moran.Moran_Rate(numer, denom, weight, @@ -143,7 +145,7 @@ def moran_local_rate(subquery, numerator, denominator, # geometries with values that are null are ignored # resulting in a collection of not as near neighbors - query = construct_neighbor_query(w_type, + query = pu.construct_neighbor_query(w_type, {"id_col": id_col, "numerator": numerator, "denominator": denominator, @@ -155,18 +157,18 @@ def moran_local_rate(subquery, numerator, denominator, result = plpy.execute(query) # if there are no neighbors, exit if len(result) == 0: - return empty_zipped_array(5) + return pu.empty_zipped_array(5) except plpy.SPIError: plpy.error('Error: areas of interest query failed, check input parameters') plpy.notice('** Query failed: "%s"' % query) plpy.notice('** Error: %s' % plpy.SPIError) - return empty_zipped_array(5) + return pu.empty_zipped_array(5) ## collect attributes - numer = get_attributes(result, 1) - denom = get_attributes(result, 2) + numer = pu.get_attributes(result, 1) + denom = pu.get_attributes(result, 2) - weight = get_weight(result, w_type, num_ngbrs) + weight = pu.get_weight(result, w_type, num_ngbrs) # calculate LISA values lisa = ps.esda.moran.Moran_Local_Rate(numer, denom, weight, @@ -191,25 +193,25 @@ def moran_local_bv(subquery, attr1, attr2, "geom_col": geom_col, "id_col": id_col} - query = construct_neighbor_query(w_type, qvals) + query = pu.construct_neighbor_query(w_type, qvals) try: result = plpy.execute(query) # if there are no neighbors, exit if len(result) == 0: - return empty_zipped_array(4) + return pu.empty_zipped_array(4) except plpy.SPIError: plpy.error("Error: areas of interest query failed, " \ "check input parameters") plpy.notice('** Query failed: "%s"' % query) - return empty_zipped_array(4) + return pu.empty_zipped_array(4) ## collect attributes - attr1_vals = get_attributes(result, 1) - attr2_vals = get_attributes(result, 2) + attr1_vals = pu.get_attributes(result, 1) + attr2_vals = pu.get_attributes(result, 2) # create weights - weight = get_weight(result, w_type, num_ngbrs) + weight = pu.get_weight(result, w_type, num_ngbrs) # calculate LISA values lisa = ps.esda.moran.Moran_Local_BV(attr1_vals, attr2_vals, weight, @@ -246,138 +248,6 @@ def map_quads(coord): else: return None -def query_attr_select(params): - """ - Create portion of SELECT statement for attributes inolved in query. - @param params: dict of information used in query (column names, - table name, etc.) - """ - - attrs = [k for k in params - if k not in ('id_col', 'geom_col', 'subquery', 'num_ngbrs')] - - template = "i.\"{%(col)s}\"::numeric As attr%(alias_num)s, " - - attr_string = "" - - for idx, val in enumerate(sorted(attrs)): - attr_string += template % {"col": val, "alias_num": idx + 1} - - return attr_string - -def query_attr_where(params): - """ - Create portion of WHERE clauses for weeding out NULL-valued geometries - """ - attrs = sorted([k for k in params - if k not in ('id_col', 'geom_col', 'subquery', 'num_ngbrs')]) - - attr_string = [] - - for attr in attrs: - attr_string.append("idx_replace.\"{%s}\" IS NOT NULL" % attr) - - if len(attrs) == 2: - attr_string.append("idx_replace.\"{%s}\" <> 0" % attrs[1]) - - out = " AND ".join(attr_string) - - return out - -def knn(params): - """SQL query for k-nearest neighbors. - @param vars: dict of values to fill template - """ - - attr_select = query_attr_select(params) - attr_where = query_attr_where(params) - - replacements = {"attr_select": attr_select, - "attr_where_i": attr_where.replace("idx_replace", "i"), - "attr_where_j": attr_where.replace("idx_replace", "j")} - - query = "SELECT " \ - "i.\"{id_col}\" As id, " \ - "%(attr_select)s" \ - "(SELECT ARRAY(SELECT j.\"{id_col}\" " \ - "FROM ({subquery}) As j " \ - "WHERE %(attr_where_j)s AND " \ - "i.\"{id_col}\" <> j.\"{id_col}\" " \ - "ORDER BY j.\"{geom_col}\" <-> i.\"{geom_col}\" ASC " \ - "LIMIT {num_ngbrs}) " \ - ") As neighbors " \ - "FROM ({subquery}) As i " \ - "WHERE " \ - "%(attr_where_i)s " \ - "ORDER BY i.\"{id_col}\" ASC;" % replacements - - return query.format(**params) - -## SQL query for finding queens neighbors (all contiguous polygons) -def queen(params): - """SQL query for queen neighbors. - @param params dict: information to fill query - """ - attr_select = query_attr_select(params) - attr_where = query_attr_where(params) - - replacements = {"attr_select": attr_select, - "attr_where_i": attr_where.replace("idx_replace", "i"), - "attr_where_j": attr_where.replace("idx_replace", "j")} - - query = "SELECT " \ - "i.\"{id_col}\" As id, " \ - "%(attr_select)s" \ - "(SELECT ARRAY(SELECT j.\"{id_col}\" " \ - "FROM ({subquery}) As j " \ - "WHERE ST_Touches(i.\"{geom_col}\", j.\"{geom_col}\") AND " \ - "%(attr_where_j)s)" \ - ") As neighbors " \ - "FROM ({subquery}) As i " \ - "WHERE " \ - "%(attr_where_i)s " \ - "ORDER BY i.\"{id_col}\" ASC;" % replacements - - return query.format(**params) - -## to add more weight methods open a ticket or pull request - -def construct_neighbor_query(w_type, query_vals): - """Return requested query. - @param w_type text: type of neighbors to calculate ('knn' or 'queen') - @param query_vals dict: values used to construct the query - """ - - if w_type == 'knn': - return knn(query_vals) - else: - return queen(query_vals) - -def get_attributes(query_res, attr_num=1): - """ - @param query_res: query results with attributes and neighbors - @param attr_num: attribute number (1, 2, ...) - """ - return np.array([x['attr' + str(attr_num)] for x in query_res], dtype=np.float) - -## Build weight object -def get_weight(query_res, w_type='queen', num_ngbrs=5): - """ - Construct PySAL weight from return value of query - @param query_res: query results with attributes and neighbors - """ - if w_type == 'knn': - row_normed_weights = [1.0 / float(num_ngbrs)] * num_ngbrs - weights = {x['id']: row_normed_weights for x in query_res} - else: - weights = {x['id']: [1.0 / len(x['neighbors'])] * len(x['neighbors']) - if len(x['neighbors']) > 0 - else [] for x in query_res} - - neighbors = {x['id']: x['neighbors'] for x in query_res} - - return ps.W(neighbors, weights) - def quad_position(quads): """ Produce Moran's I classification based of n @@ -385,20 +255,6 @@ def quad_position(quads): @param quads ndarray: an array of quads classified by 1-4 (PySAL default) Output: - @param ndarray: an array of quads classied by 'HH', 'LL', etc. + @param list: an array of quads classied by 'HH', 'LL', etc. """ - - lisa_sig = np.array([map_quads(q) for q in quads]) - - return lisa_sig - -def return_empty_zipped_array(num_nones): - """ - prepare return values for cases of empty weights objects (no neighbors) - Input: - @param num_nones int: number of columns (e.g., 4) - Output: - [(None, None, None, None)] - """ - - return [tuple([None] * num_nones)] + return [map_quads(q) for q in quads] diff --git a/src/py/crankshaft/crankshaft/pysal_utils/__init__.py b/src/py/crankshaft/crankshaft/pysal_utils/__init__.py new file mode 100644 index 0000000..835880d --- /dev/null +++ b/src/py/crankshaft/crankshaft/pysal_utils/__init__.py @@ -0,0 +1 @@ +from pysal_utils import * diff --git a/src/py/crankshaft/crankshaft/pysal_utils/pysal_utils.py b/src/py/crankshaft/crankshaft/pysal_utils/pysal_utils.py new file mode 100644 index 0000000..5482cc7 --- /dev/null +++ b/src/py/crankshaft/crankshaft/pysal_utils/pysal_utils.py @@ -0,0 +1,149 @@ +""" + Utilities module for generic PySAL functionality, mainly centered on translating queries into numpy arrays or PySAL weights objects +""" + +import numpy as np +import pysal as ps + +def construct_neighbor_query(w_type, query_vals): + """Return query (a string) used for finding neighbors + @param w_type text: type of neighbors to calculate ('knn' or 'queen') + @param query_vals dict: values used to construct the query + """ + + if w_type == 'knn': + return knn(query_vals) + else: + return queen(query_vals) + +## Build weight object +def get_weight(query_res, w_type='knn', num_ngbrs=5): + """ + Construct PySAL weight from return value of query + @param query_res: query results with attributes and neighbors + """ + if w_type == 'knn': + row_normed_weights = [1.0 / float(num_ngbrs)] * num_ngbrs + weights = {x['id']: row_normed_weights for x in query_res} + else: + weights = {x['id']: [1.0 / len(x['neighbors'])] * len(x['neighbors']) + if len(x['neighbors']) > 0 + else [] for x in query_res} + + neighbors = {x['id']: x['neighbors'] for x in query_res} + + return ps.W(neighbors, weights) + +def query_attr_select(params): + """ + Create portion of SELECT statement for attributes inolved in query. + @param params: dict of information used in query (column names, + table name, etc.) + """ + + attrs = [k for k in params + if k not in ('id_col', 'geom_col', 'subquery', 'num_ngbrs')] + + template = "i.\"{%(col)s}\"::numeric As attr%(alias_num)s, " + + attr_string = "" + + for idx, val in enumerate(sorted(attrs)): + attr_string += template % {"col": val, "alias_num": idx + 1} + + return attr_string + +def query_attr_where(params): + """ + Create portion of WHERE clauses for weeding out NULL-valued geometries + """ + attrs = sorted([k for k in params + if k not in ('id_col', 'geom_col', 'subquery', 'num_ngbrs')]) + + attr_string = [] + + for attr in attrs: + attr_string.append("idx_replace.\"{%s}\" IS NOT NULL" % attr) + + if len(attrs) == 2: + attr_string.append("idx_replace.\"{%s}\" <> 0" % attrs[1]) + + out = " AND ".join(attr_string) + + return out + +def knn(params): + """SQL query for k-nearest neighbors. + @param vars: dict of values to fill template + """ + + attr_select = query_attr_select(params) + attr_where = query_attr_where(params) + + replacements = {"attr_select": attr_select, + "attr_where_i": attr_where.replace("idx_replace", "i"), + "attr_where_j": attr_where.replace("idx_replace", "j")} + + query = "SELECT " \ + "i.\"{id_col}\" As id, " \ + "%(attr_select)s" \ + "(SELECT ARRAY(SELECT j.\"{id_col}\" " \ + "FROM ({subquery}) As j " \ + "WHERE %(attr_where_j)s AND " \ + "i.\"{id_col}\" <> j.\"{id_col}\" " \ + "ORDER BY j.\"{geom_col}\" <-> i.\"{geom_col}\" ASC " \ + "LIMIT {num_ngbrs}) " \ + ") As neighbors " \ + "FROM ({subquery}) As i " \ + "WHERE " \ + "%(attr_where_i)s " \ + "ORDER BY i.\"{id_col}\" ASC;" % replacements + + return query.format(**params) + +## SQL query for finding queens neighbors (all contiguous polygons) +def queen(params): + """SQL query for queen neighbors. + @param params dict: information to fill query + """ + attr_select = query_attr_select(params) + attr_where = query_attr_where(params) + + replacements = {"attr_select": attr_select, + "attr_where_i": attr_where.replace("idx_replace", "i"), + "attr_where_j": attr_where.replace("idx_replace", "j")} + + query = "SELECT " \ + "i.\"{id_col}\" As id, " \ + "%(attr_select)s" \ + "(SELECT ARRAY(SELECT j.\"{id_col}\" " \ + "FROM ({subquery}) As j " \ + "WHERE ST_Touches(i.\"{geom_col}\", j.\"{geom_col}\") AND " \ + "%(attr_where_j)s)" \ + ") As neighbors " \ + "FROM ({subquery}) As i " \ + "WHERE " \ + "%(attr_where_i)s " \ + "ORDER BY i.\"{id_col}\" ASC;" % replacements + + return query.format(**params) + +## to add more weight methods open a ticket or pull request + +def get_attributes(query_res, attr_num=1): + """ + @param query_res: query results with attributes and neighbors + @param attr_num: attribute number (1, 2, ...) + """ + return np.array([x['attr' + str(attr_num)] for x in query_res], dtype=np.float) + +def empty_zipped_array(num_nones): + """ + prepare return values for cases of empty weights objects (no neighbors) + Input: + @param num_nones int: number of columns (e.g., 4) + Output: + [(None, None, None, None)] + """ + + return [tuple([None] * num_nones)] diff --git a/src/py/crankshaft/test/test_clustering_moran.py b/src/py/crankshaft/test/test_clustering_moran.py index e2d2a50..399e0cb 100644 --- a/src/py/crankshaft/test/test_clustering_moran.py +++ b/src/py/crankshaft/test/test_clustering_moran.py @@ -12,6 +12,7 @@ import unittest from helper import plpy, fixture_file import crankshaft.clustering as cc +import crankshaft.pysal_utils as pu from crankshaft import random_seeds import json @@ -44,16 +45,16 @@ class MoranTest(unittest.TestCase): ans = "i.\"{attr1}\"::numeric As attr1, " \ "i.\"{attr2}\"::numeric As attr2, " - self.assertEqual(cc.query_attr_select(self.params), ans) + self.assertEqual(pu.query_attr_select(self.params), ans) def test_query_attr_where(self): - """Test query_attr_where""" + """Test pu.query_attr_where""" ans = "idx_replace.\"{attr1}\" IS NOT NULL AND " \ "idx_replace.\"{attr2}\" IS NOT NULL AND " \ "idx_replace.\"{attr2}\" <> 0" - self.assertEqual(cc.query_attr_where(self.params), ans) + self.assertEqual(pu.query_attr_where(self.params), ans) def test_knn(self): """Test knn neighbors constructor""" @@ -76,7 +77,7 @@ class MoranTest(unittest.TestCase): "i.\"jay_z\" <> 0 " \ "ORDER BY i.\"cartodb_id\" ASC;" - self.assertEqual(cc.knn(self.params), ans) + self.assertEqual(pu.knn(self.params), ans) def test_queen(self): """Test queen neighbors constructor""" @@ -90,7 +91,7 @@ class MoranTest(unittest.TestCase): "j.\"the_geom\") AND " \ "j.\"andy\" IS NOT NULL AND " \ "j.\"jay_z\" IS NOT NULL AND " \ - "j.\"jay_z\" <> 0) + "j.\"jay_z\" <> 0)" \ ") As neighbors " \ "FROM (SELECT * FROM a_list) As i " \ "WHERE i.\"andy\" IS NOT NULL AND " \ @@ -98,14 +99,14 @@ class MoranTest(unittest.TestCase): "i.\"jay_z\" <> 0 " \ "ORDER BY i.\"cartodb_id\" ASC;" - self.assertEqual(cc.queen(self.params), ans) + self.assertEqual(pu.queen(self.params), ans) def test_construct_neighbor_query(self): """Test construct_neighbor_query""" # Compare to raw knn query - self.assertEqual(cc.construct_neighbor_query('knn', self.params), - cc.knn(self.params)) + self.assertEqual(pu.construct_neighbor_query('knn', self.params), + pu.knn(self.params)) def test_get_attributes(self): """Test get_attributes""" @@ -123,8 +124,8 @@ class MoranTest(unittest.TestCase): """Test empty_zipped_array""" ans2 = [(None, None)] ans4 = [(None, None, None, None)] - self.assertEqual(cc.empty_zipped_array(2), ans2) - self.assertEqual(cc.empty_zipped_array(4), ans4) + self.assertEqual(pu.empty_zipped_array(2), ans2) + self.assertEqual(pu.empty_zipped_array(4), ans4) def test_quad_position(self): """Test lisa_sig_vals"""