# -*- coding: utf-8 -*- """Clean polygons This module has some utils and functions to clean and transform (to line, point, polygon or polygon rectangularized) the output polygons from IA analysis. """ import warnings import numpy as np from django.conf import settings from sklearn.cluster import SpectralClustering from tools import layer from tools import const from MG.tools.postgis_api import PostGis # DEACTIVATES SCIKIT WARNINGS def warn(*_, **__): """ Function to deactivate scikit warnings """ return None warnings.warn = warn ########################### CHAIKIN_SMOOTHING = "ALTER TABLE \"{schema}\".\"{table}\" " \ "ALTER COLUMN {geom_field} TYPE geometry " \ "USING ST_ChaikinSmoothing({geom_field}, 1);" RENAME_TABLE = "ALTER TABLE \"{schema}\".\"{table_old}\" RENAME TO \"{table_new}\";" def clean_to_polygon(schema, table, geom_field): """ Function to clean and transform polygon into polygon. A chaikin smoothing is applied to the geometry. Parameters ---------- schema: str Schema where the table is saved table: str Name of the table geom_field: str Field name pointing to the geometry Returns ------- success: None """ postgis_obj = PostGis('Public') postgis_obj.send_sql_command(CHAIKIN_SMOOTHING.format(schema=schema, table=table, geom_field=geom_field)) GET_AXIS = "SELECT {function_schema}.medialaxis_con_limpieza_vw('{schema}', '{input_table}'," \ " '{output_table}', 'id,', '0.00005','0.00002')" def clean_to_line(schema, table, geom_field): """ Function to transform polygons into line. An Aproximate medial axis operation is applied on the geometries Parameters ---------- schema: str Schema where the table is saved table: str Name of the table geom_field: str Field name pointing to the geometry Returns ------- success: None """ postgis_obj = PostGis('Public') intermediate_link = postgis_obj.get_unique_name(schema=schema, prefix="line_clean") postgis_obj.send_sql_command(CHAIKIN_SMOOTHING.format(schema=schema, table=table, geom_field=geom_field)) postgis_obj.send_sql_command(RENAME_TABLE.format(schema=schema, table_old=table, table_new=intermediate_link)) postgis_obj.send_sql_command(GET_AXIS.format(function_schema=postgis_obj.schema_funciones, input_table=intermediate_link, output_table=table, schema=schema)) GET_CENTROIDS = "ALTER TABLE \"{schema}\".\"{table}\" ALTER COLUMN {geom_field} TYPE geometry " \ "USING ST_Centroid({geom_field});" def clean_to_point(schema, table, geom_field): """ Function to transform polygons into point. The centroid of the geometries is calculated Parameters ---------- schema: str Schema where the table is saved table: str Name of the table geom_field: str Field name pointing to the geometry Returns ------- success: None """ postgis_obj = PostGis('Public') postgis_obj.send_sql_command(GET_CENTROIDS.format(schema=schema, table=table, geom_field=geom_field)) RENAME_AND_UNION = "CREATE TABLE \"{schema}\".\"{new_name}\" AS SELECT " \ "(ST_DUMP(ST_UNION(ST_Buffer(ST_MakeValid(geom), 0.0)))).geom as geom" \ " from \"{schema}\".{old_name};" BACKUP_TABLE = "CREATE TABLE \"{schema}\".\"{table}_backup\" AS SELECT " \ "* from \"{schema}\".\"{table}\"" COPY_TABLE = "CREATE TABLE \"{schema}\".\"{new_table}\" AS SELECT" \ " entity, (ST_Dump(ST_Union(wkb_geometry))).geom as wkb_geometry FROM \"{schema}\".\"{old_table}\"" \ " group by entity;" MAKE_VALID_GEOM = "UPDATE \"{schema}\".\"{old_table}\" " \ "SET wkb_geometry = ST_MakeValid(wkb_geometry);" def clean(schema, table, geom_field, target_geometry, new_table_name=None): """ Master function to clean polygons. Transform the geometries of `table` to `target_geometry` type. Parameters ---------- schema: str Schema where the table is saved table: str Name of the table geom_field: str Field name pointing to the geometry target_geometry: str Target geometry which transform to Returns ------- success: None """ postgis_api = PostGis('Public') if target_geometry == 'polygon_r': clean_to_polygon(schema, table, geom_field) table_new = postgis_api.get_unique_name(schema) params = { 'user': schema, 'input_layer': layer.Layer(schema, { 'protocol': 'PostgreSQL', 'type': const.VECTOR_KEY, 'ip': postgis_api.ip, 'port': postgis_api.port, 'user': postgis_api.user, 'password': postgis_api.passw, 'database_name': postgis_api.dbname, 'schema': schema, 'table_view_name': table }), 'output_layer': layer.OutLayer(schema, { 'protocol': 'PostgreSQL', 'type': const.VECTOR_KEY, 'ip': postgis_api.ip, 'port': postgis_api.port, 'user': postgis_api.user, 'password': postgis_api.passw, 'database_name': postgis_api.dbname, 'schema': schema, 'table_view_name': table_new }) } settings.MASTER_TASK_LIST.run_into_task(schema, '/GMS/vectorial/rectangularizedPolygon', params) postgis_api.send_sql_command(BACKUP_TABLE.format(schema=schema, table=table)) try: postgis_api.delete_tables(schema, table, ) postgis_api.send_sql_command(RENAME_AND_UNION.format(schema=schema, old_name=table_new, new_name=table)) except Exception: print('Warning: polygon_r type fail, trying with polygon') postgis_api.send_sql_command(RENAME_TABLE.format(schema=schema, table_new=table, table_old="{}_backup".format(table))) clean_to_polygon(schema, table, geom_field) postgis_api.delete_tables(schema, "{}_backup".format(table), ) elif target_geometry == 'polygon': postgis_api.send_sql_command(MAKE_VALID_GEOM.format(schema=schema, old_table=table)) postgis_api.send_sql_command(COPY_TABLE.format(schema=schema, old_table=table, new_table=new_table_name)) clean_to_polygon(schema, new_table_name, geom_field) elif target_geometry == 'line': postgis_api.send_sql_command(MAKE_VALID_GEOM.format(schema=schema, old_table=table)) postgis_api.send_sql_command(COPY_TABLE.format(schema=schema, old_table=table, new_table=new_table_name)) clean_to_line(schema, new_table_name, geom_field) elif target_geometry == 'point': postgis_api.send_sql_command(MAKE_VALID_GEOM.format(schema=schema, old_table=table)) postgis_api.send_sql_command(COPY_TABLE.format(schema=schema, old_table=table, new_table=new_table_name)) clean_to_point(schema, new_table_name, geom_field) SQL_SELECT_ALL = "select * from \"{schema}\".\"{res_table}\" where id = '{ent_id}';" GET_AREA = "select ST_area(geom) from \"{schema}\".\"{table}\" WHERE {id_name}={id}" GET_ENTITIES_FROM_OTB = "CREATE TABLE \"{schema}\".\"{result_table}\" "\ " AS SELECT otb.*, "\ " (ST_area(st_intersection(ent.geom2,otb.geom))/" \ "(st_area(ent.geom2)+st_area(otb.geom)-ST_area(st_intersection(ent.geom2,otb.geom)))) AS IOU," \ " ent.id," \ " ST_area(st_intersection(ent.geom2,otb.geom)) as intersect_area," \ " st_area(otb.geom) as area " \ "FROM (SELECT *, wkb_geometry as geom FROM " \ "\"{schema}\".\"{otb_table}\" WHERE \"meanB0\">0 and \"meanB1\">0 and \"meanB2\">0 ) AS otb, "\ "(select geom as geom2,* from \"{schema}\".\"{ai_table}\") AS ent "\ " WHERE ST_INTERSECTS(otb.geom, ent.geom2) and " \ "(ST_area(st_intersection(ent.geom2,otb.geom))/" \ "(st_area(ent.geom2)+st_area(otb.geom)-ST_area(st_intersection(ent.geom2,otb.geom)))) >0.01; "\ "ALTER TABLE \"{schema}\".\"{result_table}\" drop column wkb_geometry; " GET_OGC_FID = "select ogc_fid from \"{schema}\".\"{table}\"" DELETE_ROW = "delete from \"{schema}\".\"{table}\" where ogc_fid = {ogc_fid};" COUNT_ENTITIES = "SELECT COUNT(*) from \"{schema}\".\"{table}\"" SECOND_CLEAN = "CREATE TABLE \"{schema}\".\"{table_new}\" AS SELECT " \ "ST_ConcaveHull((ST_Dump(ST_Simplify" \ "(ST_Buffer(ST_UNION(ST_Buffer(geom, 0.00002)),-0.00002), 0.000008, true))).geom,0.999, false) " \ "as geom from \"{schema}\".{table_old} GROUP BY id;" def _clean_otb_geometries(ent_id, schema, otb_table, ai_table): postgis_obj = PostGis('public') result = postgis_obj.send_sql_command(SQL_SELECT_ALL.format(schema=schema, res_table=otb_table, ent_id=ent_id)) result = list(result) if len(result) == 0: return [] elif len(result) == 1: return [result[0][1]] cluster_max = 15 # 8 data_trans = [] ids_array = [] areas_data = [] for i in result: data = [float(j) / 1000 for j in i[3:9]] data.append(float(i[10])) data_trans.append(data) ids_array.append(int(i[0])) areas_data.append(i[-2:]) data_trans = np.array(data_trans) if data_trans.shape[0] >= cluster_max + 1: num_cluster = cluster_max else: num_cluster = data_trans.shape[0] - 1 means = np.mean(data_trans, axis=0) means[-1] = 1 cluster_result = {} entity_area = float(list(postgis_obj.send_sql_command(GET_AREA.format(schema=schema, table=ai_table, id_name='id', id=ent_id)))[0][0]) iou = 0 final_entities = [] for n_cluster in range(1, num_cluster + 1): clustering = SpectralClustering(n_clusters=n_cluster, assign_labels="kmeans", random_state=0).fit( data_trans) min_res = 9999999999 for cluster_index in range(clustering.n_clusters): indexes = [index for index, label in enumerate(clustering.labels_) if label == cluster_index] data_selected = data_trans[indexes, :] data_selected_means = np.mean(data_selected, axis=0) diff = data_selected_means - means diff_norm = np.linalg.norm(diff) if diff_norm < min_res: min_res = diff_norm cluster_result[n_cluster] = indexes entities = cluster_result[n_cluster] intersect_area = 0 union_area_ent = 0 for entity in entities: intersect_area += areas_data[entity][0] union_area_ent += areas_data[entity][1] iou_now = intersect_area / (entity_area + union_area_ent - intersect_area) if iou_now > iou: final_entities = [ids_array[label_index] for label_index in cluster_result[n_cluster]] iou = iou_now return final_entities def clean_otb_geometries(schema, otb_table, ai_table, table_res): """ Function to clean and select the geometries from LSMS Parameters ---------- schema otb_table ai_table table_res Returns ------- """ postgis_obj = PostGis('public') intermediate_table = postgis_obj.get_unique_name(schema, 'intermediate') postgis_obj.send_sql_command(GET_ENTITIES_FROM_OTB.format(schema=schema, result_table=intermediate_table, otb_table=otb_table, ai_table=ai_table)) entities_not_to_remove = set() total_entities = int( list(postgis_obj.send_sql_command(COUNT_ENTITIES.format(schema=schema, table=ai_table)))[0][0]) for i in range(1, total_entities + 1): for j in _clean_otb_geometries(i, schema, intermediate_table, ai_table): entities_not_to_remove.add(j) cur = list(postgis_obj.send_sql_command(GET_OGC_FID.format(schema=schema, table=intermediate_table))) all_ids = [fid[0] for fid in cur] for ent in all_ids: if ent not in entities_not_to_remove: postgis_obj.send_sql_command(DELETE_ROW.format(schema=schema, table=intermediate_table, ogc_fid=ent)) postgis_obj.send_sql_command(SECOND_CLEAN.format(schema=schema, table_new=table_res, table_old=intermediate_table)) postgis_obj.delete_tables(schema, intermediate_table,) # postgis_obj.delete_tables(schema, ai_table,) # postgis_obj.change_name_table(schema, table_res, ai_table) return table_res