Tpoisot/issue87

Project.toml

@@ -1,7 +1,7 @@
 name = "SimpleSDMLayers"
 uuid = "2c645270-77db-11e9-22c3-0f302a89c64c"
 authors = ["Timothée Poisot <[email protected]>", "Gabriel Dansereau <[email protected]>"]
-version = "0.4.8"
+version = "0.4.9"
 
 [deps]
 ArchGDAL = "c9ce4bd3-c3d5-55b8-8973-c0e20141b8c3"

src/SimpleSDMLayers.jl

@@ -23,12 +23,14 @@ include(joinpath("datasets", "sources.jl"))
 include(joinpath("datasets", "download_layer.jl"))
 export EarthEnv, WorldClim, BioClim
 
-include(joinpath("datasets", "geotiff.jl"))
 include(joinpath("datasets", "raster.jl"))
+include(joinpath("datasets", "ascii.jl"))
+include(joinpath("datasets", "geotiff.jl"))
+export geotiff
+
 include(joinpath("datasets", "worldclim.jl"))
 include(joinpath("datasets", "chelsa.jl"))
 include(joinpath("datasets", "landcover.jl"))
-include(joinpath("datasets", "ascii.jl"))
 export worldclim, bioclim, landcover
 
 include(joinpath("operations", "coarsen.jl"))

src/datasets/ascii.jl

@@ -32,7 +32,7 @@ function ascii(file::AbstractString, datatype::Type{T}=Float64) where {T <: Numb
             grid[i] = nothing
         end
     end
-    return SimpleSDMPredictor(grid, xl, xl+cs*2ncols, yl, yl+cs*2nrows)
+    return SimpleSDMPredictor(grid, xl, xl+cs*ncols, yl, yl+cs*nrows)
 end
 
 """
@@ -47,7 +47,7 @@ function ascii(layer::SimpleSDMPredictor{T}, file::AbstractString; nodata::T=con
     open(file, "w") do io 
         write(io, "ncols $(size(layer, 2))\n")
         write(io, "nrows $(size(layer, 1))\n")
-        write(io, "cellsize $(stride(layer)[1])\n")
+        write(io, "cellsize $(2*stride(layer)[1])\n")
         write(io, "xllcorner $(layer.left)\n")
         write(io, "yllcorner $(layer.bottom)\n")
         write(io, "nodata_value $(nodata)\n")

src/datasets/geotiff.jl

@@ -76,3 +76,47 @@ function geotiff(
     return LT(buffer, left_pos-0.5lon_stride, right_pos+0.5lon_stride, bottom_pos-0.5lat_stride, top_pos+0.5lat_stride)
 
 end
+
+function geotiff(layer::SimpleSDMPredictor{T}, file::AbstractString; nodata::T=convert(T, -9999)) where {T <: Number}
+    array = layer.grid
+    replace!(array, nothing => NaN)
+    array = convert(Matrix{T}, array)
+    dtype = eltype(array)
+    array_t = reverse(permutedims(array, [2, 1]); dims=2)
+    width, height = size(array_t)
+
+    # Geotransform
+    gt = zeros(Float64, 6)
+    gt[1] = layer.left
+    gt[2] = 2stride(layer, 1)
+    gt[3] = 0.0
+    gt[4] = layer.top
+    gt[5] = 0.0
+    gt[6] = -2stride(layer, 2)
+
+    # Write
+    prefix = first(split(last(splitpath(file)), '.'))
+    ArchGDAL.create(prefix,
+                driver=ArchGDAL.getdriver("MEM"),
+                width=width, height=height,
+                nbands=1, dtype=T,
+                options=["COMPRESS=LZW"]) do dataset
+
+        band = ArchGDAL.getband(dataset, 1)
+
+        # Write data to band
+        ArchGDAL.write!(band, array_t)
+
+        # Write nodata and projection info
+        ArchGDAL.setnodatavalue!(band, nodata)
+        ArchGDAL.setgeotransform!(dataset, gt)
+        ArchGDAL.setproj!(dataset, "EPSG:4326")
+
+        # Write !
+        ArchGDAL.write(dataset, file, driver=ArchGDAL.getdriver("GTiff"), options=["COMPRESS=LZW"])
+    end
+end
+
+function geotiff(layer::SimpleSDMResponse{T}, file::AbstractString; nodata::T=convert(T, -9999)) where {T <: Number}
+    geotiff(convert(SimpleSDMPredictor, layer), file; nodata=nodata)
+end

test/ascii.jl

@@ -20,6 +20,10 @@ U = SimpleSDMLayers.ascii("test.asc")
 @test S.top == U.top
 @test size(S) == size(U)
 
+geotiff(U, "test.tiff")
+@test isfile("test.tiff")
+
 rm("test.asc")
+rm("test.tiff")
 
 end