Sinuosity of Kemijoki: Difference between revisions

From wikiluntti
Line 22: Line 22:


<syntaxhighlight lang="Python">
<syntaxhighlight lang="Python">
import geopandas as gpd
finland = gpd.read_file("fi.shp")
finland = finland.to_crs(epsg=3067)
river = gpd.read_parquet("kemijoki_mainbranch.parquet")
</syntaxhighlight>
</syntaxhighlight>
Convert to points


<syntaxhighlight lang="Python">
<syntaxhighlight lang="Python">
from shapely.geometry import Point
line = river.geometry.iloc[0]
current_point = Point(line.coords[0])
points = [current_point]
</syntaxhighlight>
</syntaxhighlight>


=== ===
=== ===

Revision as of 11:44, 27 April 2026

Introduction

Use Python to estimate the sinuosity of Kemijoki.

The river was obtained from https://download.geofabrik.de/europe/finland.html but only one river (Kemijoki) with one branch is used. The extracting method is described elsewhere.

Kemijoki

Load Data Files

import geopandas as gpd

finland = gpd.read_file("fi.shp")
finland = finland.to_crs(epsg=3067)
river = gpd.read_parquet("kemijoki_mainbranch.parquet")

Convert to points

from shapely.geometry import Point

line = river.geometry.iloc[0]
current_point = Point(line.coords[0])
points = [current_point]

Theory

Douglas-Peucker simplification at different tolerances

Fractal Dimension

The fractal dimension is defined by

which gives

The corresponding plot is called a Richardson plot.

from scipy.stats import linregress
import numpy as np

#Plot everything
plt.plot( np.log(df["R"]), np.log(df["Length (km)"])) #x, y
plt.xlabel("log (R [km])")
plt.ylabel("log(Length [km])")

#Regression
Rvalue = [2, 6, 8, 10, 13]
xmin = Rvalue[0]
for xmax in Rvalue[1:]:
    mask = (df["R"] > np.exp(xmin)) & (df["R"] < np.exp(xmax))
    x = np.log(df.loc[mask, "R"].values)             # log(epsilon)
    y = np.log(df.loc[mask, "Length (km)"].values)   # log(length)
    slope, intercept, _, _, _ = linregress(x, y)
    x_line = np.linspace(x[0], x[-1], 100)
    y_line = slope * x_line + intercept
    plt.plot(x_line, y_line, color='red', label="Fit line")
    plt.text(x[0], y[0], 'd='+ str(round(1-slope, 2)))

    print( xmin, xmax )
    xmin = xmax

plt.title("log-log Length vs R")
plt.show()

References

https://www.tud.ttu.ee/web/dmitri.kartofelev/YFX1560/Loeng_14no_vid.pdf

https://preview-www.nature.com/articles/s41598-021-85405-0.pdf