diff --git a/src/sf/bbox.jl b/src/sf/bbox.jl index ea60262..424e1bd 100644 --- a/src/sf/bbox.jl +++ b/src/sf/bbox.jl @@ -24,6 +24,12 @@ function bbox2cellsize(b::bbox, size) end # 两个至少提供一个 +""" + bbox2dims(b::bbox; size=nothing, cellsize=nothing, reverse_lat=true) + +# Arguments +- `cellsize`: celly为负,则lat倒序 +""" function bbox2dims(b::bbox; size=nothing, cellsize=nothing, reverse_lat=true) if size !== nothing && cellsize === nothing cellsize = bbox2cellsize(b, size) diff --git a/src/sf/st_dims.jl b/src/sf/st_dims.jl index ef07a8d..b8929b7 100644 --- a/src/sf/st_dims.jl +++ b/src/sf/st_dims.jl @@ -48,6 +48,8 @@ function st_cellsize(lon::AbstractVector, lat::AbstractVector) lon[2] - lon[1], lat[2] - lat[1] # cellx, celly end +st_cellsize(f::AbstractString) = gdalinfo(f)["cellsize"] + # function st_cellsize(ra::AbstractRaster) # # x, y = st_dims(r) # # x[2] - x[1], y[2] - y[1] diff --git a/src/sf/st_mosaic.jl b/src/sf/st_mosaic.jl index 704469d..5bc6f91 100644 --- a/src/sf/st_mosaic.jl +++ b/src/sf/st_mosaic.jl @@ -32,7 +32,10 @@ function st_mosaic(rs::Vector{SpatRaster{T,N}}; kw...) where {T,N} end -function merge_var!(R::AbstractArray, f; var=nothing, box::bbox) +function merge_var!(R::AbstractArray, f; + var=nothing, box::bbox, cellsize=nothing) + + isnothing(cellsize) && (cellsize = st_cellsize(f)) b = st_bbox(f) bands = bandnames(f) nband = length(bands) @@ -40,24 +43,29 @@ function merge_var!(R::AbstractArray, f; var=nothing, box::bbox) ntime = length(inds_var) println("Reading data ...") - @time A = read_gdal(f, inds_var) + A = read_gdal(f, inds_var) ilon, ilat = bbox_overlap(b, box; cellsize) R[ilon, ilat, 1:ntime] .= A + nothing end function merge_var(fs; vars=nothing, var=nothing, + progress=true, box::bbox=bbox(-180, -60, 180, 90)) f = fs[1] - cellsize = gdalinfo(f)["cellsize"][1] + cellsize = st_cellsize(f) lon, lat = bbox2dims(box; cellsize) nlon, nlat = length(lon), length(lat) + bands = bandnames(f) ntime = isnothing(vars) ? length(bands) : length(bands) / length(vars) R = zeros(Float32, nlon, nlat, ntime) - @showprogress for f in fs - merge_var!(R, f; var, box) + p = Progress(ntime) + for f in fs + progress && next!(p) + merge_var!(R, f; var, box, cellsize) end R end diff --git a/test/sf/test-st_mosaic.jl b/test/sf/test-st_mosaic.jl index 3d4ea93..fa2cde6 100644 --- a/test/sf/test-st_mosaic.jl +++ b/test/sf/test-st_mosaic.jl @@ -1,10 +1,21 @@ -@testset "st_mosaic" begin +using Test, Ipaper, Ipaper.sf, ArchGDAL + +# box = st_bbox(st_bbox.(fs)) +@testset "st_mosaic" begin bands = string.(1:4) - r1 = rast(rand(4, 4, 4), bbox(-180.0, -60.0, 180.0, -30.0); bands) r2 = rast(rand(4, 4, 4), bbox(-180.0, -30.0, 180.0, 0.0); bands) - + r1 = rast(rand(4, 4, 4), bbox(-180.0, -60.0, 180.0, -30.0); bands) + rs = [r1, r2] r_big = st_mosaic(rs) - @test st_bbox(r_big) == bbox(-180.0, -60.0, 180.0, 0.0) - # r_big + @test st_bbox(r_big) == bbox(-180.0, -60.0, 180.0, 0.0) + + ## test for `merge_var` + fs = ["r1.tif", "r2.tif"] + write_gdal(r1, fs[1]) + write_gdal(r2, fs[2]) + + R = merge_var(fs; box) + @test cat(r2.A, r1.A, dims=2) ≈ R + rm.(fs) end