Skip to content
Permalink

Comparing changes

This is a direct comparison between two commits made in this repository or its related repositories. View the default comparison for this range or learn more about diff comparisons.

Open a pull request

Create a new pull request by comparing changes across two branches. If you need to, you can also . Learn more about diff comparisons here.
base repository: yeesian/ArchGDAL.jl
Failed to load repositories. Confirm that selected base ref is valid, then try again.
Loading
base: b113963532af1dc902dce04a1b877c79fcaaf6e8
Choose a base ref
..
head repository: yeesian/ArchGDAL.jl
Failed to load repositories. Confirm that selected head ref is valid, then try again.
Loading
compare: d860b69636e52df3f2586b694286bb927177617d
Choose a head ref
Showing with 65 additions and 15 deletions.
  1. +1 −0 docs/Project.toml
  2. +9 −0 docs/src/projections.md
  3. +1 −1 src/dataset.jl
  4. +12 −6 src/geointerface.jl
  5. +5 −5 src/ogr/geometry.jl
  6. +3 −1 src/spatialref.jl
  7. +9 −1 test/test_geometry.jl
  8. +25 −1 test/test_geos_operations.jl
1 change: 1 addition & 0 deletions docs/Project.toml
Original file line number Diff line number Diff line change
@@ -4,6 +4,7 @@ DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
DiskArrays = "3c3547ce-8d99-4f5e-a174-61eb10b00ae3"
Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
GDAL = "add2ef01-049f-52c4-9ee2-e494f65e021a"
GeoFormatTypes = "68eda718-8dee-11e9-39e7-89f7f65f511f"
Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"

