Permalink
Browse files

Add support for gdal utilities (#75)

* add support for gdal utilities

closes #74

* check for usage error

* comment out gdalwarp test that requires projection for now

* add test for gdalwarp

* type the default empty options vector to String

and add a test which uses the default.
It was already working fine, though this perhaps gives users a better hint

* fix copied gdaldem docstring

and consistenly use imperative
  • Loading branch information...
yeesian authored and visr committed Jan 30, 2019
1 parent fce734c commit 9cc11c3f3caa57b57680725a356b24db8ff069c8
Showing with 356 additions and 1 deletion.
  1. +1 −0 src/ArchGDAL.jl
  2. +3 −1 src/context.jl
  3. +234 −0 src/utilities.jl
  4. +1 −0 test/runtests.jl
  5. +117 −0 test/test_gdalutilities.jl
@@ -20,6 +20,7 @@ module ArchGDAL
include("ogr/featuredefn.jl")
include("ogr/fielddefn.jl")
include("ogr/styletable.jl")
include("utilities.jl")
include("context.jl")
include("base/iterators.jl")
include("base/display.jl")
@@ -69,7 +69,9 @@ for gdalfunc in (
:createmultipolygon, :createmultipolygon_noholes, :createpoint,
:createpolygon, :createRAT, :createstylemanager, :createstyletable,
:createstyletool, :delaunaytriangulation, :difference, :forceto,
:fromGML, :fromJSON, :fromWKB, :fromWKT, :getcurvegeom, :getfeature,
:fromGML, :fromJSON, :fromWKB, :fromWKT, :gdalbuildvrt, :gdaldem,
:gdalgrid, :gdalnearblack, :gdalrasterize, :gdaltranslate,
:gdalvectortranslate, :gdalwarp, :getcurvegeom, :getfeature,
:getlineargeom, :getpart, :intersection, :importEPSG, :importEPSGA,
:importESRI, :importPROJ4, :importWKT, :importXML, :importURL,
:newspatialref, :nextfeature, :pointalongline, :pointonsurface,
@@ -0,0 +1,234 @@
"Throw an error if the pointer to int variable indicates an error."
function _failsafe(usage_error::Ref{Cint})
# follows https://github.com/JuliaGeo/GDAL.jl/blob/017bf6b8492dcd2186ced297076b283c3591d798/src/error.jl#L31-L37
if usage_error[] != 0 && GDAL.getlasterrortype() in GDAL.throw_class
throw(GDAL.GDALError())
end
end

"""
List various information about a GDAL supported raster dataset.
### Parameters
* **dataset**: The source dataset.
* **options**: List of options (potentially including filename and open
options). The accepted options are the ones of the gdalinfo utility.
### Returns
String corresponding to the information about the raster dataset.
"""
function gdalinfo(dataset::Dataset, options = String[])
options = GDAL.infooptionsnew(options, C_NULL)
result = GDAL.info(dataset.ptr, options)
GDAL.infooptionsfree(options)
return result
end

"""
Convert raster data between different formats.
### Parameters
* **dataset**: The dataset to be translated.
* **options**: List of options (potentially including filename and open
options). The accepted options are the ones of the gdal_translate utility.
### Returns
The output dataset.
"""
function unsafe_gdaltranslate(
dataset::Dataset,
options = String[];
dest = "/vsimem/tmp"
)
options = GDAL.translateoptionsnew(options, C_NULL)
usage_error = Ref{Cint}()
result = GDAL.translate(dest, dataset.ptr, options, usage_error)
GDAL.translateoptionsfree(options)
_failsafe(usage_error)
return Dataset(result)
end

"""
Image reprojection and warping function.
### Parameters
* **datasets**: The list of input datasets.
* **options**: List of options (potentially including filename and open
options). The accepted options are the ones of the gdalwarp utility.
### Returns
The output dataset.
"""
function unsafe_gdalwarp(
datasets::Vector{Dataset},
options = String[];
dest = "/vsimem/tmp"
)
options = GDAL.warpappoptionsnew(options, C_NULL)
usage_error = Ref{Cint}()
result = GDAL.warp(dest, Ptr{GDAL.GDALDatasetH}(C_NULL), length(datasets),
[ds.ptr for ds in datasets], options, usage_error)
GDAL.warpappoptionsfree(options)
_failsafe(usage_error)
return Dataset(result)
end

"""
Convert vector data between file formats.
### Parameters
* **datasets**: The list of input datasets (only 1 supported currently).
* **options**: List of options (potentially including filename and open
options). The accepted options are the ones of the ogr2ogr utility.
### Returns
The output dataset.
"""
function unsafe_gdalvectortranslate(
datasets::Vector{Dataset},
options = String[];
dest = "/vsimem/tmp"
)
options = GDAL.vectortranslateoptionsnew(options, C_NULL)
usage_error = Ref{Cint}()
result = GDAL.vectortranslate(dest, Ptr{GDAL.GDALDatasetH}(C_NULL),
length(datasets), [ds.ptr for ds in datasets], options, usage_error)
GDAL.vectortranslateoptionsfree(options)
_failsafe(usage_error)
return Dataset(result)
end

"""
Tools to analyze and visualize DEMs.
### Parameters
* **dataset**: The source dataset.
* **pszProcessing**: the processing to apply (one of "hillshade", "slope",
"aspect", "color-relief", "TRI", "TPI", "Roughness").
* **options**: List of options (potentially including filename and open
options). The accepted options are the ones of the gdaldem utility.
# Keyword Arguments
* **colorfile**: color file (mandatory for "color-relief" processing,
should be empty otherwise).
### Returns
The output dataset.
"""
function unsafe_gdaldem(
dataset::Dataset,
processing::String,
options = String[];
dest = "/vsimem/tmp",
colorfile = C_NULL
)
if processing == "color-relief"
@assert colorfile != C_NULL
end
options = GDAL.demprocessingoptionsnew(options, C_NULL)
usage_error = Ref{Cint}()
result = GDAL.demprocessing(dest, dataset.ptr, processing, colorfile,
options, usage_error)
GDAL.demprocessingoptionsfree(options)
_failsafe(usage_error)
return Dataset(result)
end

"""
Convert nearly black/white borders to exact value.
### Parameters
* **dataset**: The source dataset.
* **options**: List of options (potentially including filename and open
options). The accepted options are the ones of the nearblack utility.
### Returns
The output dataset.
"""
function unsafe_gdalnearblack(
dataset::Dataset,
options = String[];
dest = "/vsimem/tmp"
)
options = GDAL.nearblackoptionsnew(options, C_NULL)
usage_error = Ref{Cint}()
result = GDAL.nearblack(dest, Ptr{GDAL.GDALDatasetH}(C_NULL), dataset.ptr,
options, usage_error)
GDAL.nearblackoptionsfree(options)
_failsafe(usage_error)
return Dataset(result)
end

"""
Create a raster from the scattered data.
### Parameters
* **dataset**: The source dataset.
* **options**: List of options (potentially including filename and open
options). The accepted options are the ones of the gdal_grid utility.
### Returns
The output dataset.
"""
function unsafe_gdalgrid(
dataset::Dataset,
options = String[];
dest = "/vsimem/tmp"
)
options = GDAL.gridoptionsnew(options, C_NULL)
usage_error = Ref{Cint}()
result = GDAL.grid(dest, dataset.ptr, options, usage_error)
GDAL.gridoptionsfree(options)
_failsafe(usage_error)
return Dataset(result)
end

"""
Burn vector geometries into a raster.
### Parameters
* **dataset**: The source dataset.
* **options**: List of options (potentially including filename and open
options). The accepted options are the ones of the gdal_rasterize utility.
### Returns
The output dataset.
"""
function unsafe_gdalrasterize(
dataset::Dataset,
options = String[];
dest = "/vsimem/tmp"
)
options = GDAL.rasterizeoptionsnew(options, C_NULL)
usage_error = Ref{Cint}()
result = GDAL.rasterize(dest, Ptr{GDAL.GDALDatasetH}(C_NULL), dataset.ptr,
options, usage_error)
GDAL.rasterizeoptionsfree(options)
_failsafe(usage_error)
return Dataset(result)
end

"""
Build a VRT from a list of datasets.
### Parameters
* **datasets**: The list of input datasets.
* **options**: List of options (potentially including filename and open
options). The accepted options are the ones of the gdalbuildvrt utility.
### Returns
The output dataset.
"""
function unsafe_gdalbuildvrt(
datasets::Vector{Dataset},
options = String[];
dest = "/vsimem/tmp"
)
options = GDAL.buildvrtoptionsnew(options, C_NULL)
usage_error = Ref{Cint}()
result = GDAL.buildvrt(dest, length(datasets), [ds.ptr for ds in datasets],
C_NULL, options, usage_error)
GDAL.buildvrtoptionsfree(options)
_failsafe(usage_error)
return Dataset(result)
end
@@ -80,6 +80,7 @@ end
include("test_dataset.jl")
include("test_rasterband.jl")
include("test_rasterio.jl")
include("test_gdalutilities.jl")
include("test_rasterattrtable.jl")
include("test_ospy_examples.jl")
include("test_geos_operations.jl")
@@ -0,0 +1,117 @@
using ArchGDAL, GDAL; AG = ArchGDAL
using Test

AG.registerdrivers() do
AG.read("data/utmsmall.tif") do ds_small
@testset "GDAL Error" begin
@test_throws GDAL.GDALError AG.gdalinfo(ds_small, ["-novalidoption"])
@test_throws GDAL.GDALError AG.unsafe_gdaltranslate(ds_small, ["-novalidoption"])
@test_throws GDAL.GDALError AG.unsafe_gdalbuildvrt([ds_small], ["-novalidoption"])
@test_throws GDAL.GDALError AG.unsafe_gdaldem(ds_small, "hillshade", ["-novalidoption"])
@test_throws GDAL.GDALError AG.unsafe_gdalnearblack(ds_small, ["-novalidoption"])
@test_throws GDAL.GDALError AG.unsafe_gdalwarp([ds_small], ["-novalidoption"])
end

@testset "GDAL Info" begin
infostr = AG.gdalinfo(ds_small, ["-checksum"])
@test occursin("Checksum=50054", infostr)
info_default = AG.gdalinfo(ds_small)
@test occursin("Driver: GTiff/GeoTIFF", info_default)
end

AG.gdaltranslate(ds_small, # resample to a 5×5 ascii grid
["-of","AAIGrid","-r","cubic","-tr","1200","1200"]
) do ds_tiny
@testset "GDAL Translate" begin
@test AG.read(ds_tiny, 1) == [128 171 127 93 83;
126 164 148 114 101;
161 175 177 164 140;
185 206 205 172 128;
193 205 209 181 122]
end

@testset "GDAL Build VRT" begin
AG.gdalbuildvrt([ds_tiny]) do ds_vrt
@test AG.read(ds_vrt, 1) == [128 171 127 93 83;
126 164 148 114 101;
161 175 177 164 140;
185 206 205 172 128;
193 205 209 181 122]
end
end

@testset "GDAL DEM Processing" begin
AG.gdaldem(ds_tiny, "hillshade", ["-of","AAIGrid"]) do ds_dempr
@test AG.read(ds_dempr, 1) == [ 0 0 0 0 0;
0 183 180 181 0;
0 184 182 181 0;
0 183 181 177 0;
0 0 0 0 0]
end
end

@testset "GDAL Near Black" begin
AG.gdalnearblack(ds_tiny, ["-of","GTiff","-color","0"]) do ds_nearblack
@test AG.read(ds_nearblack, 1) == [ 0 0 0 0 0;
0 0 0 0 0;
0 0 177 0 0;
0 0 0 0 0;
0 0 0 0 0]
end
end
end

# cannot reproject file on AppVeyor yet
# GDALError (CE_Failure, code 4):
# Unable to open EPSG support file gcs.csv. Try setting the
# GDAL_DATA environment variable to point to the directory
# containing EPSG csv files.
# @testset "GDAL Warp" begin
# AG.gdalwarp([ds_small], ["-of","MEM","-t_srs","EPSG:4326"]) do ds_warped
# @test AG.width(ds_small) == 100
# @test AG.height(ds_small) == 100
# @test AG.width(ds_warped) == 109
# @test AG.height(ds_warped) == 91
# end
# end
@testset "GDAL Warp" begin
AG.gdalwarp([ds_small], ["-of","MEM"]) do ds_warped
@test AG.width(ds_small) == 100
@test AG.height(ds_small) == 100
@test AG.shortname(AG.getdriver(ds_small)) == "GTiff"
@test AG.width(ds_warped) == 100
@test AG.height(ds_warped) == 100
@test AG.shortname(AG.getdriver(ds_warped)) == "MEM"
end
end
end

AG.read("data/point.geojson") do ds_point
@testset "GDAL Grid" begin
AG.gdalgrid(ds_point, ["-of","MEM","-outsize","3",
"10","-txe","100","100.3","-tye","0","0.1"]) do ds_grid
@test AG.getgeotransform(ds_grid) ≈ [100.0,0.1,0.0,0.0,0.0,0.01]
end
end

@testset "GDAL Rasterize" begin
AG.gdalrasterize(ds_point, ["-of","MEM","-tr","0.05","0.05"]) do ds_rasterize
@test AG.getgeotransform(ds_rasterize) ≈ [99.975,0.05,0.0,0.1143,0.0,-0.05]
end
end

@testset "GDAL Vector Translate" begin
AG.gdalvectortranslate([ds_point], ["-f","CSV","-lco",
"GEOMETRY=AS_XY"], dest = "data/point.csv") do ds_csv
end
@test replace(read("data/point.csv", String), "\r" => "") == """
X,Y,FID,pointname
100,0,2,point-a
100.2785,0.0893,3,point-b
100,0,0,a
100.2785,0.0893,3,b
"""
rm("data/point.csv")
end
end
end

0 comments on commit 9cc11c3

Please sign in to comment.