diff --git a/Project.toml b/Project.toml
index 1173f1e1c..34e9264cc 100644
--- a/Project.toml
+++ b/Project.toml
@@ -28,6 +28,7 @@ Setfield = "efcf1570-3423-57d1-acb7-fd33fddbac46"
 ArchGDAL = "c9ce4bd3-c3d5-55b8-8973-c0e20141b8c3"
 CoordinateTransformations = "150eb455-5306-5404-9cee-2592286d6298"
 HDF5 = "f67ccb44-e63f-5c2f-98bd-6dc0ccc4ba2f"
+KernelDensity = "5ab0869b-81aa-558d-bb23-cbf5423bbe9b"
 Makie = "ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a"
 NCDatasets = "85f8d34a-cbdd-5861-8df4-14fed0d494ab"
 RasterDataSources = "3cb90ccd-e1b6-4867-9617-4276c8b2ca36"
@@ -36,6 +37,7 @@ RasterDataSources = "3cb90ccd-e1b6-4867-9617-4276c8b2ca36"
 RastersArchGDALExt = "ArchGDAL"
 RastersCoordinateTransformationsExt = "CoordinateTransformations"
 RastersHDF5Ext = "HDF5"
+RastersKernelDensityExt = "KernelDensity"
 RastersMakieExt = "Makie"
 RastersNCDatasetsExt = "NCDatasets"
 RastersRasterDataSourcesExt = "RasterDataSources"
@@ -54,6 +56,7 @@ Flatten = "0.4"
 GeoFormatTypes = "0.4"
 GeoInterface = "1"
 HDF5 = "0.14, 0.15, 0.16"
+KernelDensity = "0.6"
 Makie = "0.19, 0.20"
 Missings = "0.4, 1"
 NCDatasets = "0.10, 0.11, 0.12, 0.13"
@@ -74,6 +77,7 @@ CoordinateTransformations = "150eb455-5306-5404-9cee-2592286d6298"
 DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
 GeometryBasics = "5c1252a2-5f33-56bf-86c9-59e7332b4326"
 HDF5 = "f67ccb44-e63f-5c2f-98bd-6dc0ccc4ba2f"
+KernelDensity = "5ab0869b-81aa-558d-bb23-cbf5423bbe9b"
 Makie = "ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a"
 NCDatasets = "85f8d34a-cbdd-5861-8df4-14fed0d494ab"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
@@ -84,4 +88,4 @@ Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["Aqua", "ArchGDAL", "CFTime", "CoordinateTransformations", "DataFrames", "GeometryBasics", "HDF5", "NCDatasets", "Plots", "RasterDataSources", "SafeTestsets", "Shapefile", "Statistics", "Test"]
+test = ["Aqua", "ArchGDAL", "CFTime", "CoordinateTransformations", "DataFrames", "GeometryBasics", "HDF5", "KernelDensity", "NCDatasets", "Plots", "RasterDataSources", "SafeTestsets", "Shapefile", "Statistics", "Test"]
diff --git a/ext/RastersKernelDensityExt/RastersKernelDensityExt.jl b/ext/RastersKernelDensityExt/RastersKernelDensityExt.jl
new file mode 100644
index 000000000..3085bff37
--- /dev/null
+++ b/ext/RastersKernelDensityExt/RastersKernelDensityExt.jl
@@ -0,0 +1,82 @@
+module RastersKernelDensityExt
+
+using Extents,
+      GeoInterface,
+      KernelDensity,
+      Rasters
+
+using Rasters: Tables, maybeshiftlocus, Center, lookup, _extent2dims
+
+const GI = GeoInterface
+
+const VectorTuple = Tuple{<:AbstractVector{<:Number},<:AbstractVector{<:Number}}
+
+function Rasters.kerneldensity(geom; geomcolumn=nothing, kw...)
+    if Tables.istable(typeof(geom))
+        if geomcolumn isa Tuple
+            # Use the columns directly and skip _to_vectors
+            data = (Tables.getcolumn(geom, geomcolumn[1]), Tables.getcolumn(geom, geomcolumn[2]))
+            return _kerneldensity(data; kw...)
+        elseif isnothing(geomcolumn)
+            geomcolumn = GI.geometrycolumns(geom)[1]
+        end
+        geom = GI.getcolumn(geom, geomcolumn)
+    end
+    return _kerneldensity(_to_vectors(geom); kw...)
+end
+# Also accept the regular KernelDensity.jl data formats
+Rasters.kerneldensity(data::VectorTuple; kw...) = _kerneldensity(data; kw...)
+Rasters.kerneldensity(data::AbstractMatrix{<:Number}; kw...) =
+    _kerneldensity((data[:, 1], data[:, 2]); kw...)
+
+function _kerneldensity(data::VectorTuple;
+    to=nothing, size=nothing, res=nothing, crs=nothing, kw...
+)
+    # Get the extent from `data` if there is no `to` keyword
+    if isnothing(to)
+        x, y = map(extrema, data)
+        to = Extents.Extent(; X=x, Y=y)
+    end
+    # Convert input keywords to Dimensions
+    ds = _extent2dims(to; size, res, crs)
+
+    # KernelDensity needs midpoints
+    midpoint_ranges = map(parent ∘ lookup, maybeshiftlocus(Center(), ds))
+
+    # Create a bivariate kernel density estimator
+    bkde = KernelDensity.kde(data, midpoint_ranges; kw...)
+
+    # Return a raster
+    return Raster(bkde.density, ds)
+end
+
+_to_vectors(geoms) = _to_vectors(GI.trait(geoms), geoms)
+function _to_vectors(::Nothing, geoms::AbstractArray)
+    if all(GI.geomtrait(g) isa PointTrait for g in geoms)
+        return _fill_vectors(geoms, length(geoms))
+    elseif all(GI.geomtrait(g) isa MultiPointTrait for g in geoms)
+        n = sum(GI.npoint(g) for g in geoms)
+        points = Iterators.flatten(GI.getpoint(g) for g in geoms)
+        return _fill_vectors(points, n)
+    end
+end
+function _to_vectors(::GI.MultiPointTrait, geom)
+    n = GI.npoint(geom)
+    points = Iterators.flatten(GI.getpoint(g) for g in geoms)
+    return _fill_vectors(points, n)
+end
+_to_vectors(::PointTrait, geom) = [GI.x(geom), GI.y(geom)]
+function _to_vectors(trait, geom)
+    throw(ArgumentError("$trait not supported: only `MultiPointTrait` `PointTrait` and `MultiPointTrait` are supported"))
+end
+
+function _fill_vectors(points, n)
+    vec1 = Vector{Float64}(undef, n)
+    vec2 = Vector{Float64}(undef, n)
+    for (i, p) in enumerate(points)
+        vec1[i], vec2[i] = GI.x(p), GI.y(p)
+    end
+    return (vec1, vec2)
+end
+
+end
diff --git a/src/Rasters.jl b/src/Rasters.jl
index 64bc8131c..818edd1bc 100644
--- a/src/Rasters.jl
+++ b/src/Rasters.jl
@@ -59,7 +59,7 @@ export missingval, boolmask, missingmask, replace_missing, replace_missing!,
        aggregate, aggregate!, disaggregate, disaggregate!, mask, mask!,
        resample, warp, zonal, crop, extend, trim, slice, combine, points,
        classify, classify!, mosaic, mosaic!, extract, rasterize, rasterize!,
