Merge branch 'master' of github.com:swiss-territorial-data-lab/detect…

…or-interface into nh-hr-procedure

tools/prediction-thresholding/README.md

@@ -2,8 +2,9 @@
 
 This script allows to extract the prediction out of the detector _GeoJSON_ based on a series of provided threshold values.
 
-The prediction are first filtered based on their polygons' ares value, then through score value. Following, with the input of a Digital Elevation Model an elevation thresholding is processed. Finally the predicition polygons are joined (spatial join) with a distance criteria. 
+The first step going on is a clustering of the centroids of every prediction polygon. This is used as a way to maintain the scores for a further unary union of the polygons, that then uses the cluster value assigned as an aggregation method. This allows the removal of lower scores in a smarter way, *i.e.*, by maintaining the integrity of the final polygon. Then, predictions are filtered based on the polygons' areas value. Following, with the input of a Digital Elevation Model an elevation thresholding is processed. Finally, the predictions are filtered based on score. This score is calculated as the maximum score obtained from the polygons intersecting the merged polygons.
 
+It is important to mention that this sequence was modified on the 26th of November 2021. Predictions weren't once aesthetic, especially with smaller tiles. As the old procedure was used for the delivery of results before, it is for now maintaned in the past-filter branch of this repository.
 ## Usage
 
 The script expects then a prediction _GeoJSON_ file with all geometries containing a _score_ value normalised in _[0,1]_ and the filtering threshold :

tools/prediction-thresholding/src/prediction-thresholding.py

@@ -4,7 +4,7 @@
 #  prediction-thresholding
 #
 #      Nils Hamel - [email protected]
-#      Huriel Reichel - huriel.ruan@gmail.com
+#      Huriel Reichel - huriel.reichel@protonmail.com
 #      Alessandro Cerioni
 #      Copyright (c) 2020-2022 Republic and Canton of Geneva
 #
@@ -20,55 +20,52 @@
 #
 #  You should have received a copy of the GNU General Public License
 #  along with this program.  If not, see <http://www.gnu.org/licenses/>.
+
 import geopandas as gpd
-import rasterio
 import pandas as pd
+import rasterio
 import argparse
-import sys
+from sklearn.cluster import KMeans
 
+# argument parser
 pm_argparse = argparse.ArgumentParser()
 
-pm_argparse.add_argument( '-i', '--input', type=str  ,     help='input geojson path' )
-pm_argparse.add_argument( '-d', '--dem',   type=str  ,     help='input DEM path' )
-pm_argparse.add_argument( '-o', '--output',type=str  ,     help='output geojson path' )
-pm_argparse.add_argument( '-a', '--area' , type=float  ,     default = 1728.,  help='area threshold. Default to 1728' )
-pm_argparse.add_argument( '-s', '--score', type=float,     default = 0.9  , help='score threshold. Default to 0.9' )
-pm_argparse.add_argument( '-e', '--elevation', type=float, default = 1155 , help='elevation threshold. Default to 1155' )
-pm_argparse.add_argument( '--distance', type=float, default = 8, help="distance for union. Default to 8")
+pm_argparse.add_argument( '-i', '--input',     type=str,                         help='input geojson path' )
+pm_argparse.add_argument( '-d', '--dem',       type=str,                         help='input DEM path' )
+pm_argparse.add_argument( '-o', '--output',    type=str,                         help='output geojson path' )
+pm_argparse.add_argument( '-a', '--area' ,     type=float,     default = 1728.,  help='area threshold. Default to 1728' )
+pm_argparse.add_argument( '-s', '--score',     type=float,     default = 0.96  ,  help='score threshold. Default to 0.96' )
+pm_argparse.add_argument( '-e', '--elevation', type=float,     default = 1155. , help='elevation threshold. Default to 1155' )
+pm_argparse.add_argument( '--distance',        type=float,     default = 8,      help="distance for union. Default to 8")
 
 pm_args = pm_argparse.parse_args()
 
+# import predictions GeoJSON
 input = gpd.read_file(pm_args.input)
 input = input.to_crs(2056)
-
 total = len(input)
-input["area"] = input['geometry'].area/ 10**6
 
-input = input[input.area > pm_args.area]
-ta = len(input)
-ar = total - ta
-print(str(ar) + " predictions were removed by area threshold")
+# Centroid of every prediction polygon
+centroids = gpd.GeoDataFrame()
+centroids.geometry = input.representative_point()
+
+# KMeans Unsupervised Learning
+centroids = pd.DataFrame({'x': centroids.geometry.x, 'y': centroids.geometry.y})
+k = int( ( len(input) / 3 ) + 1 )
+cluster = KMeans(n_clusters=k, algorithm = 'auto', random_state = 1)
+model = cluster.fit(centroids)
+labels = model.predict(centroids)
+print("KMeans algorithm computed with k = " + str(k))
+
+# Dissolve and Aggregate
+input['cluster'] = labels
+input = input.dissolve(by = 'cluster', aggfunc = 'max')
+total = len(input)
 
+# filter by score
 input = input[input['score'] > pm_args.score]
-ts = len(input)
-sr = ta - ts
-print(str(sr) + " predictions were removed by score threshold")
-
-input['x'] = input.centroid.x
-input['y'] = input.centroid.y
-
-r = rasterio.open(pm_args.dem)
-
-row, col = r.index(input.x,input.y)
-values = r.read(1)[row,col]
-
-input['elev'] = values
-
-input = input[input['elev'] < pm_args.elevation]
-te = len(input)
-se = ts - te
-print(str(se) + " predictions were removed by elevation threshold")
-print(str(te) + " predictions left")
+sc = len(input)
+print(str(total - sc) + " predictions removed by score threshold")
 
 # Create empty data frame
 geo_merge = gpd.GeoDataFrame()
@@ -77,11 +74,23 @@
 geo_merge = input.buffer( +pm_args.distance, resolution = 2 )
 geo_merge = geo_merge.geometry.unary_union
 geo_merge = gpd.GeoDataFrame(geometry=[geo_merge], crs = input.crs )
-geo_merge = geo_merge.explode().reset_index(drop=True)
+geo_merge = geo_merge.explode(index_parts = True).reset_index(drop=True)
 geo_merge = geo_merge.buffer( -pm_args.distance, resolution = 2 )
-
 td = len(geo_merge)
-print(str(td) + " predictions left after joining resulting polygons")
+print(str(sc - td) + " difference to clustered predictions after union")
 
-geo_merge.to_file(pm_args.output, driver='GeoJSON')
+geo_merge = geo_merge[geo_merge.area > pm_args.area]
+ta = len(geo_merge)
+print(str(td - ta) + " predictions removed by area threshold")
+
+r = rasterio.open(pm_args.dem)
+row, col = r.index(geo_merge.centroid.x, geo_merge.centroid.y)
+values = r.read(1)[row,col]
+geo_merge.elev = values
+geo_merge = geo_merge[geo_merge.elev < pm_args.elevation]
+te = len(geo_merge)
+print(str(ta - te) + " predictions removed by elevation threshold")
 
+print(str(te) + " predictions left")
+
+geo_merge.to_file(pm_args.output, driver='GeoJSON')