Add possibility to save as .tif with compression

[glaroc]

We are producing a very large number of Circuitscape outputs and the .asc format is really not a good format for storing large rasters. It seems like it would be quite straightforward to add the possibility to save the outputs as .tif with LZW lossless compression using GDAL.jl. This would create files that are much smaller and more optimized for GIS viewing.

[ViralBShah]

Is GDAL now available through BB?

[ViralBShah]

A PR would be much appreciated, since we otherwise have a long list of things to address in Circuitscape and may not get around to this before JuliaCon.

[ranjanan]

GDAL is available as an artifact: https://github.com/JuliaBinaryWrappers/GDAL_jll.jl

[glaroc]

I'll try to do a PR, but I can't guarantee I'll have time soon. I'll try!

[ranjanan]

@glaroc I'm happy to help along the way. Feel free to ping me here or on Slack if you have any questions.

[vlandau]

If this is figured out, might make it low-hanging fruit to also allow tif files as input, which would help make workflows more efficient. +1 from me!

[glaroc]

I think the easiest method would be to use the GeoArrays package to load the tif files, then copy the CRS from the input to create GeoArrays for the outputs and save them as tifs again with GeoArrays.write(). This workflow should be quite straightforward.

[glaroc]

Note that this should also speed up i/o quite a bit, especially for large files. I just did a test by importing the same file as asc and geotiff:

@time ga2=readdlm("test.asc",Float32,skipstart=6)
55.028136 seconds (491.58 M allocations: 13.272 GiB, 11.95% gc time)
9524×13197 Array{Float32,2}:

@time ga = GeoArrays.read("test.tif")
3.126899 seconds (96 allocations: 988.900 MiB)
13197×9524×1 GeoArray{Float32}:

[vlandau]

Is there an easy way of converting ga to a 2D Float64 array? If so it would be as easy as doing the file read as you have it, then doing a conversion, and leaving everything else in Circuitscape as-is. All of the functions in Circuitscape (and likely those from many of its dependencies) don't have methods for GeoArrays.

[glaroc]

I'm not a Julia expert, but it should be as simple as:

@time ga2=Array{Float32}(ga)[:,:,1]'
25.763421 seconds (628.45 M allocations: 12.174 GiB, 15.61% gc time)
9524×13197 LinearAlgebra.Adjoint{Float32,Array{Float32,2}}:

There might be a better, more efficient way of doing it...?

[ViralBShah]

That seems to kill a lot of the speed. Any idea why the conversion is so slow? Also may need a copy or something, since this is still an adjoint. I wonder if permutedims is faster.

[glaroc]

I'm not sure why this is so slow. I guess it's just going from a 2x2xN array to a 2x2 array?

A simple solution is just to extract the first band. permutedims seems to work.

@time ga2=permutedims(ga[:,:,1])
22.075344 seconds (491.36 M allocations: 8.266 GiB, 13.30% gc time)
9524×13197 Array{Float32,2}

[vlandau]

I'm 100% for adding support for tifs, so this is just playing devil's advocate...

One potential downside to accommodating tifs is that it requires users to have GDAL, which is a rather large dependency. Does Package.add("GeoArrays") also add GDAL to the host, or does it need to be installed manually? A manual install might prove a challenge for some users (but it could be easy enough to just link to instructions from the Circuitscape docs).

Many Circuitscape users are already "spatial" people, so many would probably already have GDAL installed for other applications (such as GIS packages in R). The added dependency might not really be a big deal at all -- and users probably won't mind the slight hassle of installing GDAL and accommodating the space it takes up if it means they no longer need to worry about file format conversion.

[vlandau]

A GDAL dep could also open up a whole world of additional possibilities for new Circuitscape features.

[glaroc]

@vlandau Since GDAL is available as a Binary Build (see BB discussion above), GDAL doesn't need to be compiled and installed separately. Just need to run

Pkg.add("GDAL_jll")

[vlandau]

Awesome. I wasn't sure what a BB was -- I'm still a bit of a programming newbie! :)

[ViralBShah]

Yes, it will all come automatically. We definitely don't want to use a GDAL on people's computers. With BB (BinaryBuilder), we provide precompiled binaries for all architectures. See: https://julialang.org/blog/2019/11/artifacts/

In fact GDAL_jll will just be a Circuitscape depenency. So people will get it automatically.

[ViralBShah]

Awesome. I wasn't sure what a BB was -- I'm still a bit of a programming newbie! :)

BB is a very Julia specific acronym and a very recent one that only package authors mostly know about! We're going to get a website and more information on it out eventually...

[glaroc]

Just published a PR that implements this using GeoArrays. Seems to work well with my files, but I didn't do in-depth testing.

Edit: I just cleaned up the PR and resubmitted it.

#227