-       coverage, coverage!, setcrs, setmappedcrs, smapseries, cellsize
+       coverage, coverage!, setcrs, setmappedcrs, smapseries, cellsize, kerneldensity
 export crs, mappedcrs, mappedindex, mappedbounds, projectedindex, projectedbounds
 export reproject, convertlookup
 
diff --git a/src/extensions.jl b/src/extensions.jl
index 4dff7411e..bee5ce9b9 100644
--- a/src/extensions.jl
+++ b/src/extensions.jl
@@ -195,6 +195,41 @@ $EXPERIMENTAL
 """
 cellsize(args...; kw...) = throw_extension_error(cellsize, "ArchGDAL", :RastersArchGDALExt, args)
 
+
+"""
+    kerneldensity(geoms; kw...)
+
+Calculate the kernel density of the points in `geoms` using KernelDensity.jl.
+
+Returns a `Raster`.
+
+## Keywords
+
+- `geoms`: an iterable of GeoInterface.jl compatible points or multipoints,
+    or a Tables.jl compatible table with a geometry column. KernelDensity.jl
+    compatible 2-column `Matrix` or `Tuple` of `Vectors` are also supported.
+$TO_KEYWORD
+$RES_KEYWORD
+$SIZE_KEYWORD
+$CRS_KEYWORD
+$GEOMCOLUMN_KEYWORD
+- `kernel`: a distribution from Distributions.jl, default is `Normal()`.
+- `bandwidth` bandwidth to be passed to the kernel density estimator.
+- `weights`: weights to be passed to the kernel density estimator.
+
+## Example
+
+```julia
+using Rasters, KernelDensity
+points = zip(rand(Normal(), 100), rand(Normal(), 100))
+
+kerneldensity(points; res=0.1, weights=rand(100))
+```
+
+$EXPERIMENTAL
+"""
+kerneldensity(args...; kw...) = throw_extension_error(warp, "KernelDensity", :RastersKernelDensityExt, args)
+
 # Other shared stubs
 function layerkeys end
 function smapseries end
diff --git a/src/methods/rasterize.jl b/src/methods/rasterize.jl
index 916274eb8..7bd81e46c 100644
--- a/src/methods/rasterize.jl
+++ b/src/methods/rasterize.jl
@@ -382,6 +382,7 @@ These are detected automatically from `data` where possible.
 
 $GEOM_KEYWORDS
 $RASTERIZE_KEYWORDS
+$GEOMCOLUMN_KEYWORD
 $FILENAME_KEYWORD
 $SUFFIX_KEYWORD
 
diff --git a/src/methods/shared_docstrings.jl b/src/methods/shared_docstrings.jl
index f184534c4..fc2c30860 100644
--- a/src/methods/shared_docstrings.jl
+++ b/src/methods/shared_docstrings.jl
@@ -57,3 +57,9 @@ const PROGRESS_KEYWORD = """
 const VERBOSE_KEYWORD = """
 - `vebose`: whether to print messages about potential problems. `true` by default.
 """
