Raster writing via ArchGDAL.copy suddenly broken?

[vlandau]

I'm noticing that raster writing is not working for my package, Omniscape.jl, when using ArchGDAL. I think it has something to do with a dependency of ArchGDAL being updated? I'm not 100% sure. It was working fine before, and it also seems to work when running from a top-level function in my package using a specific manifest.toml; however, when I try to run write_raster directly (see below, and also here), or use my package with a fresh install of ArchGDAL, it fails then too (even in the docker image with the specific manifest.toml).

I'm trying to write a tif raster (by copying an in-memory raster) using LZW compression. Instead of writing the values as it was before, it is writing an empty raster (all of it is "NoData")

Code to reproduce the issue is below. Anything that you notice right away that seems wrong?

Thanks very much for any assistance.

using Pkg; Pkg.add("ArchGDAL")
using ArchGDAL

function write_raster(fn_prefix::AbstractString,
                      array,
                      wkt::AbstractString,
                      transform,
                      file_format::String)
    # transponse array back to columns by rows
    array_t = permutedims(array, [2, 1])

    width, height = size(array_t)

    # Define extension and driver based in file_format
    file_format == "tif" ? (ext = ".tif"; driver = "GTiff") :
            (ext = ".asc"; driver = "AAIGrid")

    file_format == "tif" ? (options = ["COMPRESS=LZW"]) :
                           (options = [])

    # Append file extention to filename
    fn = string(fn_prefix, ext)

    # Create raster in memory
    ArchGDAL.create(fn_prefix,
                    driver = ArchGDAL.getdriver("MEM"),
                    width = width,
                    height = height,
                    nbands = 1,
                    dtype = eltype(array_t),
                    options = options) do dataset
        band = ArchGDAL.getband(dataset, 1)
        # Write data to band
        ArchGDAL.write!(band, array_t)

        # Write nodata and projection info
        ArchGDAL.setnodatavalue!(band, -9999.0)
        ArchGDAL.setgeotransform!(dataset, transform)
        ArchGDAL.setproj!(dataset, wkt)

        # Copy memory object to disk (necessary because ArchGDAL.create
        # does not support creation of ASCII rasters)
        ArchGDAL.copy(dataset,
                      filename = fn,
                      driver = ArchGDAL.getdriver(driver),
                      options = options)
    end
end

wkt = "PROJCS[\"NAD_1983_2011_StatePlane_Vermont_FIPS_4400\",GEOGCS[\"NAD83(2011)\"," *
    "DATUM[\"NAD83_National_Spatial_Reference_System_2011\",SPHEROID[\"GRS 1980\",6378137," *
    "298.257222101004,AUTHORITY[\"EPSG\",\"7019\"]],AUTHORITY[\"EPSG\",\"1116\"]]," *
    "PRIMEM[\"Greenwich\",0],UNIT[\"degree\",0.0174532925199433," *
    "AUTHORITY[\"EPSG\",\"9122\"]],AUTHORITY[\"EPSG\",\"6318\"]]," *
    "PROJECTION[\"Transverse_Mercator\"],PARAMETER[\"latitude_of_origin\",42.5]," *
    "PARAMETER[\"central_meridian\",-72.5],PARAMETER[\"scale_factor\",0.999964286]," *
    "PARAMETER[\"false_easting\",500000],PARAMETER[\"false_northing\",0]," *
    "UNIT[\"metre\",1,AUTHORITY[\"EPSG\",\"9001\"]],AXIS[\"Easting\",EAST]," *
    "AXIS[\"Northing\",NORTH],AUTHORITY[\"EPSG\",\"6589\"]]"

transform = [224746.09399999958, 800.0, 0.0, 556380.4152000006, 0.0, -800.0]
array = [    1.2     14.2    42.1    55  1.76   4.3  -9999.0  -9999.0  -9999.0  -9999.0;
             1.2     14.2    42.1    55  1.76   4.3  -9999.0  -9999.0  -9999.0  -9999.0;
         -9999.0  -9999.0    -1.2   1.2  14.2  42.1       55     1.76      4.3  -9999.0;
         -9999.0  -9999.0     4.5   1.2  14.2  42.1       55     1.76      4.3  -9999.0;
         -9999.0  -9999.0   4.124   1.2  14.2  42.1       55     1.76      4.3  -9999.0;
         -9999.0  -9999.0      15   1.2  14.2  42.1       55     1.76      4.3  -9999.0;
         -9999.0  -9999.0    33.2   1.2  14.2  42.1       55     1.76      4.3  -9999.0;
         -9999.0  -9999.0       0   1.2  14.2  42.1       55     1.76      4.3  -9999.0;
         -9999.0       42     1.2  14.2  42.1    55     1.76      4.3  -9999.0  -9999.0;
             1.2     14.2    42.1    55  1.76   4.3  -9999.0  -9999.0  -9999.0  -9999.0]

