statistics.py

import seaborn as sns
import json

f = open("settings.json", "r")
settings = json.load(f)

def intersecting_points():
    cm_file = open(cm_paths[0])
    cm = cityjson.reader(file=cm_file, ignore_duplicate_keys=True)
    
    for pc_path in pc_paths:
        infile = laspy.file.File(pc_path)
        
        scale = infile.header.scale

        points = np.vstack([infile.X * scale[0], infile.Y * scale[1]]).transpose()
        tree = spatial.KDTree(points)
        
        intersects = {}
        
        for co_id in cm.j["CityObjects"]:
            try:
                print(i)
                intersect_count = 0
                co_type = cm.j["CityObjects"][co_id]["type"]
                if co_type not in intersects.keys():
                    intersects[co_type] = []
                
                triangles = []
                
                for geom in cm.j["CityObjects"][co_id]["geometry"]:
                    geoms_to_points(geom["boundaries"], triangles, cm)
    
                triangles_2d = np.array(triangles)[:,:,:2]
                
                shapely_triangles = [geometry.Polygon(t_2d) for t_2d in triangles_2d]
    
                polygon = ops.unary_union(shapely_triangles)
                    
                bbox = polygon.bounds
                bbox_centre = [(bbox[2] - bbox[0]) / 2 + bbox[0], (bbox[3] - bbox[1]) / 2 + bbox[1]]
                bbox_extent = max(bbox[2] - bbox[0], bbox[3] - bbox[1])
                
                candidates = tree.query_ball_point(bbox_centre, bbox_extent)
                    
                for p_i in candidates:
                    p = geometry.Point(points[p_i])
                    if p.intersects(polygon):
                        intersect_count += 1
                        
                intersects[co_type].append(intersect_count)

            except:
                continue

#intersecting_points()
#terrain_diff(cm_paths)
#print(set(d["type"]))
#df = pd.DataFrame(d)
#df["difference"] = abs(df["cm1"] - df["cm2"])

# visualise relationship between object height and height difference
#plt.scatter(x=df["cm1"], y=df["difference"], s=5)
#plt.xlabel("Object height")
#plt.ylabel("Height difference")

#print("Mean difference = " + str(df["difference"].mean()))
#print("Max difference = " + str(df["difference"].max()))

#print(df.groupby(['type'])["difference"].mean())
#print(df.groupby(['type'])["difference"].max())
                


def figure_appearance():
    sns.set(font='Franklin Gothic Book',
            rc={
     'axes.axisbelow': False,
     'axes.edgecolor': 'lightgrey',
     'axes.facecolor': 'None',
     'axes.grid': False,
     'axes.labelcolor': 'dimgrey',
     'axes.spines.right': False,
     'axes.spines.top': False,
     'figure.facecolor': 'white',
     'lines.solid_capstyle': 'round',
     'patch.edgecolor': 'w',
     'patch.force_edgecolor': True,
     'text.color': 'dimgrey',
     'xtick.bottom': False,
     'xtick.color': 'dimgrey',
     'xtick.direction': 'out',
     'xtick.top': False,
     'ytick.color': 'dimgrey',
     'ytick.direction': 'out',
     'ytick.left': False,
     'ytick.right': False})
     
    sns.set_context("notebook", rc={"font.size":12,
                                    "axes.titlesize":16,
                                    "axes.labelsize":14})
    
    df_input = {"cm1_height": cm1_height, "cm2_height": cm2_height, "difference": difference}
    
    df = pd.DataFrame(df_input)
    
    CB91_Blue = '#2CBDFE'
    CB91_Green = '#47DBCD'
    CB91_Pink = '#F3A0F2'
    CB91_Purple = '#9D2EC5'
    CB91_Violet = '#661D98'
    CB91_Amber = '#F5B14C'

plt.scatter(x=df["cm1_height"], y=abs(df["difference"]), s=5, color=CB91_Blue)
plt.xlabel("Building height (dense image matching)")
plt.ylabel("Absolute height difference")
plt.show()

print("Mean difference = " + str(df["difference"].mean()))
print("Max difference = " + str(df["difference"].max()))

labels = ["0-0.1", "0.1-0.5", "0.5-1", "1-1.5", "1.5-2", ">=2"]

bins = [0, 0.2, 1, 2, 3, 4, 5]
binned = pd.cut(df["difference"], bins, include_lowest=True)
ax = binned.value_counts(sort=False, normalize=True).plot.bar(rot=0, color=CB91_Blue)
#ax.set_xticklabels([str(c.left) + " - " + str(c.right) for c in binned.cat.categories])
plt.xticks(np.arange(len(labels)), labels)
plt.ylabel("Frequency (normalised)")
plt.xlabel("Absolute height difference (m)")
plt.title("Height difference per building")
#plt.show()

plt.savefig('buildings.pdf')
plt.savefig('buildings.jpg', dpi=200)

data = {}
bins = [0, 0.1, 0.5, 1, 1.5, 2, 10]
#labels = ["[0, 0.1>", "[0.1, 0.5>", "[0.5, 1>", "[1, 1.5>", "[1.5, 2>", ">2="]
labels = ["0-0.1", "0.1-0.5", "0.5-1", "1-1.5", "1.5-2", ">=2"]
xs = np.arange(len(labels))  # the label locations
width = 0.25  # the width of the bars
xs = [n - width * len(data.keys()) for n in xs]

colors = [CB91_Blue, CB91_Green, CB91_Pink, CB91_Purple, CB91_Violet, CB91_Amber]

for root, dirs, files in os.walk("."):
    for file in files:
        if "_diff" in file and file[-3:] == "asc":
            fn = file.split("_")[0]
            print(fn)
            r = rasterio.open(file)
            b = r.read(1)
            i = np.where(b != -99)
            data[fn] = b[i]
            print(fn)
            print("Mean difference = " + str(data[fn].mean()))
            print("Max difference = " + str(data[fn].max()))

fig, ax = plt.subplots()


for i, (k,v) in enumerate(data.items()):
    binned = pd.cut(v, bins, include_lowest=True).value_counts().values

    # percentages
    binned = [(n / sum(binned)) for n in binned]
    x = [n + width * i - width for n in xs]
    ax.bar(x, binned, width, color=colors[i], label=k)

plt.xticks(np.arange(len(labels)), labels)
plt.ylabel("Frequency (normalised)")
plt.xlabel("Absolute height difference (m)")
plt.title("Height difference per raster cell (1m²)")
plt.legend()
#plt.show()
plt.savefig('other_features.pdf')
plt.savefig('other_features.jpg', dpi=200)