From b13ded9383c1f6c52a60d9a8ee94bdf6516e442d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Timoth=C3=A9e=20Poisot?= Date: Tue, 23 Mar 2021 19:33:20 -0400 Subject: [PATCH 1/3] ASCII grid size is two small Fixes #87 --- Project.toml | 2 +- src/datasets/ascii.jl | 4 ++-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/Project.toml b/Project.toml index 26d1b470..083fa7f9 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "SimpleSDMLayers" uuid = "2c645270-77db-11e9-22c3-0f302a89c64c" authors = ["Timothée Poisot ", "Gabriel Dansereau "] -version = "0.4.8" +version = "0.4.9" [deps] ArchGDAL = "c9ce4bd3-c3d5-55b8-8973-c0e20141b8c3" diff --git a/src/datasets/ascii.jl b/src/datasets/ascii.jl index 917e59be..c792f229 100644 --- a/src/datasets/ascii.jl +++ b/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") From 680867639b05486fd7649493947252c1cb4c752e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Timoth=C3=A9e=20Poisot?= Date: Tue, 23 Mar 2021 20:03:39 -0400 Subject: [PATCH 2/3] =?UTF-8?q?=F0=9F=8E=89=20write=20geotiff?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- src/SimpleSDMLayers.jl | 6 ++++-- src/datasets/geotiff.jl | 39 +++++++++++++++++++++++++++++++++++++++ 2 files changed, 43 insertions(+), 2 deletions(-) diff --git a/src/SimpleSDMLayers.jl b/src/SimpleSDMLayers.jl index 0c2d1ae2..32e2ac2d 100644 --- a/src/SimpleSDMLayers.jl +++ b/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")) diff --git a/src/datasets/geotiff.jl b/src/datasets/geotiff.jl index 8af0056d..0c7af2ff 100644 --- a/src/datasets/geotiff.jl +++ b/src/datasets/geotiff.jl @@ -76,3 +76,42 @@ 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] = bc.left + gt[2] = 2stride(bc, 1) + gt[3] = 0.0 + gt[4] = bc.top + gt[5] = 0.0 + gt[6] = -2stride(bc, 2) + + # Write + ArchGDAL.create(fn_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, -9999) + ArchGDAL.setgeotransform!(dataset, gt) + ArchGDAL.setproj!(dataset, "EPSG:4326") + + # Write ! + ArchGDAL.write(dataset, file, driver=ArchGDAL.getdriver(GTiff), options=["COMPRESS=LZW"]) + end +end \ No newline at end of file From 971c2d5c931a8dc24aaa4697af70d8c96c1c09fd Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Timoth=C3=A9e=20Poisot?= Date: Tue, 23 Mar 2021 20:11:09 -0400 Subject: [PATCH 3/3] =?UTF-8?q?=E2=9C=85=20tests=20for=20tiff?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- src/datasets/geotiff.jl | 21 +++++++++++++-------- test/ascii.jl | 4 ++++ 2 files changed, 17 insertions(+), 8 deletions(-) diff --git a/src/datasets/geotiff.jl b/src/datasets/geotiff.jl index 0c7af2ff..058e3156 100644 --- a/src/datasets/geotiff.jl +++ b/src/datasets/geotiff.jl @@ -77,7 +77,7 @@ function geotiff( end -function geotiff(layer::SimpleSDMPredictor{T}, file::AbstractString; nodata::T=convert(T, -9999)) where {T <: Number}) +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) @@ -87,15 +87,16 @@ function geotiff(layer::SimpleSDMPredictor{T}, file::AbstractString; nodata::T=c # Geotransform gt = zeros(Float64, 6) - gt[1] = bc.left - gt[2] = 2stride(bc, 1) + gt[1] = layer.left + gt[2] = 2stride(layer, 1) gt[3] = 0.0 - gt[4] = bc.top + gt[4] = layer.top gt[5] = 0.0 - gt[6] = -2stride(bc, 2) + gt[6] = -2stride(layer, 2) # Write - ArchGDAL.create(fn_prefix, + prefix = first(split(last(splitpath(file)), '.')) + ArchGDAL.create(prefix, driver=ArchGDAL.getdriver("MEM"), width=width, height=height, nbands=1, dtype=T, @@ -107,11 +108,15 @@ function geotiff(layer::SimpleSDMPredictor{T}, file::AbstractString; nodata::T=c ArchGDAL.write!(band, array_t) # Write nodata and projection info - ArchGDAL.setnodatavalue!(band, -9999) + 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"]) + 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 \ No newline at end of file diff --git a/test/ascii.jl b/test/ascii.jl index 2442bbff..7787dffa 100644 --- a/test/ascii.jl +++ b/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 \ No newline at end of file