fn_prefix = "test_output"
file_format = "tif"

write_raster(fn_prefix,
             array,
             wkt,
             transform,
             file_format)

[vlandau]

I'll also note that writing as ASCII, (with file_format = "asc") works just fine still.

[visr]

Could this be the same issue as evetion/GeoArrays.jl#26? What versions of GDAL and it's dependencies are you using?

[vlandau]

It seems unrelated AFAICT. Here's the manifest.

[vlandau]

Were you able to reproduce the issue with the code I posted?

[visr]

Hmm yeah that manifest just has recent versions. And if I run your code I also get all missings.

But isn't the problem with how the raster is written? I don't have a fix now, but if I run it interactively

ds = write_raster(fn_prefix,
             array,
             wkt,
             transform,
             file_format)

AG.read(ds, 1)

I read out the same values that are put in. Is the problem perhaps that you don't close the dataset so the values are never flushed to disk? Returning an open dataset is probably not what you want write_raster to do.

[vlandau]

Interesting. I had assumed that the connection would be closed upon exiting the do block, similar to how:

open(file) do f
  # do some stuff
end

will automatically close the connection.

Do I need to add a line in write_raster to force that?

[vlandau]

But the open files issue would make sense given some other issues I've been having with trying to implement these functions in another related package. There are tons of tests, and I had been getting too many open files errors, so your hypothesis is in line with that.

[vlandau]

You had it right!! I needed to do ArchGDAL.create, then ArchGDAL.copy, then ArchGDAL.write!() so the data would be written to the file location on the disk.

My new function is below and it seems to be working. Thank you so much for the help! You're a life saver. As I said above, this will also likely resolve an issue that I was having in another package! Still seems strange that it was working okay with the AAIGrid driver? Either way, glad that it's working now!

function write_raster(fn_prefix::AbstractString,
                      array,
                      wkt::AbstractString,
                      transform,
                      file_format::String)
    # transponse array back to columns by rows
    array_t = permutedims(array, [2, 1])

    width, height = size(array_t)

    # Define extension and driver based in file_format
    file_format == "tif" ? (ext = ".tif"; driver = "GTiff") :
            (ext = ".asc"; driver = "AAIGrid")

    file_format == "tif" ? (options = ["COMPRESS=LZW"]) :
                           (options = [])

    # Append file extention to filename
    fn = string(fn_prefix, ext)

    # Create raster in memory
    ArchGDAL.create(fn_prefix,
                    driver = ArchGDAL.getdriver("MEM"),
                    width = width,
                    height = height,
                    nbands = 1,
                    dtype = eltype(array_t),
                    options = options) do ds

        # Copy memory object to disk (necessary because ArchGDAL.create
        # does not support creation of ASCII rasters)
        ArchGDAL.copy(ds,
                      filename = fn,
                      driver = ArchGDAL.getdriver(driver),
                      options = options) do dataset
            band = ArchGDAL.getband(dataset, 1)
            # Write data to band
            ArchGDAL.write!(band, array_t)

            # Write nodata and projection info
            ArchGDAL.setnodatavalue!(band, -9999.0)
            ArchGDAL.setgeotransform!(dataset, transform)
            ArchGDAL.setproj!(dataset, wkt)
        end
    end
end

[visr]

Ah great that it's working now. We definitely need more examples like this in the documentation.

I had assumed that the connection would be closed upon exiting the do block

That's right, but the copy at the end created a new one and was not created with a do block, yielding an IDataset that was not closed.

[vlandau]

Great! Thanks again for the help! Presumably as it is now all of the datasets should be closed properly?

[vlandau]

Hm, now it is working when I run it as in the code above, but when running it inside my package, I'm getting another error:
band 1: Write operation not permitted on dataset opened in read-only mode

Any ideas for a quick fix? Thanks again.

[visr]

Hmm are you sure that in your package it is not still dealing with an old read only opened dataset? What if you change the filename for instance?

[vlandau]

Yeah, the failure is actually occurring with Pkg.test. The error is triggered by write!. Here's the whole stack trace:

ERROR: LoadError: GDALError (CE_Failure, code 1):
    test3/cum_currmap.asc, band 1: Write operation not permitted on dataset opened in read-only mode

