Fix gdalwarp memory usage. (#172)

evetion · Dec 14, 2024 · 206d76b · 206d76b
1 parent 8fe698c
commit 206d76b
Show file tree

Hide file tree

Showing 4 changed files with 42 additions and 40 deletions.
diff --git a/CHANGELOG.md b/CHANGELOG.md
@@ -7,6 +7,12 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 ## [Unreleased]
 
+## [0.9.3] - 2024-12-14
+
+- Fixed increasing memory usage of gdalwarp over time.
+- Improved coding style, using more do blocks.
+- Fixed `nthreads` not being correctly set in `warp`.
+
 ## [0.9.2] - 2024-12-8
 
 ### Fixed

diff --git 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.9.2"
+version = "0.9.3"
 
 
 [deps]

diff --git a/src/io.jl b/src/io.jl
@@ -20,10 +20,9 @@ function read(fn::AbstractString; masked::Bool=true, band=nothing)
     return ga
 end
 
-function GeoArray(ds::ArchGDAL.Dataset)
-    dataset = ArchGDAL.RasterDataset(ds)
-    ga = GeoArray(dataset)
-    ArchGDAL.destroy(ds)
+function GeoArray(ds::ArchGDAL.AbstractDataset)
+    ga = GeoArray(ArchGDAL.RasterDataset(ds))
+    ArchGDAL.destroy(ds)  # masked is true
     return ga
 end
 
@@ -69,6 +68,7 @@ function GeoArray(dataset::ArchGDAL.RasterDataset, masked=true, band=nothing)
             else
                 @warn "Unknown/unsupported mask."
             end
+            ArchGDAL.destroy(band)
         end
     end
 
@@ -115,23 +115,23 @@ function write(fn::AbstractString, ga::GeoArray; nodata::Union{Nothing,Number}=n
 
         ArchGDAL.create(fn, driver=driver, width=w, height=h, nbands=b, dtype=dtype, options=stringlist(options)) do dataset
             for i = 1:b
-                band = ArchGDAL.getband(dataset, i)
-                ArchGDAL.write!(band, data[:, :, i])
-                !isnothing(nodata) && ArchGDAL.GDAL.gdalsetrasternodatavalue(band, nodata)
-                !isnothing(bandnames[i]) && ArchGDAL.GDAL.gdalsetdescription(band, bandnames[i])
+                ArchGDAL.getband(dataset, i) do band
+                    ArchGDAL.write!(band, data[:, :, i])
+                    !isnothing(nodata) && ArchGDAL.GDAL.gdalsetrasternodatavalue(band, nodata)
+                    !isnothing(bandnames[i]) && ArchGDAL.GDAL.gdalsetdescription(band, bandnames[i])
+                end
             end
 
             # Set geotransform and crs
             gt = affine_to_geotransform(ga.f)
             ArchGDAL.GDAL.gdalsetgeotransform(dataset, gt)
             ArchGDAL.GDAL.gdalsetprojection(dataset, GFT.val(ga.crs))
             setmetadata(dataset, ga.metadata)
-
         end
     elseif cancopy
-        dataset = ArchGDAL.Dataset(ga::GeoArray, nodata, bandnames)
-        ArchGDAL.copy(dataset, filename=fn, driver=driver, options=stringlist(options))
-        ArchGDAL.destroy(dataset)
+        ArchGDAL.Dataset(ga::GeoArray, nodata, bandnames) do dataset
+            ArchGDAL.copy(dataset, filename=fn, driver=driver, options=stringlist(options))
+        end
     else
         @error "Cannot create file with $shortname driver."
     end
@@ -143,7 +143,7 @@ write!(args...) = write(args...)
 write(fn, ga, nodata=nothing, shortname=find_shortname(fn), options=Dict{String,String}()) = write(fn, ga; nodata=nodata, shortname=shortname, options=options)
 
 
-function ArchGDAL.Dataset(ga::GeoArray, nodata=nothing, bandnames=nothing)
+function ArchGDAL.Dataset(f, ga::GeoArray, nodata=nothing, bandnames=nothing)
     w, h, b = _size(ga)
 
     if isnothing(bandnames)
@@ -154,25 +154,24 @@ function ArchGDAL.Dataset(ga::GeoArray, nodata=nothing, bandnames=nothing)
 
     data, dtype, nodata = prep(ga, nodata)
 
-    dataset = ArchGDAL.create(string("/vsimem/$(gensym())"), driver=ArchGDAL.getdriver("MEM"), width=w, height=h, nbands=b, dtype=dtype)
-    for i = 1:b
-        band = ArchGDAL.getband(dataset, i)
-        ArchGDAL.write!(band, data[:, :, i])
-        !isnothing(nodata) && ArchGDAL.GDAL.gdalsetrasternodatavalue(band, nodata)
-        !isnothing(bandnames[i]) && ArchGDAL.GDAL.gdalsetdescription(band, bandnames[i])
-    end
+    ArchGDAL.create(string("/vsimem/$(gensym())"), driver=ArchGDAL.getdriver("MEM"), width=w, height=h, nbands=b, dtype=dtype) do dataset
+        for i = 1:b
+            ArchGDAL.getband(dataset, i) do band
+                ArchGDAL.write!(band, data[:, :, i])
+                !isnothing(nodata) && ArchGDAL.GDAL.gdalsetrasternodatavalue(band, nodata)
+                !isnothing(bandnames[i]) && ArchGDAL.GDAL.gdalsetdescription(band, bandnames[i])
+            end
+        end
 
-    # Set geotransform and crs
-    gt = affine_to_geotransform(ga.f)
-    ArchGDAL.GDAL.gdalsetgeotransform(dataset, gt)
-    ArchGDAL.GDAL.gdalsetprojection(dataset, GFT.val(ga.crs))
-    setmetadata(dataset, ga.metadata)
-    return dataset
+        # Set geotransform and crs
+        gt = affine_to_geotransform(ga.f)
+        ArchGDAL.GDAL.gdalsetgeotransform(dataset, gt)
+        ArchGDAL.GDAL.gdalsetprojection(dataset, GFT.val(ga.crs))
+        setmetadata(dataset, ga.metadata)
+        f(dataset)
+    end
 end
 
-ArchGDAL.RasterDataset(ga::GeoArray) = ArchGDAL.RasterDataset(ArchGDAL.Dataset(ga))
-
-
 function prep(ga, nodata=nothing)
     need_nodata = false
     need_convert = false

diff --git a/src/operations.jl b/src/operations.jl
@@ -69,21 +69,19 @@ julia> epsg!(ga, 4326)
 julia> ga2 = GeoArrays.warp(ga, Dict("t_srs" => "EPSG:4326+3855"))
 ```
 """
-function warp(ga::GeoArray, options::Dict{String}, dest="/vsimem/$(gensym())")
-    dataset = ArchGDAL.Dataset(ga)
+function warp(ga::GeoArray, options::Dict{String}, dest="/vsimem/jlgeoarraystempwarp")
     noptions = Dict{String,Any}()
     warpdefaults!(noptions)
     merge!(noptions, options)
-    ds = [dataset]
-    nga = ArchGDAL.gdalwarp(ds, warpstringlist(options); dest
-    ) do warped
-        GeoArray(warped)
+    nga = ArchGDAL.Dataset(ga) do dataset
+        ArchGDAL.gdalwarp([dataset], warpstringlist(options); dest) do warped
+            GeoArray(warped)
+        end
     end
-    ArchGDAL.destroy(dataset)
     return nga
 end
 
-function warp(ga::GeoArray, like::GeoArray, options::Dict{String}=Dict{String,Any}(), dest="/vsimem/$(gensym())")
+function warp(ga::GeoArray, like::GeoArray, options::Dict{String}=Dict{String,Any}(), dest="/vsimem/jlgeoarraystempwarp")
     noptions = warpoptions(like)
     warpdefaults!(noptions)
     merge!(noptions, options)
@@ -112,6 +110,5 @@ function warpoptions(ga::GeoArray)::Dict{String,Any}
 end
 
 function warpdefaults!(d::Dict)
-    get!(d, "wo", Dict("NUM_THREADS" => string(Threads.nthreads)))
-    get!(d, "multi", "")
+    get!(d, "wo", Dict("NUM_THREADS" => string(Threads.nthreads())))
 end