Polygonize a Raster

Specific Domains Geo
1

Hello there,
i know there is

ArchGDAL.polygonize|>methods
 [1] polygonize(f::Function, args...; kwargs...)
     @ 
 [2] polygonize(geom::ArchGDAL.AbstractGeometry)

and one could use PyCall or RCall to polygonize a raster.
Currently i use

gdal_polygonize_path = Conda.script_dir(Conda.ROOTENV)*"/gdal_polygonize.py"
cmd = pipeline(`python $(gdal_polygonize_path) $(input_raster_path) -f $(output_format) $(output_file_path)`)
run(cmd)

I was wondering if this works with pure Julia, preferably converting a Rasters.jl Raster to a VectorFormat…

Thank you

1 Like
2

My experience with gdalpoligonize (that you can acesa also from GMT.jl as “poliginize”) is that if not used with mask images (just 0s & 1s) is that it generates LOTS of polygons.

3

Hi @consumere

The new GeometryOps.jl has a polygonize method, but we haven’t added an extension there or to Rasters.jl to make it work easily on a Raster.

The method you want is:

polygonize(xs, ys, A::AbstractMatrix; minpoints=10)

Probably you can do:

using GeometryOps
polygonize(lookup(rast, X), lookup(rast, Y), rast)

And that will mostly work. But not its not properly handling Rasters.jl Intervals, e.g. with gdal files the lookup values are the start or the interval, not the middle.

We should add a GeometryOpsRastersExt extension to Rasters.jl to make this work correctly with just:

using GeometryOps
polygonize(rast)

I do work on both of these packages, but we haven’t gotten around to these integrations yet. Issue here:

4

Thank you @Raf for the quick response.

works on a AbstractMatrix but not on my data…
477Ă—285 Raster{Int32,2} Rasters.NoKW() extent: Extent(X = (525700.0, 621100.0), Y = (5.5461e6, 5.6031e6)) missingval: -9999 crs: EPSG:25832
i just get Any[] without an error

5

Hmm I’m not sure how that line can work on a non-DimensionalData.jl compatible AbstractMatrix…

It would be best to post a MWE. You need to include all coda and all data, generating rasters or downloading them in the script so we can copy, paste and it just runs without any other work.

Otherwise we can only guess at your problem.

6

You are right, i quoted the wrong line.
But polygonize seems to fail on negative cellvalues:

using Rasters
import GeometryOps
xlon = Rasters.X(525700.0:100:621100.0)
xlat = Rasters.Y(5.5461e6:100:5.6031e6)
rast = fill(-1, xlon, xlat)
pol = GeometryOps.polygonize(lookup(rast, X), lookup(rast, Y), rast)
#negative rastervalues return Any[]
rast = fill(10, xlon, xlat)
pol = GeometryOps.polygonize(lookup(rast, X), lookup(rast, Y), rast)
#gives 1-element Vector{GeometryBasics.  ...

how would i store the output to a GeoJson?
Thx

7

This works for writing. Unfortunately GeoJSON doesn’t handle a Vector of geometries, so wrap pol with a GeometryCollection

GeoJSON.write("test.json", GeometryOps.GeoInterface.GeometryCollection(pol))

We should with make GeoJSON handle a Vector, or make GeometryOps return a GeometryCollection or just MultiPolygon from polygonize. Its very early days for GeometryOps. It should be accurate and its insanely fast, but workflows like this are not polished yet.

Then read will get you a GeometryCollection back:

julia> GeoJSON.read("test.json")
GeometryCollection with 1 2D geometries

Issue here:

8

The negative cell value thing is not good:

polygonise is probably the least polished and tested method in GeometryOps, this is a good reminder to make it better. Ideally we would allow multiple algorithms too.

See what we can do, and thanks for all the feedback.

1 Like
9

Ok I just realised I have a branch from months ago that fixes half of this already! I will PR soon.

The output will be writeable to GeoJSON directly. Rasters.jl intergration will come after that.

10

That PR is merged and released - and we have tested that GeometryOps.polygonize handles Rasters correctly (as well as any array with custom axes).

11

I haven’t gotten around to adding the extension to Rasters to make it “just work” on a raster. But it mostly works if you get the lookups manually.

12

It seems to work here, but is there something missing there?