diff --git a/Project.toml b/Project.toml
index 6f20c7e5f..2b992fe01 100644
--- a/Project.toml
+++ b/Project.toml
@@ -12,6 +12,7 @@ InteractiveUtils = "b77e0a4c-d291-57a0-90e8-8db25a27a240"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 SortTileRecursiveTree = "746ee33f-1797-42c2-866d-db2fce69d14d"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
+StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
 Tables = "bd369af6-aec1-5ad0-b16a-f7cc5008161c"
 
 [weakdeps]
@@ -35,6 +36,7 @@ LinearAlgebra = "1"
 Proj = "1"
 SortTileRecursiveTree = "0.1"
 Statistics = "1"
+StatsBase = "0.34"
 Tables = "1"
 julia = "1.9"
 
diff --git a/src/GeometryOps.jl b/src/GeometryOps.jl
index 1aced0e82..424de9d35 100644
--- a/src/GeometryOps.jl
+++ b/src/GeometryOps.jl
@@ -6,6 +6,7 @@ using GeoInterface
 using GeometryBasics
 import Tables
 using LinearAlgebra, Statistics
+import StatsBase
 import GeometryBasics.StaticArrays
 import ExactPredicates
 import Base.@kwdef
@@ -47,6 +48,7 @@ include("methods/geom_relations/intersects.jl")
 include("methods/geom_relations/overlaps.jl")
 include("methods/geom_relations/touches.jl")
 include("methods/geom_relations/within.jl")
+include("methods/sample.jl")
 include("methods/orientation.jl")
 include("methods/polygonize.jl")
 
diff --git a/src/methods/sample.jl b/src/methods/sample.jl
new file mode 100644
index 000000000..5ac61b2fb
--- /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