diff --git a/Project.toml b/Project.toml
index f6d5737842..0a4ec11072 100644
--- a/Project.toml
+++ b/Project.toml
@@ -19,6 +19,7 @@ Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 SortTileRecursiveTree = "746ee33f-1797-42c2-866d-db2fce69d14d"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
+StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
 Tables = "bd369af6-aec1-5ad0-b16a-f7cc5008161c"
 
 [weakdeps]
@@ -57,6 +58,7 @@ Random = "1"
 SortTileRecursiveTree = "0.1.2"
 StaticArrays = "1"
 Statistics = "1"
+StatsBase = "0.34"
 TGGeometry = "0.1"
 Tables = "1"
 julia = "1.10"
diff --git a/src/GeometryOps.jl b/src/GeometryOps.jl
index 2000e0d915..e8f25dcdb3 100644
--- a/src/GeometryOps.jl
+++ b/src/GeometryOps.jl
@@ -27,6 +27,7 @@ using StaticArrays
 import Tables, DataAPI
 import StaticArrays
 import DelaunayTriangulation # for convex hull and triangulation
+import StatsBase
 import ExactPredicates
 import Base.@kwdef
 import GeoInterface.Extents: Extents
@@ -85,6 +86,7 @@ include("methods/geom_relations/overlaps.jl")
 include("methods/geom_relations/touches.jl")
 include("methods/geom_relations/within.jl")
 include("methods/geom_relations/common.jl")
+include("methods/sample.jl")
 include("methods/orientation.jl")
 include("methods/polygonize.jl")
 include("methods/voronoi.jl")
diff --git a/src/methods/sample.jl b/src/methods/sample.jl
new file mode 100644
index 0000000000..5ac61b2fb3
--- /dev/null
+++ b/src/methods/sample.jl
@@ -0,0 +1,60 @@
+#=
+# Sample
+
+```@example sample
+import GeoInterface as GI, GeometryOps as GO
+p1 = GI.Polygon([[[-55965.680060140774, -31588.16072168928], [-55956.50771556479, -31478.09258677756], [-31577.548550575284, -6897.015828572996], [-15286.184961223798, -15386.952072224134], [-9074.387601621409, -27468.20712382156], [-8183.4538916097845, -31040.003969070774], [-27011.85123029944, -38229.02388009402], [-54954.72822634951, -32258.9734800704], [-55965.680060140774, -31588.16072168928]]])
+points = GO.sample(p1, 100)
+using CairoMakie
+f, a, p = poly(p1)
+scatter!(a, points)
+f
+```
+
+=#
+
+export sample, UniformSampling
+
+struct UniformSampling
+end
+
+application_level(::UniformSampling) = TraitTarget(GI.MultiPolygonTrait(), GI.MultiLineStringTrait(), GI.MultiPointTrait(), GI.PolygonTrait(), GI.LineStringTrait())
+
+function sample(geom, n::Int)
+    return sample(UniformSampling(), geom, n)
+end
+
+function sample(alg, geom, n)
+    return apply(x -> _sample(alg, GI.trait(x), x, n), application_level(alg), geom)
+end
+
+function _sample(alg::UniformSampling, ::Union{GI.PolygonTrait, GI.MultiPolygonTrait}, geom, n)
+    (; X, Y) = GI.extent(geom)
+    points = fill((0.0, 0.0), n)
+    i = 1
+    while i <= n
+        x = rand() * (X[2] - X[1]) + X[1]
+        y = rand() * (Y[2] - Y[1]) + Y[1]
+        if contains(geom, (x, y))
+            points[i] = (x, y)
+            i += 1
+        end
+    end
+    return points
+end
+
+function _sample(alg::UniformSampling, ::GI.LineStringTrait, geom, n)
+    edges = to_edges(geom)
+    edge_lengths = map(splat(distance), edges)
+    # normalize the vector
+    edge_probabilities = edge_lengths ./ sum(edge_lengths)
+    edge_idxs = 1:length(edges)
+    return map(1:n) do _
+        edge_idx = sample(edge_idxs, edge_probabilities)
+        x1, y1 = edges[edge_idx][1]
+        x2, y2 = edges[edge_idx][2]
+        distance = edge_lengths[edge_idx]
+        t = rand() * distance
+        (x1 + t * (x2 - x1), y1 + t * (y2 - y1))
+    end
+end
\ No newline at end of file