+
+const GEOMCOLUMN_KEYWORD = """
+- `geomcolumn`: `Symbol` to manually select the column the geometries are in
+    when `data` is a Tables.jl compatible table, or a tuple of `Symbol` for columns
+    of point coordinates.
+"""
diff --git a/src/utils.jl b/src/utils.jl
index 3280fab59..7699c6cc2 100644
--- a/src/utils.jl
+++ b/src/utils.jl
@@ -128,15 +128,14 @@ function _without_mapped_crs(f, st::AbstractRasterStack, mappedcrs::GeoFormat)
     return x
 end
 
-function _extent2dims(to; size=nothing, res=nothing, crs=nothing, kw...) 
-    _extent2dims(to, size, res, crs)
-end
-function _extent2dims(to::Extents.Extent, size::Nothing, res::Nothing, crs)
-    isnothing(res) && throw(ArgumentError("Pass either `size` or `res` keywords or a `Tuple` of `Dimension`s for `to`."))
-end
-function _extent2dims(to::Extents.Extent, size, res, crs)
-    isnothing(res) || _size_and_res_error()
-end
+_extent2dims(to; size=nothing, res=nothing, crs=nothing, kw...) = _extent2dims(to, size, res, crs)
+_extent2dims(to::DimTuple, size::Nothing, res::Nothing, crs) = setcrs(to, crs)
+_extent2dims(to::DimTuple, size::Nothing, res::Nothing, crs::Nothing) = to
+_extent2dims(to, size::Nothing, res::Nothing, crs) = _extent2dims(dims(to), size, res, crs)
+_extent2dims(to, size, res, crs) = _extent2dims(Extents.extent(to), size, res, crs)
+_extent2dims(to::Extents.Extent, size, res, crs) = _size_and_res_error()
+_extent2dims(to::Union{Nothing,Extents.Extent}, size::Nothing, res::Nothing, crs) =
+    throw(ArgumentError("Pass either `size` or `res` keywords or a `Tuple` of `Dimension`s for `to`."))
 function _extent2dims(to::Extents.Extent{K}, size::Nothing, res::Real, crs) where K
     tuple_res = ntuple(_ -> res, length(K))
     _extent2dims(to, size, tuple_res, crs)
diff --git a/test/kerneldensity.jl b/test/kerneldensity.jl
new file mode 100644
index 000000000..6e589f723
--- /dev/null
+++ b/test/kerneldensity.jl
@@ -0,0 +1,36 @@
+using Test, Rasters, KernelDensity, Extents
+
+points = collect(zip(rand(100) .* 10, rand(100)))
+# Using just res, detecting extent from points
+res_rast = kerneldensity(points; res=(0.1, 0.01))
+# This doesn't exactly match, someting is wrong with `_extent2dims`
+@test_broken map(step, dims(res_rast)) == (0.1, 0.01)
+# Using just size, detecting extent from points
+size_rast = kerneldensity(points; size=250)
+@test size(size_rast) == (250, 250)
+# Using a defined extent with res
+to = Extents.Extent(X=(0, 10), Y=(0, 5))
+ext_res_rast = kerneldensity(points; to, res=0.05)
+@test size(ext_rast) == (200, 100)
+@test map(step, dims(ext_rast)) == (0.05, 0.05)
+# Using a defined extent with size
+to = Extents.Extent(X=(0, 10), Y=(0, 5))
+ext_res_rast = kerneldensity(points; to, size=(200, 100))
+@test size(ext_rast) == (200, 100)
+@test map(step, dims(ext_rast)) == (0.05, 0.05)
+# Using a defined DimArray
+to = rand(X(0:0.01:10), Y(0:0.001:1))
+template_rast = kerneldensity(points; to)
+@test size(to) == size(template_rast)
+@test map(step, dims(template_rast)) == (0.01, 0.001)
+
+# Test plots
+# using GLMakie
+# p = Makie.heatmap(res_rast)
+# Makie.scatter!(points)
+# p = Makie.heatmap(size_rast)
+# Makie.scatter!(points)
+# p = Makie.heatmap(ext_rast)
+# Makie.scatter!(points)
+# p = Makie.heatmap(template_rast)
+# Makie.scatter!(points)
diff --git a/test/runtests.jl b/test/runtests.jl
index eaaf8d944..6e50b771a 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -23,7 +23,7 @@ end
 @time @safetestset "reproject" begin include("reproject.jl") end
 @time @safetestset "warp" begin include("warp.jl") end
 @time @safetestset "resample" begin include("resample.jl") end
-
+@time @safetestset "kernel density" begin include("kerneldensity.jl") end
 
 # Only test SMAP locally for now, also RasterDataSources because CI dowloads keep breaking
 @time @safetestset "ncdatasets" begin include("sources/ncdatasets.jl") end