# Polygonize a Raster

**URL:** <https://discourse.julialang.org/t/polygonize-a-raster/114189>\
**Category:** Geo\
**Created:** [May 13, 2024, 7:42am UTC](https://discourse.julialang.org/t/polygonize-a-raster/114189 "2024-05-13T07:42:09Z")\
**Posts on this page:** 12\
**Page:** 1

<div class="post-metadata">

**Author:** ![consumere](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/consumere/32/209031_2.png) [@consumere](https://discourse.julialang.org/u/consumere)\
**Post date:** [May 13, 2024, 7:42am UTC](https://discourse.julialang.org/t/polygonize-a-raster/114189/1 "2024-05-13T07:42:09Z")

</div>

Hello there,  
i know there is

```julia
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

```julia
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

---

<div class="post-metadata">

**Author:** ![joa-quim](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/joa-quim/32/227_2.png) [@joa-quim](https://discourse.julialang.org/u/joa-quim)\
**Post date:** [May 13, 2024, 8:26am UTC](https://discourse.julialang.org/t/polygonize-a-raster/114189/2 "2024-05-13T08:26:59Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![Raf](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/raf/32/3383_2.png) [@Raf](https://discourse.julialang.org/u/Raf)\
**Post date:** [May 13, 2024, 9:00am UTC](https://discourse.julialang.org/t/polygonize-a-raster/114189/3 "2024-05-13T09:00:40Z")

</div>

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:

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

```

Probably you can do:

```julia
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:

```julia
using GeometryOps
polygonize(rast)

```

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

> <https://github.com/rafaqz/Rasters.jl/issues/668>
>
> So that these methods know about the dimensions and crs of the \`Raster\`.
> 
> We c…ould also just add a hard dependency and start using \`apply\` and other useful things internally here.
> 
> I think its best that Rasters dependends on/extends GeometryOps, and not the other way around? @asinghvi17 @skygering ?

---

<div class="post-metadata">

**Author:** ![consumere](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/consumere/32/209031_2.png) [@consumere](https://discourse.julialang.org/u/consumere)\
**Post date:** [May 13, 2024, 10:02am UTC](https://discourse.julialang.org/t/polygonize-a-raster/114189/4 "2024-05-13T10:02:56Z")

</div>

Thank you @Raf for the quick response.

> [@Raf](#):
>
> `polygonize(lookup(rast, X), lookup(rast, Y), rast)`

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

---

<div class="post-metadata">

**Author:** ![Raf](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/raf/32/3383_2.png) [@Raf](https://discourse.julialang.org/u/Raf)\
**Post date:** [May 13, 2024, 10:58am UTC](https://discourse.julialang.org/t/polygonize-a-raster/114189/5 "2024-05-13T10:58:55Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![consumere](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/consumere/32/209031_2.png) [@consumere](https://discourse.julialang.org/u/consumere)\
**Post date:** [May 13, 2024, 12:06pm UTC](https://discourse.julialang.org/t/polygonize-a-raster/114189/6 "2024-05-13T12:06:15Z")

</div>

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

```julia
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

---

<div class="post-metadata">

**Author:** ![Raf](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/raf/32/3383_2.png) [@Raf](https://discourse.julialang.org/u/Raf)\
**Post date:** [May 13, 2024, 12:23pm UTC](https://discourse.julialang.org/t/polygonize-a-raster/114189/7 "2024-05-13T12:23:33Z")

</div>

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

```julia
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](https://github.com/JuliaGeo/GeometryOps.jl?tab=readme-ov-file#performance-comparison-to-other-packages), but workflows like this are not polished yet.

Then read will get you a GeometryCollection back:

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

```

Issue here:

> <https://github.com/JuliaGeo/GeometryOps.jl/issues/139>
>
> See: https://discourse.julialang.org/t/polygonize-a-raster/114189/7

> <https://github.com/JuliaGeo/GeoJSON.jl/issues/92>
>
> It would be nice if arrays/iterables of geometries would just write as a \`Geomet…ryCollection\` without having to explicitly wrap them as that, as its a common way to hold and manipulate geometries you are working on.

---

<div class="post-metadata">

**Author:** ![Raf](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/raf/32/3383_2.png) [@Raf](https://discourse.julialang.org/u/Raf)\
**Post date:** [May 13, 2024, 12:30pm UTC](https://discourse.julialang.org/t/polygonize-a-raster/114189/8 "2024-05-13T12:30:13Z")

</div>

The negative cell value thing is not good:

> <https://github.com/JuliaGeo/GeometryOps.jl/issues/140>
>
> MWE:
> 
> \`\`\`julia
> using Rasters
> import GeometryOps
> xlon = X(525700.0:100:62110…0.0)
> xlat = 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. ...
> \`\`\`

`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.

---

<div class="post-metadata">

**Author:** ![Raf](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/raf/32/3383_2.png) [@Raf](https://discourse.julialang.org/u/Raf)\
**Post date:** [May 13, 2024, 12:58pm UTC](https://discourse.julialang.org/t/polygonize-a-raster/114189/9 "2024-05-13T12:58:50Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![asinghvi17](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/asinghvi17/32/8272_2.png) [@asinghvi17](https://discourse.julialang.org/u/asinghvi17)\
**Post date:** [June 16, 2024, 1:51pm UTC](https://discourse.julialang.org/t/polygonize-a-raster/114189/10 "2024-06-16T13:51:54Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![Raf](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/raf/32/3383_2.png) [@Raf](https://discourse.julialang.org/u/Raf)\
**Post date:** [June 16, 2024, 2:03pm UTC](https://discourse.julialang.org/t/polygonize-a-raster/114189/11 "2024-06-16T14:03:23Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![asinghvi17](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/asinghvi17/32/8272_2.png) [@asinghvi17](https://discourse.julialang.org/u/asinghvi17)\
**Post date:** [June 18, 2024, 6:42pm UTC](https://discourse.julialang.org/t/polygonize-a-raster/114189/12 "2024-06-18T18:42:53Z")

</div>

It seems to work [here](https://github.com/JuliaGeo/GeometryOps.jl/blob/0d15e9c3f57e6180c6af86ea5343721505f83766/test/methods/polygonize.jl#L61-L68), but is there something missing there?