[compat]
9 changes: 9 additions & 0 deletions docs/src/projections.md
Original file line number Diff line number Diff line change
@@ -62,6 +62,15 @@ ArchGDAL.createcoordtrans(source, target) do transform
end
```

```@example projections
using DataFrames
import GeoFormatTypes as GFT
coords = zip(rand(10), rand(10))
df = DataFrame(geom=ArchGDAL.createpoint.(coords), name="test");
df.geom = ArchGDAL.reproject(df.geom, GFT.EPSG(4326), GFT.EPSG(28992))
```

## Reprojecting from a layer
```@setup projections
# Getting vector data
2 changes: 1 addition & 1 deletion src/dataset.jl
Original file line number Diff line number Diff line change
@@ -800,7 +800,7 @@ function unsafe_executesql(
dataset::AbstractDataset,
query::AbstractString;
dialect::AbstractString = "",
spatialfilter::Geometry = Geometry(C_NULL),
spatialfilter::AbstractGeometry = Geometry(C_NULL),
)::FeatureLayer
return FeatureLayer(
GDAL.gdaldatasetexecutesql(dataset, query, spatialfilter, dialect),
18 changes: 12 additions & 6 deletions src/geointerface.jl
Original file line number Diff line number Diff line change
@@ -1,5 +1,7 @@
using Extents

geointerface_geomtype(::GeoInterface.AbstractGeometryTrait) = IGeometry

const lookup_method = Dict{DataType,Function}(
GeoInterface.PointTrait => createpoint,
GeoInterface.MultiPointTrait => createmultipoint,
@@ -280,16 +282,20 @@ let pointtypes = (wkbPoint, wkbPoint25D, wkbPointM, wkbPointZM),

function GeoInterface.convert(
::Type{T},
type::GeometryTraits,
trait::GeometryTraits,
geom,
) where {T<:IGeometry}
f = get(lookup_method, typeof(type), nothing)
isnothing(f) && error(
"Cannot convert an object of $(typeof(geom)) with the $(typeof(type)) trait (yet). Please report an issue.",
)
return f(GeoInterface.coordinates(geom))
f = get(lookup_method, typeof(trait), nothing)
isnothing(f) && _convert_error(geom, trait)
coords = GeoInterface.coordinates(geom)
@show typeof(coords)
return f(coords)
end

@noinline _convert_error(geom, trait) = error(
"Cannot convert an object of $(typeof(geom)) with the $(typeof(trait)) trait (yet). Please report an issue.",
)

function GeoInterface.geomtrait(
geom::Union{map(T -> AbstractGeometry{T}, pointtypes)...},
)
10 changes: 5 additions & 5 deletions src/ogr/geometry.jl
Original file line number Diff line number Diff line change
@@ -1,10 +1,11 @@
# Local convert method converts any GeoInterface geometry
# to the equivalent ArchGDAL geometry
# to the equivalent ArchGDAL geometry.
# `geom` can be points, so this is a hot path that needs to be fast.
function to_gdal(geom)
GeoInterface.isgeometry(geom) || _not_a_geom_error()
trait = GeoInterface.geomtrait(geom)
typ = geointerface_geomtype(trait)
GeoInterface.convert(type, trait, geom)
return GeoInterface.convert(typ, trait, geom)
end

@noinline _not_a_geom_error() =
@@ -942,11 +943,11 @@ multisurfaces (multipolygons).
"""
pointonsurface(geom::AbstractGeometry)::IGeometry =
IGeometry(GDAL.ogr_g_pointonsurface(geom))
pointonsurface(g1, g2) = pointonsurface(to_gdal(g1), to_gdal(g2))
pointonsurface(g) = pointonsurface(to_gdal(g))

unsafe_pointonsurface(geom::AbstractGeometry)::Geometry =
Geometry(GDAL.ogr_g_pointonsurface(geom))
unsafe_pointonsurface(g1, g2) = unsafe_pointonsurface(to_gdal(g1), to_gdal(g2))
unsafe_pointonsurface(g) = unsafe_pointonsurface(to_gdal(g1))

"""
difference(g1, g2)
@@ -1697,7 +1698,6 @@ MULTICURVE or MULTISURFACE in it, by approximating curve geometries.
* `options`: options as a null-terminated list of strings or NULL.
See OGRGeometryFactory::curveToLineString() for valid options.
"""

function lineargeom(
geom::AbstractGeometry,
stepsize::Real = 0;
4 changes: 3 additions & 1 deletion src/spatialref.jl
Original file line number Diff line number Diff line change
@@ -177,8 +177,10 @@ function reproject(
)
if GI.isgeometry(geom)
reproject(to_gdal(geom), sourcecrs, targetcrs; kwargs...)
elseif geom isa AbstractArray
elseif geom isa AbstractArray && (length(geom) > 0) && GI.isgeometry(first(geom))
reproject(to_gdal.(geom), Ref(sourcecrs), Ref(targetcrs); kwargs...)
else
throw(ArgumentError("geom is not a GeoInterface compatible geometry"))
end
end

10 changes: 9 additions & 1 deletion test/test_geometry.jl
Original file line number Diff line number Diff line change
@@ -1006,6 +1006,14 @@ import GeoFormatTypes as GFT

geom = MyLine()
ag_geom = GI.convert(AG.IGeometry, geom)
GI.coordinates(ag_geom) == [[1, 2], [1, 2]]
@test ag_geom isa AG.IGeometry{AG.wkbLineString}
@test GI.coordinates(ag_geom) ==
GI.coordinates(ag_geom) ==
[[1, 2], [1, 2]]
ag_geom1 = GI.convert(AG, geom)
@test ag_geom1 isa AG.IGeometry{AG.wkbLineString}
@test GI.coordinates(ag_geom1) ==
GI.coordinates(ag_geom) ==
[[1, 2], [1, 2]]
end
end
26 changes: 25 additions & 1 deletion test/test_geos_operations.jl
Original file line number Diff line number Diff line change
@@ -1,6 +1,9 @@
using Test
import ArchGDAL as AG
import GeoInterface
using Extents
const GI = GeoInterface


@testset "test_geos_operations.jl" begin
function equivalent_to_wkt(geom::AG.Geometry, wkt::String)
@@ -122,11 +125,15 @@ using Extents
@test AG.toWKT(result) == wkt2
end
@test AG.toWKT(f(geom)) == wkt2
if parentmodule(f) != GeoInterface && GI.ngeom(geom) != 0
wrapped_geom = GI.convert(GI, geom)
@test AG.toWKT(f(wrapped_geom)) == wkt2
end
end
end

@testset "Centroid" begin
test_method(AG.centroid, "POINT(10 0)", "POINT (10 0)")
test_method_gi(AG.centroid, "POINT(10 0)", "POINT (10 0)")
test_method(AG.centroid, "LINESTRING(0 0, 10 0)", "POINT (5 0)")
test_method(
AG.centroid,
@@ -164,6 +171,11 @@ using Extents
@test AG.toWKT(result) == wkt3
end
@test AG.toWKT(f(geom1, geom2)) == wkt3
if parentmodule(f) != GeoInterface && GI.ngeom(geom1) != 0 && GI.ngeom(geom2) != 0
wrapped_geom1 = GI.convert(GI, geom1)
wrapped_geom2 = GI.convert(GI, geom2)
@test AG.toWKT(f(wrapped_geom1, wrapped_geom2)) == wkt3
end
end
end
end
@@ -176,14 +188,26 @@ using Extents
AG.fromWKT(wkt1) do geom1
AG.fromWKT(wkt2) do geom2
@test AG.toWKT(f(geom1, geom2)) == wkt3
if parentmodule(f) != GeoInterface && GI.ngeom(geom1) != 0 && GI.ngeom(geom2) != 0
wrapped_geom1 = GI.convert(GI, geom1)
wrapped_geom2 = GI.convert(GI, geom2)
@test AG.toWKT(f(wrapped_geom1, wrapped_geom2)) == wkt3
end
end
end
end

function test_predicate(f::Function, wkt1, wkt2, result::Bool)
AG.fromWKT(wkt1) do geom1
AG.fromWKT(wkt2) do geom2
# Test GDAL geoms
@test f(geom1, geom2) == result
# Test GeoInterface geoms
if parentmodule(f) != GeoInterface && GI.ngeom(geom1) != 0 && GI.ngeom(geom2) != 0
wrapped_geom1 = GI.convert(GI, geom1)
wrapped_geom2 = GI.convert(GI, geom2)
@test f(geom1, geom2) == result
end
end
end
end