Stacktrace:
 [1] gdaljl_errorhandler(::GDAL.CPLErr, ::Int32, ::Cstring) at /root/.julia/packages/GDAL/cghbO/src/error.jl:36
 [2] rasterio! at /root/.julia/packages/ArchGDAL/JhHWz/src/raster/rasterio.jl:457 [inlined]
 [3] rasterio! at /root/.julia/packages/ArchGDAL/JhHWz/src/raster/rasterio.jl:455 [inlined]
 [4] rasterio!(::ArchGDAL.IRasterBand, ::Array{Float64,2}, ::GDAL.GDALRWFlag, ::Int64, ::Int64) at /root/.julia/packages/ArchGDAL/JhHWz/src/raster/rasterio.jl:103 (repeats 2 times)
 [5] write! at /root/.julia/packages/ArchGDAL/JhHWz/src/raster/rasterio.jl:175 [inlined]
 [6] (::Omniscape.var"#10#12"{String,Array{Float64,1},Array{Float64,2}})(::ArchGDAL.Dataset) at /home/omniscape/src/io.jl:72
 [7] copy(::Omniscape.var"#10#12"{String,Array{Float64,1},Array{Float64,2}}, ::ArchGDAL.Dataset; kwargs::Base.Iterators.Pairs{Symbol,Any,Tuple{Symbol,Symbol,Symbol},NamedTuple{(:filename, :driver, :options),Tuple{String,ArchGDAL.Driver,Array{Any,1}}}}) at /root/.julia/packages/ArchGDAL/JhHWz/src/context.jl:215
 [8] (::Omniscape.var"#9#11"{String,Array{Float64,1},Array{Float64,2},String})(::ArchGDAL.Dataset) at /home/omniscape/src/io.jl:66
 [9] create(::Omniscape.var"#9#11"{String,Array{Float64,1},Array{Float64,2},String}, ::String; kwargs::Base.Iterators.Pairs{Symbol,Any,NTuple{6,Symbol},NamedTuple{(:driver, :width, :height, :nbands, :dtype, :options),Tuple{ArchGDAL.Driver,Int64,Int64,Int64,DataType,Array{Any,1}}}}) at /root/.julia/packages/ArchGDAL/JhHWz/src/context.jl:215
 [10] write_raster(::String, ::Array{Float64,2}, ::String, ::Array{Float64,1}, ::String) at /home/omniscape/src/io.jl:56
 [11] run_omniscape(::String) at /home/omniscape/src/run_omniscape.jl:376
 [12] top-level scope at /home/omniscape/test/runtests.jl:86
 [13] include(::String) at ./client.jl:439
 [14] top-level scope at none:6

[visr]

Wait so for the GTiff it is fine to ArchGDAL.write!(band, array_t) to the dataset that copy creates, but not for AAIGrid? And what if you instead write the array to the MEM dataset like before? Or does the data then not end up in the file?

[vlandau]

It does appear to work with GTiff. Only the tests that attempt to write with AAIGrid trigger the error. I confirmed by adjusting all tests to write with GTiff, and the error wasn't triggered that time.

I think this might be the sweet spot? It seems to be working.

function write_raster(fn_prefix::AbstractString,
                      array,
                      wkt::AbstractString,
                      transform,
                      file_format::String)
    # transponse array back to columns by rows
    array_t = permutedims(array, [2, 1])

    width, height = size(array_t)

    # Define extension and driver based in file_format
    file_format == "tif" ? (ext = ".tif"; driver = "GTiff") :
            (ext = ".asc"; driver = "AAIGrid")

    file_format == "tif" ? (options = ["COMPRESS=LZW"]) :
                           (options = [])

    # Append file extention to filename
    fn = string(fn_prefix, ext)

    # Create raster in memory
    ArchGDAL.create(fn_prefix,
                    driver = ArchGDAL.getdriver("MEM"),
                    width = width,
                    height = height,
                    nbands = 1,
                    dtype = eltype(array_t),
                    options = options) do dataset
        band = ArchGDAL.getband(dataset, 1)
        # Write data to band
        ArchGDAL.write!(band, array_t)

        # Write nodata and projection info
        ArchGDAL.setnodatavalue!(band, -9999.0)
        ArchGDAL.setgeotransform!(dataset, transform)
        ArchGDAL.setproj!(dataset, wkt)
        # Copy memory object to disk (necessary because ArchGDAL.create
        # does not support creation of ASCII rasters)
        ArchGDAL.destroy(ArchGDAL.copy(dataset,
                      filename = fn,
                      driver = ArchGDAL.getdriver(driver),
                      options = options))
    end
    nothing
end

[vlandau]

Please feel free to copy any of the code above to use in examples if/as you see fit.

We definitely need more examples like this in the documentation.

[rafaqz]

@vlandau thanks for posting these, I'm setting up GeoData.jl to write all formats based on extension, and this is perfect.

[rafaqz]

BTW do you know if there is much overhead to using this for tif, as compared to just writing to disk?

[vlandau]

Not sure about overhead -- it seems pretty efficient for the most part, but I suppose there could be overhead from using GDAL as a middleman.