speed of GeoInterface.extent could be improved for large collections

[alex-s-gardner]

Load packages

using GeoInterface
using GeometryOps

Download file

# for convenience I've temporarily put the files here
fn = "02_rgi60_WesternCanadaUS"
url2shp =  "https://its-live-data.s3.amazonaws.com/test/$fn.zip"

Downloads.download(url2shp, "$fn.zip")
run(`unzip $fn`)

Load polygons

glacier = Shapefile.Handle("$fn/$fn.shp")

Find extents

here's where it gets interesting

calculate extents on raw file

@time GeoInterface.extent(glacier)
0.000008 seconds (1 allocation: 64 bytes)
Extent(X = (-133.73235307, -105.60821024699995), Y = (36.38625038500004, 65.15664046000005), Z = (0.0, 0.0))

I guess this is just grabbing the extents from the header... that's the only place I can see extent info in glacier .

now retroject and get extents [requires https://github.com/JuliaGeo/Proj.jl/pull/97]

epsg = GeoFormatTypes.EPSG("EPSG:32609")
glacier_rep = GeometryOps.reproject(glacier; target_crs=epsg);

@time GeoInterface.extent(glacier_rep)
  3.575146 seconds (14.00 M allocations: 714.567 MiB, 11.20% gc time)
Extent(X = (273581.3643637521, 2.5071260739734517e6), Y = (4.0779523613324095e6, 7.23421013737856e6))

3.5s seems like a very long time for this operation. I suspect that the Extents if being calculated on each feature then the union of extents is being taken.. if this is the case it would probably be better to load all feature nodes in as points then take the extrema once.

[rafaqz]

Yeah, theres something wrong with that time.

Try ProfileView on it and see its all type instability, and what function its in.

As an aside, reproject should probably update the extent as it runs.

[rafaqz]

Turns out this is already fixed on master we just forgot to register it.

Should be 50x faster but it looks like it still isn't completely type stable.

[rafaqz]

Fixed in the latest release.