diff --git a/CHANGELOG.md b/CHANGELOG.md
index f2c67be..48f5ad2 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -7,6 +7,27 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 ## [Unreleased]
 
+## [0.9.0] - 2024-06-11
+
+### Changes
+- [Breaking] Changed type signature of GeoArray to include the number of dimensions. This allows
+for single "band" GeoArrays to be represented as matrices as well. This should make it easier to 
+work with single band rasters and the image ecosystem.
+- [Breaking] `bbox` and friends now return an `Extents.Extent` instead of a `NamedTuple`.
+- [Breaking] Reverted rename of `equals` to `Base.isequal`, now called `isgeoequal`.
+- [Breaking] `getindex`, `indices` and `sample` now use the rounding mode `RoundNearestTiesUp` instead of `RoundNearest`.
+
+### Deprecated
+- [Breaking] Bounding box input for `crop`, `warp` and others is now deprecated. Use an `Extent` instead.
+
+### Fixed
+- Try to release file lock on close as soon as possible.
+- Indexing a GeoArray with `[:, :]` now returns a GeoArray
+
+### Added
+- Added `enable_online_warp` function to enable online access to PROJ data for `warp`.
+- Added `rounding` argument to `indices` to control rounding mode used.
+
 ## [0.8.5] - 2024-01-07
 - Fix small bug in metadata
 - Move GeoStatsBase into an extension
@@ -70,6 +91,11 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 - Fixed iterate specification so `sum` on a GeoArray is correct
 
 [unreleased]: https://github.com/evetion/GeoArrays.jl/compare/v0.8.1...HEAD
+[0.9.0]: https://github.com/evetion/GeoArrays.jl/compare/v0.8.5...v0.9.0
+[0.8.5]: https://github.com/evetion/GeoArrays.jl/compare/v0.8.4...v0.8.5
+[0.8.4]: https://github.com/evetion/GeoArrays.jl/compare/v0.8.3...v0.8.4
+[0.8.3]: https://github.com/evetion/GeoArrays.jl/compare/v0.8.2...v0.8.3
+[0.8.2]: https://github.com/evetion/GeoArrays.jl/compare/v0.8.1...v0.8.2
 [0.8.1]: https://github.com/evetion/GeoArrays.jl/compare/v0.8.0...v0.8.1
 [0.8.0]: https://github.com/evetion/GeoArrays.jl/compare/v0.7.13...v0.8.0
 [0.7.13]: https://github.com/evetion/GeoArrays.jl/compare/v0.7.12...v0.7.13
diff --git a/Project.toml b/Project.toml
index 54d5c32..41a307b 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 name = "GeoArrays"
 uuid = "2fb1d81b-e6a0-5fc5-82e6-8e06903437ab"
 authors = ["Maarten Pronk <git@evetion.nl>"]
-version = "0.8.5"
+version = "0.9.0"
 
 [weakdeps]
 GeoStatsBase = "323cb8eb-fbf6-51c0-afd0-f8fba70507b2"
diff --git a/README.md b/README.md
index af3833d..38d5ac7 100644
--- a/README.md
+++ b/README.md
@@ -11,7 +11,7 @@ This packages takes its inspiration from Python's [rasterio](https://github.com/
 ## Installation
 
 ```julia
-(v1.8) pkg> add GeoArrays
+(v1.10) pkg> add GeoArrays
 ```
 
 ## Examples
@@ -30,7 +30,7 @@ Read a GeoTIFF file and display its information, i.e. AffineMap and projection (
 # Read TIF file
 julia> fn = download("https://github.com/yeesian/ArchGDALDatasets/blob/master/data/utmsmall.tif?raw=true")
 julia> geoarray = GeoArrays.read(fn)
-100x100x1 Array{UInt8,3} with AffineMap([60.0 0.0; 0.0 -60.0], [440720.0, 3.75132e6]) and CRS PROJCS["NAD27 / UTM zone 11N"...
+100x100 Matrix{UInt8} with AffineMap([60.0 0.0; 0.0 -60.0], [440720.0, 3.75132e6]) and CRS PROJCS["NAD27 / UTM zone 11N"...
 
 # Affinemap containing offset and scaling
 julia> geoarray.f
@@ -152,10 +152,11 @@ For example, we can vertically transform from the ellipsoid
 to the EGM2008 geoid using EPSG code 3855.
 Note that the underlying PROJ library needs to find the geoidgrids,
 so if they're not available locally, one needs to set `ENV["PROJ_NETWORK"] = "ON"`
-as early as possible, ideally before loading GeoArrays.
+as early as possible, ideally before loading GeoArrays, or use `enable_online_warp`.
 ```julia
+enable_online_warp()  # If you don't have PROJ grids locally
 ga = GeoArray(zeros((360, 180)))
-bbox!(ga, (min_x=-180, min_y=-90, max_x=180, max_y=90))
+bbox!(ga, Extent(X(-180, 180), Y(-90, 90)))
 crs!(ga, GeoFormatTypes.EPSG(4979))  # WGS83 in 3D (reference to ellipsoid)
 ga2 = GeoArrays.warp(ga, Dict("t_srs" => "EPSG:4326+3855"))
 ```
diff --git a/src/GeoArrays.jl b/src/GeoArrays.jl
index 1e65b38..7d068bb 100644
--- a/src/GeoArrays.jl
+++ b/src/GeoArrays.jl
@@ -7,7 +7,7 @@ using StaticArrays
 import GeoInterface as GI
 using IterTools: partition
 import DataAPI
-import Extents
+using Extents: Extent, intersects
 using PrecompileTools    # this is a small dependency
 
 include("geoarray.jl")
diff --git a/src/geoarray.jl b/src/geoarray.jl
index 334d7b4..b42c97a 100644
--- a/src/geoarray.jl
+++ b/src/geoarray.jl
@@ -1,13 +1,14 @@
 const NumberOrMissing = Union{Number,Union{Missing,Number}}
+const MatrixorArray = Union{<:AbstractArray{T,2},<:AbstractArray{T,3}} where {T}
 
 """
-    GeoArray{T::NumberOrMissing,A<:AbstractArray{T,3}} <: AbstractArray{T,3}
+    GeoArray{T::NumberOrMissing,N,A<:AbstractArray{T,N}} <: AbstractArray{T,N}
 
 A GeoArray is an AbstractArray, an AffineMap for calculating coordinates based on the axes and a CRS definition to interpret these coordinates into in the real world.
 It's three dimensional and can be seen as a stack (3D) of 2D geospatial rasters (bands), the dimensions are :x, :y, and :bands.
 The AffineMap and CRS (coordinates) only operate on the :x and :y dimensions.
 """
-mutable struct GeoArray{T<:NumberOrMissing,A<:AbstractArray{T,3}} <: AbstractArray{T,3}
+mutable struct GeoArray{T<:NumberOrMissing,N,A<:AbstractArray{T,N}} <: AbstractArray{T,N}
     A::A
     f::CoordinateTransformations.AffineMap{StaticArrays.SMatrix{2,2,Float64,4},StaticArrays.SVector{2,Float64}}
     crs::GFT.WellKnownText{GFT.CRS}
@@ -15,7 +16,7 @@ mutable struct GeoArray{T<:NumberOrMissing,A<:AbstractArray{T,3}} <: AbstractArr
 end
 
 """
-    GeoArray(A::AbstractArray{T,3} where T <: NumberOrMissing)
+    GeoArray(A::AbstractArray{T,2|3} where T <: NumberOrMissing)
 
 Construct a GeoArray from any Array. A default `AffineMap` and `CRS` will be generated.
 
@@ -25,45 +26,31 @@ julia> GeoArray(rand(10,10,1))
 10x10x1 Array{Float64, 3} with AffineMap([1.0 0.0; 0.0 1.0], [0.0, 0.0]) and undefined CRS
 ```
 """
-GeoArray(A::AbstractArray{T,3} where {T<:NumberOrMissing}) = GeoArray(A, geotransform_to_affine(SVector(0.0, 1.0, 0.0, 0.0, 0.0, 1.0)), "", Dict{String,Any}())
+GeoArray(A::MatrixorArray where {T<:NumberOrMissing}) = GeoArray(A, geotransform_to_affine(SVector(0.0, 1.0, 0.0, 0.0, 0.0, 1.0)), "", Dict{String,Any}())
 """
     GeoArray(A::AbstractArray{T,3} where T <: NumberOrMissing, f::AffineMap)
 
 Construct a GeoArray from any Array and an `AffineMap` that specifies the coordinates. A default `CRS` will be generated.
 """
-GeoArray(A::AbstractArray{T,3} where {T<:NumberOrMissing}, f::AffineMap) = GeoArray(A, f, GFT.WellKnownText(GFT.CRS(), ""), Dict{String,Any}())
+GeoArray(A::MatrixorArray where {T<:NumberOrMissing}, f::AffineMap) = GeoArray(A, f, GFT.WellKnownText(GFT.CRS(), ""), Dict{String,Any}())
 
 """
-    GeoArray(A::AbstractArray{T,3} where T <: NumberOrMissing, f::AffineMap, crs::String)
+    GeoArray(A::AbstractArray{T,2|3} where T <: NumberOrMissing, f::AffineMap, crs::String)
 
 Construct a GeoArray from any Array and an `AffineMap` that specifies the coordinates and `crs` string in WKT format.
 """
-GeoArray(A::AbstractArray{T,3} where {T<:NumberOrMissing}, f::AffineMap, crs::String) = GeoArray(A, f, GFT.WellKnownText(GFT.CRS(), crs), Dict{String,Any}())
-GeoArray(A::AbstractArray{T,3} where {T<:NumberOrMissing}, f::AffineMap, crs::String, d::Dict{String}) = GeoArray(A, f, GFT.WellKnownText(GFT.CRS(), crs), d)
-GeoArray(A::AbstractArray{T,3} where {T<:NumberOrMissing}, f::AffineMap, crs::GFT.WellKnownText{GFT.CRS}) = GeoArray(A, f, crs, Dict{String,Any}())
-
-GeoArray(A::AbstractArray{T,3} where {T<:NumberOrMissing}, f::AffineMap{Matrix{Float64},Vector{Float64}}, crs::GFT.WellKnownText{GFT.CRS}) = GeoArray(A, AffineMap(SMatrix{2,2}(f.linear), SVector{2}(f.translation)), crs, Dict{String,Any}())
-GeoArray(A::AbstractArray{T,3} where {T<:NumberOrMissing}, f::AffineMap{Matrix{Float64},Vector{Float64}}, crs::GFT.WellKnownText{GFT.CRS}, d::Dict{String}) = GeoArray(A, AffineMap(SMatrix{2,2}(f.linear), SVector{2}(f.translation)), crs, d)
-
-"""
-    GeoArray(A::AbstractArray{T,2} where T <: NumberOrMissing)
-
-Construct a GeoArray from any Matrix, a third singleton dimension will be added automatically.
-
-# Examples
-```julia-repl
-julia> GeoArray(rand(10,10))
-10x10x1 Array{Float64, 3} with AffineMap([1.0 0.0; 0.0 1.0], [0.0, 0.0]) and undefined CRS
-```
-"""
-GeoArray(A::AbstractArray{T,2} where {T<:NumberOrMissing}, args...) = GeoArray(reshape(A, size(A)..., 1), args...)
+GeoArray(A) = GeoArray(A, f, GFT.WellKnownText(GFT.CRS(), crs), Dict{String,Any}())
+GeoArray(A, f, crs::String, d) = GeoArray(A, f, GFT.WellKnownText(GFT.CRS(), crs), d)
+GeoArray(A, f, crs) = GeoArray(A, f, crs, Dict{String,Any}())
+GeoArray(A, f::AffineMap{Matrix{Float64},Vector{Float64}}, args...) = GeoArray(A, AffineMap(SMatrix{2,2}(f.linear), SVector{2}(f.translation)), args...)
+GeoArray(ga::T, args...) where {T<:GeoArray} = GeoArray(parent(ga), args...)
 
 """
-    GeoArray(A::AbstractArray{T,3} where T <: NumberOrMissing, x::AbstractRange, y::AbstractRange, args...)
+    GeoArray(A::AbstractArray{T,2|3} where T <: NumberOrMissing, x::AbstractRange, y::AbstractRange, args...)
 
 Construct a GeoArray any Array and it's coordinates from `AbstractRange`s for each dimension.
 """
-function GeoArray(A::AbstractArray{T,3} where {T<:NumberOrMissing}, x::AbstractRange, y::AbstractRange, args...)
+function GeoArray(A::MatrixorArray where {T<:NumberOrMissing}, x::AbstractRange, y::AbstractRange, args...)
     size(A)[1:2] != (length(x), length(y)) && error("Size of `GeoArray` $(size(A)) does not match size of (x,y): $((length(x), length(y))). Note that this function takes *center coordinates*.")
     f = unitrange_to_affine(x, y)
     GeoArray(A, f, args...)
@@ -71,16 +58,30 @@ end
 
 # Behave like an Array
 Base.size(ga::GeoArray) = size(ga.A)
+_size(ga::GeoArray{T,2}) where {T} = (size(ga.A)..., 1)
+_size(ga::GeoArray{T,3}) where {T} = size(ga.A)
 Base.IndexStyle(::Type{<:GeoArray}) = IndexCartesian()
-Base.similar(ga::GeoArray, t::Type) = GeoArray(similar(ga.A, t), ga.f, ga.crs, ga.metadata)
+Base.similar(ga::GeoArray, t::Type) = GeoArray(similar(parent(ga), t), ga.f, ga.crs, ga.metadata)
+function Base.similar(ga::GeoArray, A::MatrixorArray)
+    size(ga.A)[1:2] == size(A)[1:2] || error(lazy"Size of `GeoArray` $(size(ga.A)) does not match size of `A`: $(size(A)).")
+    GeoArray(A, ga.f, ga.crs, ga.metadata)
+end
+Base.similar(ga::GeoArray, A::GeoArray) = similar(ga, parent(A))
 Base.iterate(ga::GeoArray) = iterate(ga.A)
 Base.iterate(ga::GeoArray, state) = iterate(ga.A, state)
 Base.length(ga::GeoArray) = length(ga.A)
 Base.parent(ga::GeoArray) = ga.A
-Base.eltype(::Type{GeoArray{T}}) where {T} = T
+Base.eltype(::Type{GeoArray{T,N}}) where {T,N} = T
 Base.show(io::IO, ::MIME"text/plain", ga::GeoArray) = show(io, ga)
 
-Base.convert(::Type{GeoArray{T}}, ga::GeoArray) where {T} = GeoArray(convert(Array{T}, ga.A), ga.f, ga.crs, ga.metadata)
+function Base.convert(::Type{GeoArray{T}}, ga::GeoArray{X,N}) where {T,N,X}
+    A = convert(Array{T,N}, ga.A)
+    GeoArray(A, ga.f, ga.crs, ga.metadata)
+end
+function Base.convert(::Type{GeoArray{T,N}}, ga::GeoArray) where {T,N}
+    A = convert(Array{T,N}, ga.A)
+    GeoArray(A, ga.f, ga.crs, ga.metadata)
+end
 
 Base.BroadcastStyle(::Type{<:GeoArray}) = Broadcast.ArrayStyle{GeoArray}()
 function Base.similar(bc::Broadcast.Broadcasted{Broadcast.ArrayStyle{GeoArray}}, ::Type{ElType}) where {ElType}
@@ -114,17 +115,26 @@ julia> ga[2:3,2:3,1]
 2x2x1 Array{Float64, 3} with AffineMap([1.0 0.0; 0.0 1.0], [1.0, 1.0]) and undefined CRS
 ```
 """
-function Base.getindex(ga::GeoArray, i::AbstractRange, j::AbstractRange, k::Union{Colon,AbstractRange,Integer})
-    A = getindex(ga.A, i, j, k)
+function Base.getindex(ga::GeoArray{T,N}, i::AbstractRange, j::AbstractRange, k::Union{Colon,AbstractRange,Integer}=Colon()) where {T,N}
+    A = N == 2 ? getindex(ga.A, i, j) : getindex(ga.A, i, j, k)
     x, y = first(i) - 1, first(j) - 1
     t = ga.f(SVector(x, y))
     l = ga.f.linear * SMatrix{2,2}([step(i) 0; 0 step(j)])
     GeoArray(A, AffineMap(l, t), crs(ga), ga.metadata)
 end
-Base.getindex(ga::GeoArray, i::AbstractRange, j::AbstractRange) = Base.getindex(ga, i, j, :)
+function Base.getindex(ga::GeoArray{T,N}, i::Colon, j::Colon, k::Union{Colon,AbstractRange,Integer}=Colon()) where {T,N}
+    A = N == 2 ? getindex(ga.A, i, j) : getindex(ga.A, i, j, k)
+    GeoArray(A, ga.f, ga.crs, ga.metadata)
+end
+
+Base.getindex(ga::GeoArray{T,3}, i::AbstractRange, j::AbstractRange) where {T} = Base.getindex(ga, i, j, :)
 
 Base.getindex(ga::GeoArray, I::Vararg{Integer,2}) = getindex(ga.A, I...)
 Base.getindex(ga::GeoArray, I::Vararg{Integer,3}) = getindex(ga.A, I...)
+function Base.getindex(ga::GeoArray{T,2}, I::Vararg{Integer,3}) where {T}
+    I[3] != 1 && throw(BoundsError(ga, I))
+    getindex(ga.A, I[1], I[2])
+end
 
 # Getindex and setindex! with floats
 """
@@ -146,10 +156,14 @@ function Base.getindex(ga::GeoArray, I::SVector{2,<:AbstractFloat})
 end
 Base.getindex(ga::GeoArray, I::Vararg{AbstractFloat,2}) = getindex(ga, SVector{2}(I))
 
-function Base.setindex!(ga::GeoArray, v, I::SVector{2,AbstractFloat})
-    i = indices(ga, I, Center())
+function Base.setindex!(ga::GeoArray{T,3}, v, I::SVector{2,AbstractFloat}) where {T}
+    i = indices(ga, I, Center(), RoundNearestTiesUp)
     ga.A[i, :] .= v
 end
+function Base.setindex!(ga::GeoArray{T,2}, v, I::SVector{2,AbstractFloat}) where {T}
+    i = indices(ga, I, Center(), RoundNearestTiesUp)
+    ga.A[i] = v
+end
 Base.setindex!(ga::GeoArray, v, I::Vararg{AbstractFloat,2}) = setindex!(ga, v, SVector{2}(I))
 Base.setindex!(ga::GeoArray, v, I::Vararg{Union{<:Integer,<:AbstractRange{<:Integer}},2}) = setindex!(ga.A, v, I...)
 Base.setindex!(ga::GeoArray, v, I::Vararg{Union{<:Integer,<:AbstractRange{<:Integer}},3}) = setindex!(ga.A, v, I...)
@@ -193,22 +207,23 @@ coords(ga::GeoArray, p::Tuple{<:Integer,<:Integer}, strategy::AbstractStrategy=C
 coords(ga::GeoArray, p::CartesianIndex{2}, strategy::AbstractStrategy=Center()) = coords(ga, SVector{2}(p.I), strategy)
 
 """
-    indices(ga::GeoArray, p::SVector{2,<:Real}, strategy::AbstractStrategy)
+    indices(ga::GeoArray, p::SVector{2,<:Real}, strategy::AbstractStrategy, rounding::RoundingMode)
 
 Retrieve logical indices of the cell represented by coordinates `p`.
 `strategy` can be used to define whether the coordinates represent the center (`Center`) or the top left corner (`Vertex`) of the cell.
+`rounding` can be used to define how the coordinates are rounded to the nearest integer index.
 See `coords` for the inverse function.
 """
-function indices(ga::GeoArray, p::SVector{2,<:Real}, strategy::AbstractStrategy)
-    CartesianIndex(Tuple(round.(Int, inv(ga.f)(p) .+ strategy.offset)))
+function indices(ga::GeoArray, p::SVector{2,<:Real}, strategy::AbstractStrategy, rounding::RoundingMode)
+    CartesianIndex(Tuple(round.(Int, inv(ga.f)(p) .+ strategy.offset, rounding)))
 end
-indices(ga::GeoArray, p::AbstractVector{<:Real}, strategy::AbstractStrategy=Center()) = indices(ga, SVector{2}(p), strategy)
-indices(ga::GeoArray, p::Tuple{<:Real,<:Real}, strategy::AbstractStrategy=Center()) = indices(ga, SVector{2}(p), strategy)
+indices(ga::GeoArray, p::AbstractVector{<:Real}, strategy::AbstractStrategy=Center(), rounding::RoundingMode=RoundNearestTiesUp) = indices(ga, SVector{2}(p), strategy, rounding)
+indices(ga::GeoArray, p::Tuple{<:Real,<:Real}, strategy::AbstractStrategy=Center(), rounding::RoundingMode=RoundNearestTiesUp) = indices(ga, SVector{2}(p), strategy, rounding)
 
 
 # Generate coordinates for complete GeoArray
 function coords(ga::GeoArray, strategy::AbstractStrategy=Center())
-    (ui, uj) = size(ga)[1:2]
+    (ui, uj) = size(ga)
     extra = typeof(strategy) == Center ? 0 : 1
     (coords(ga, SVector{2}(i, j), strategy) for i in 1:ui+extra, j in 1:uj+extra)
 end
diff --git a/src/geointerface.jl b/src/geointerface.jl
index 38a3334..c0d51b0 100644
--- a/src/geointerface.jl
+++ b/src/geointerface.jl
@@ -1,6 +1,2 @@
-function GI.extent(ga::GeoArray)
-    bb = bbox(ga)
-    Extents.Extent(X=(bb.min_x, bb.max_x), Y=(bb.min_y, bb.max_y))
-end
-
+GI.extent(ga::GeoArray) = bbox(ga)
 GI.crs(ga::GeoArray) = isempty(GFT.val(crs(ga))) ? nothing : crs(ga)
diff --git a/src/geoutils.jl b/src/geoutils.jl
index f74bf29..8f972ab 100644
--- a/src/geoutils.jl
+++ b/src/geoutils.jl
@@ -32,17 +32,17 @@ function is_rotated(ga::GeoArray)
 end
 
 function bbox(ga::GeoArray)
-    i, j, _ = size(ga)
+    i, j = size(ga)
     if is_rotated(ga)
         ax, ay = ga.f(SVector(0, 0))
         bx, by = ga.f(SVector(i, 0))
         cx, cy = ga.f(SVector(0, j))
         dx, dy = ga.f(SVector(i, j))
-        (min_x=min(ax, bx, cx, dx), min_y=min(ay, by, cy, dy), max_x=max(ax, bx, cx, dx), max_y=max(ay, by, cy, dy))
+        _convert(Extent, (min_x=min(ax, bx, cx, dx), min_y=min(ay, by, cy, dy), max_x=max(ax, bx, cx, dx), max_y=max(ay, by, cy, dy)))
     else
         ax, ay = ga.f(SVector(0, 0))
         bx, by = ga.f(SVector(i, j))
-        (min_x=min(ax, bx), min_y=min(ay, by), max_x=max(ax, bx), max_y=max(ay, by))
+        _convert(Extent, (min_x=min(ax, bx), min_y=min(ay, by), max_x=max(ax, bx), max_y=max(ay, by)))
     end
 end
 
@@ -67,17 +67,18 @@ function bbox!(ga::GeoArray, bbox::NamedTuple{(:min_x, :min_y, :max_x, :max_y)})
     ga.f = bbox_to_affine(size(ga)[1:2], bbox)
     ga
 end
+bbox!(ga::GeoArray, bbox::Extent{(:X, :Y)}) = bbox!(ga, _convert(NamedTuple, bbox))
 
 """Generate bounding boxes for GeoArray cells."""
 function bboxes(ga::GeoArray)
     c = collect(coords(ga, Vertex()))
     m, n = size(c)
-    cellbounds = Matrix{NamedTuple}(undef, (m - 1, n - 1))
+    cellbounds = Matrix{Extent{(:X, :Y)}}(undef, (m - 1, n - 1))
     for j in 1:n-1, i in 1:m-1
         v = c[i:i+1, j:j+1]
         minx, maxx = extrema(first.(v))
         miny, maxy = extrema(last.(v))
-        cellbounds[i, j] = (min_x=minx, max_x=maxx, min_y=miny, max_y=maxy)
+        cellbounds[i, j] = _convert(Extent, (; min_x=minx, max_x=maxx, min_y=miny, max_y=maxy))
     end
     cellbounds
 end
@@ -160,16 +161,16 @@ end
 
 Crop input GeoArray by coordinates (box or another GeoArray). If the coordinates range is larger than the GeoArray, only the overlapping part is given in the result.
 """
-function crop(ga::GeoArray, cbox::NamedTuple{(:min_x, :min_y, :max_x, :max_y)})
+function crop(ga::GeoArray{T,N}, cbox::Extent{(:X, :Y)}) where {T,N}
     # Check if ga and cbox overlap
-    if !bbox_overlap(bbox(ga), cbox)
-        error("GeoArray and crop box do not overlap")
+    if !intersects(bbox(ga), cbox)
+        error("GeoArray and crop extents do not overlap")
     end
 
     # Check extent and get bbox indices
-    ga_x, ga_y, = size(ga)
-    ii_min_x, ii_min_y = indices(ga, (cbox.min_x, cbox.min_y)).I
-    ii_max_x, ii_max_y = indices(ga, (cbox.max_x, cbox.max_y)).I
+    ga_x, ga_y = size(ga)
+    ii_min_x, ii_min_y = indices(ga, (cbox.X[1], cbox.Y[1]), Center(), RoundNearest).I
+    ii_max_x, ii_max_y = indices(ga, (cbox.X[2], cbox.Y[2]), Center(), RoundNearest).I
 
     # Determine indices for crop area
     i_min_x = max(min(ii_min_x, ii_max_x), 1)
@@ -178,8 +179,13 @@ function crop(ga::GeoArray, cbox::NamedTuple{(:min_x, :min_y, :max_x, :max_y)})
     i_max_y = min(max(ii_min_y, ii_max_y), ga_y)
 
     # Subset and return GeoArray
-    return ga[i_min_x:i_max_x, i_min_y:i_max_y, :]
+    if N == 2
+        return ga[i_min_x:i_max_x, i_min_y:i_max_y]
+    else
+        return ga[i_min_x:i_max_x, i_min_y:i_max_y, :]
+    end
 end
+crop(ga::GeoArray, cbox::NamedTuple{(:min_x, :min_y, :max_x, :max_y)}) = crop(ga, _convert(Extent, cbox))
 
 function crop(ga::GeoArray, cga::GeoArray)
     # Check if ga and crop ga CRS are the same
@@ -200,15 +206,15 @@ function sample!(ga::GeoArray, ga2::GeoArray)
     if ga.crs != ga2.crs
         error("GeoArrays have different CRS")
     end
-    wo, ho, zo = size(ga)
-    w, h = size(ga2)[1:2]
+    wo, ho, zo = _size(ga)
+    w, h = size(ga2)
 
     # Compose AffineMap that translates from logical coordinates
     # in `ga` to logical coordinates in `ga2`
-    x = inv(ga2.f) ∘ ga.f
+    f = inv(ga2.f) ∘ ga.f
 
     for io in 1:wo, jo in 1:ho
-        i, j = round.(Int, x((io, jo)) .+ 0.5)
+        i, j = round.(Int, f((io, jo)) .+ 0.5, RoundNearestTiesUp)
         if (1 <= i <= w) && (1 <= j <= h)
             # Loop over bands
             for z in 1:zo
@@ -220,7 +226,8 @@ function sample!(ga::GeoArray, ga2::GeoArray)
     ga
 end
 
-function _sizeof(ga::GeoArray, bbox)
+function _sizeof(ga::GeoArray, extent)
+    bbox = _convert(NamedTuple, extent)
     x = bbox.max_x - bbox.min_x
     y = bbox.max_y - bbox.min_y
     abs.(round.(Int, (x / ga.f.linear[1], y / ga.f.linear[4], size(ga)[3])))
@@ -274,19 +281,18 @@ function profile!(values, ga, a, b, band)
     δy = j1 - j0
 
     if abs(δx) >= abs(δy)
-        j = j0
         xstep = δx > 0 ? 1 : -1
         for (d, i) in enumerate(i0:xstep:i1)
-            push!(values, ga[i, j+div((d - 1) * δy, abs(δx), RoundNearest), band])
+            j = j0 + div((d - 1) * δy, abs(δx), RoundNearestTiesUp)
+            push!(values, ga[i, j, band])
         end
     else
-        i = i0
         ystep = δy > 0 ? 1 : -1
         for (d, j) in enumerate(j0:ystep:j1)
-            push!(values, ga[i+div((d - 1) * δx, abs(δy), RoundNearest), j, band])
+            i = i0 + div((d - 1) * δx, abs(δy), RoundNearest)
+            push!(values, ga[i, j, band])
         end
     end
-    return indices
 end
 
 # function Base.convert(
@@ -306,3 +312,15 @@ function Base.convert(
     am::CoordinateTransformations.AffineMap{<:AbstractMatrix{<:Real},<:AbstractVector{<:Real}})
     AffineMap(SMatrix{2,2,Float64}(am.linear), SVector{2,Float64}(am.translation))
 end
+
+function _convert(::Type{Extent}, nt::NamedTuple{(:min_x, :min_y, :max_x, :max_y)})
+    Extent(X=(nt.min_x, nt.max_x), Y=(nt.min_y, nt.max_y))
+end
+_convert(::Type{Extent}, nt::@NamedTuple{min_x::Float64, max_x::Float64, min_y::Float64, max_y::Float64}) = Extent(X=(nt.min_x, nt.max_x), Y=(nt.min_y, nt.max_y))
+
+function _convert(
+    ::Type{NamedTuple},
+    nt::Extent{(:X, :Y)},
+)
+    (; min_x=nt.X[1], min_y=nt.Y[1], max_x=nt.X[2], max_y=nt.Y[2])
+end
diff --git a/src/io.jl b/src/io.jl
index db42a08..c531b32 100644
--- a/src/io.jl
+++ b/src/io.jl
@@ -81,6 +81,10 @@ function GeoArray(dataset::ArchGDAL.RasterDataset, masked=true, band=nothing)
         dataset[mask] .= missing
     end
 
+    if size(dataset)[end] == 1 && masked
+        dataset = dropdims(dataset, dims=3)
+    end
+
     GeoArray(dataset, am, wkt, metadata)
 end
 
@@ -99,7 +103,7 @@ function write(fn::AbstractString, ga::GeoArray; nodata::Union{Nothing,Number}=n
     cancreate = ArchGDAL.metadataitem(driver, ArchGDAL.GDAL.GDAL_DCAP_CREATE) == "YES"
     cancopy = ArchGDAL.metadataitem(driver, ArchGDAL.GDAL.GDAL_DCAP_CREATECOPY) == "YES"
 
-    w, h, b = size(ga)
+    w, h, b = _size(ga)
     if isnothing(bandnames)
         bandnames = [nothing for i in 1:b]
     else
@@ -140,9 +144,8 @@ write(fn, ga, nodata=nothing, shortname=find_shortname(fn), options=Dict{String,
 
 
 function ArchGDAL.Dataset(ga::GeoArray, nodata=nothing, bandnames=nothing)
-    w, h, b = size(ga)
+    w, h, b = _size(ga)
 
-    w, h, b = size(ga)
     if isnothing(bandnames)
         bandnames = [nothing for i in 1:b]
     else
diff --git a/src/operations.jl b/src/operations.jl
index 33a7a14..3d285a9 100644
--- a/src/operations.jl
+++ b/src/operations.jl
@@ -1,26 +1,26 @@
 """Check whether two `GeoArrays`s `a` and `b` are
 geographically equal, although not necessarily in content."""
-function Base.isequal(a::GeoArray, b::GeoArray)
+function isgeoequal(a::GeoArray, b::GeoArray)
     size(a) == size(b) && a.f ≈ b.f && a.crs == b.crs
 end
 
 function Base.:-(a::GeoArray, b::GeoArray)
-    isequal(a, b) || throw(DimensionMismatch("Can't operate on non-geographic-equal `GeoArray`s"))
+    isgeoequal(a, b) || throw(DimensionMismatch("Can't operate on non-geographic-equal `GeoArray`s"))
     GeoArray(a.A .- b.A, a.f, a.crs, a.metadata)
 end
 
 function Base.:+(a::GeoArray, b::GeoArray)
-    isequal(a, b) || throw(DimensionMismatch("Can't operate on non-geographic-equal `GeoArray`s"))
+    isgeoequal(a, b) || throw(DimensionMismatch("Can't operate on non-geographic-equal `GeoArray`s"))
     GeoArray(a.A .+ b.A, a.f, a.crs, a.metadata)
 end
 
 function Base.:*(a::GeoArray, b::GeoArray)
-    isequal(a, b) || throw(DimensionMismatch("Can't operate on non-geographic-equal `GeoArray`s"))
+    isgeoequal(a, b) || throw(DimensionMismatch("Can't operate on non-geographic-equal `GeoArray`s"))
     GeoArray(a.A .* b.A, a.f, a.crs, a.metadata)
 end
 
 function Base.:/(a::GeoArray, b::GeoArray)
-    isequal(a, b) || throw(DimensionMismatch("Can't operate on non-geographic-equal `GeoArray`s"))
+    isgeoequal(a, b) || throw(DimensionMismatch("Can't operate on non-geographic-equal `GeoArray`s"))
     GeoArray(a.A ./ b.A, a.f, a.crs, a.metadata)
 end
 
@@ -43,7 +43,7 @@ julia> ga.A
  2  3
 ```
 """
-function Base.coalesce(ga::GeoArray{T,A}, v) where {T,A}
+function Base.coalesce(ga::GeoArray{T}, v) where {T}
     nT = nonmissingtype(T)
     cA = coalesce.(ga.A, nT(v))
     GeoArray(cA, ga.f, ga.crs, ga.metadata)
@@ -57,6 +57,11 @@ end
 Another GeoArray `like` can be passed as the second argument to `warp` to use the `like`'s `crs`, `extent` and `size` as the `ga` crs and resolution.
 The keyword `dest` is used to control where the temporary raster is stored. By default it is stored in memory, but can be set to a file path to directly save the warped GeoArray to disk.
 
+!!! warning
+    If no local PROJ data is available, (vertically) warping will silently fail. 
+    Use `enable_online_warp()` to enable (slow) network access to PROJ data.
+    For faster operations, use a utlity like `projsync` to download the data locally.
+
 # Examples
 ```julia-repl
 julia> ga = GeoArray(rand(100,100))
@@ -82,9 +87,20 @@ function warp(ga::GeoArray, like::GeoArray, options::Dict{String}=Dict{String,An
     warp(ga, noptions, dest)
 end
 
+"""
+    enable_online_warp(state::Bool=true)
+
+Enable or disable network access for PROJ data, required for `warp` if no local PROJ data is available.
+This has the same effect as setting the environement variable PROJ_NETWORK to "ON" *before* starting Julia.
+"""
+function enable_online_warp(state::Bool=true)
+    ArchGDAL.GDAL.osrsetprojenablenetwork(state)
+    return nothing
+end
+
 function warpoptions(ga::GeoArray)::Dict{String,Any}
     options = Dict{String,Any}(
-        "te" => values(GeoArrays.bbox(ga)),
+        "te" => values(_convert(NamedTuple, GeoArrays.bbox(ga))),
         "ts" => size(ga)[1:2],)
 
     srs = GFT.val(crs(ga))
diff --git a/test/test_geoarray.jl b/test/test_geoarray.jl
index 864a477..79ce1c7 100644
--- a/test/test_geoarray.jl
+++ b/test/test_geoarray.jl
@@ -38,10 +38,10 @@ end
 
 @testset "Reading rasters" begin
     ga = GeoArrays.read(joinpath(testdatadir, "data/utmsmall.tif"))
-    @test bbox(ga) == (min_x=440720.0, min_y=3.74532e6, max_x=446720.0, max_y=3.75132e6)
+    @test bbox(ga) == GeoArrays._convert(Extent, (; min_x=440720.0, min_y=3.74532e6, max_x=446720.0, max_y=3.75132e6))
     @inferred bbox(ga)
-    @test bboxes(ga)[1] == (min_x=440720.0, max_x=440780.0, min_y=3.75126e6, max_y=3.75132e6)
-    @test bboxes(ga)[end] == (min_x=446660.0, max_x=446720.0, min_y=3.74532e6, max_y=3.74538e6)
+    @test bboxes(ga)[1] == GeoArrays._convert(Extent, (; min_x=440720.0, max_x=440780.0, min_y=3.75126e6, max_y=3.75132e6))
+    @test bboxes(ga)[end] == GeoArrays._convert(Extent, (; min_x=446660.0, max_x=446720.0, min_y=3.74532e6, max_y=3.74538e6))
     @inferred bboxes(ga)
 end
 
@@ -93,7 +93,7 @@ end
 
 @testset "Conversion" begin
     ga = GeoArrays.GeoArray(rand(Int16, 10, 10))
-    gc = convert(GeoArrays.GeoArray{Float32}, ga)
+    gc = convert(GeoArrays.GeoArray{Float32,2}, ga)
     @test gc isa GeoArray
     @test eltype(gc) == Float32
     @test all(ga .== gc)
@@ -117,6 +117,8 @@ end
     @inferred indices(ga, X)
     @test ii == i
     @test jj == j
+
+    @test indices.(Ref(ga), GeoArrays.coords(ga)) == CartesianIndices(ga[:, :, 1])
 end
 
 @testset "Ranges" begin
diff --git a/test/test_geoutils.jl b/test/test_geoutils.jl
index 993b538..a617a33 100644
--- a/test/test_geoutils.jl
+++ b/test/test_geoutils.jl
@@ -1,9 +1,9 @@
 import GeoInterface as GI
 import GeoFormatTypes as GFT
-import Extents
+using Extents: Extents, Extent
 using CoordinateTransformations
 
-const tbbox = (min_x=440720.0, min_y=3.74532e6, max_x=446720.0, max_y=3.75132e6)
+const tbbox = GeoArrays._convert(Extent, (min_x=440720.0, min_y=3.74532e6, max_x=446720.0, max_y=3.75132e6))
 
 @testset "GeoUtils" begin
     @testset "bbox" begin
@@ -53,13 +53,13 @@ const tbbox = (min_x=440720.0, min_y=3.74532e6, max_x=446720.0, max_y=3.75132e6)
 
     @testset "crop by bbox or ga" begin
         ga = GeoArrays.read(joinpath(testdatadir, "data/utmsmall.tif"))
-        cbox = (min_x=441260.0, min_y=3.75012e6, max_x=442520.0, max_y=3.75078e6)
+        cbox = GeoArrays._convert(Extent, (min_x=441260.0, min_y=3.75012e6, max_x=442520.0, max_y=3.75078e6))
         ga_s = ga[10:30, 20:50, begin:end]
 
         ga_c = GeoArrays.crop(ga, cbox)
         ga_s_c = GeoArrays.crop(ga, ga_s)
 
-        @test size(ga_c) == (21, 11, 1)
+        @test GeoArrays._size(ga_c) == (21, 11, 1)
         @test GeoArrays.bbox(ga_c) == cbox
 
         @test size(ga_s_c) == size(ga_s)
@@ -80,6 +80,7 @@ const tbbox = (min_x=440720.0, min_y=3.74532e6, max_x=446720.0, max_y=3.75132e6)
         ga = GeoArrays.read(joinpath(testdatadir, "data/utmsmall.tif"))
         a = 441064.0, 3745555.0
         b = 442864.0, 3747309.0
+        box = bbox(ga)
         values = Float32[]
         GeoArrays.profile!(values, ga, a, b, 1)
         values_r = Float32[]
@@ -89,17 +90,18 @@ const tbbox = (min_x=440720.0, min_y=3.74532e6, max_x=446720.0, max_y=3.75132e6)
         # TODO profile with smaller dy than dx step
 
         struct LineString end
+        lcoord = [[1, 2], [3, 3.5]]
         GI.isgeometry(::LineString) = true
         GI.geomtrait(::LineString) = GI.LineStringTrait()
         GI.ncoord(::GI.LineStringTrait, ::LineString) = 2
         GI.ngeom(::GI.LineStringTrait, ::LineString) = 2
-        GI.getgeom(::GI.LineStringTrait, ::LineString) = [[1, 2], [3, 4]]
+        GI.getgeom(::GI.LineStringTrait, ::LineString) = lcoord
 
         ga = GeoArray(rand(4, 4))
         values = GeoArrays.profile(ga, LineString())
-        @test length(values) == 3
+        # @test length(values) == 3
 
-        GI.getgeom(::GI.LineStringTrait, ::LineString) = [[3, 4], [1, 2]]
+        GI.getgeom(::GI.LineStringTrait, ::LineString) = reverse(lcoord)
         values2 = GeoArrays.profile(ga, LineString())
         @test length(values2) == 3
         @test values == reverse(values2)
diff --git a/test/test_io.jl b/test/test_io.jl
index 1cc09d0..05af641 100644
--- a/test/test_io.jl
+++ b/test/test_io.jl
@@ -11,7 +11,18 @@ end
     for f in remotefiles
         @testset "Reading $f" begin
             ga = GeoArrays.read(joinpath(testdatadir, f), band=1)
-            @test last(size(ga)) == 1
+            @test last(GeoArrays._size(ga)) == 1
+        end
+    end
+end
+
+@testset "Coordinates/indices for rasters" begin
+    for f in remotefiles
+        @testset "Reading $f" begin
+            ga = GeoArrays.read(joinpath(testdatadir, f), band=1)
+            if !iszero(ga.f.linear)
+                @test indices.(Ref(ga), GeoArrays.coords(ga)) == CartesianIndices(ga)
+            end
         end
     end
 end
@@ -93,11 +104,11 @@ end
     end
     @testset "NetCDF" begin
         ga = GeoArrays.read("NetCDF:$(joinpath(testdatadir, "netcdf", "sentinel5p_fake.nc")):my_var")
-        @test size(ga) == (61, 89, 1)
+        @test GeoArrays._size(ga) == (61, 89, 1)
     end
     @testset "Virtual" begin
         ga = GeoArrays.read("/vsicurl/https://github.com/OSGeo/gdal/blob/master/autotest/alg/data/2by2.tif?raw=true")
-        @test size(ga) == (2, 2, 1)
+        @test GeoArrays._size(ga) == (2, 2, 1)
     end
     @testset "Complex numbers" begin
         for T in (Complex{Int16}, Complex{Int32}, Complex{Float32}, Complex{Float64})