# Rasterizing feature data

**URL:** <https://discourse.julialang.org/t/rasterizing-feature-data/59594>\
**Category:** Geo\
**Created:** [April 19, 2021, 1:53pm UTC](https://discourse.julialang.org/t/rasterizing-feature-data/59594 "2021-04-19T13:53:07Z")\
**Posts on this page:** 18\
**Page:** 1

<div class="post-metadata">

**Author:** ![mihalybaci](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mihalybaci/32/13528_2.png) [@mihalybaci](https://discourse.julialang.org/u/mihalybaci)\
**Post date:** [April 19, 2021, 1:53pm UTC](https://discourse.julialang.org/t/rasterizing-feature-data/59594/1 "2021-04-19T13:53:07Z")

</div>

I have a set of WKT strings that I would like to convert into 256x256 raster images (masks, really). I read them into geometries (specifically, multipolygons) just fine with `ArchGDAL.fomWKT`, but `ArchGDAL.unsafe_rasterize` (the only rasterize function I could find) fails because they are not `ArchGDAL.Dataset` types. Is there a way to convert the `IGeometry` type to `Dataset`? Or, is there a better way to rasterize feature data?

---

<div class="post-metadata">

**Author:** ![fabiangans](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/fabiangans/32/2624_2.png) [@fabiangans](https://discourse.julialang.org/u/fabiangans)\
**Post date:** [April 19, 2021, 3:51pm UTC](https://discourse.julialang.org/t/rasterizing-feature-data/59594/2 "2021-04-19T15:51:28Z")

</div>

Unfortunately I can not answer your ArchGDAL question, but I once coded a small `rasterize` function in Julia which operates on MultiPolygons, here is the link: [ESDL.jl/Shapes.jl at master · esa-esdl/ESDL.jl · GitHub](https://github.com/esa-esdl/ESDL.jl/blob/master/src/Shapes.jl#L104-L158) Maybe it works for your use case…

---

<div class="post-metadata">

**Author:** ![mihalybaci](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mihalybaci/32/13528_2.png) [@mihalybaci](https://discourse.julialang.org/u/mihalybaci)\
**Post date:** [April 19, 2021, 4:07pm UTC](https://discourse.julialang.org/t/rasterizing-feature-data/59594/3 "2021-04-19T16:07:23Z")

</div>

hmmm… it very well could work. The only problem is that it looks like I would need to convert the WKT strings to the `Shapefile.Polygon` type somehow.

---

<div class="post-metadata">

**Author:** ![fabiangans](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/fabiangans/32/2624_2.png) [@fabiangans](https://discourse.julialang.org/u/fabiangans)\
**Post date:** [April 19, 2021, 4:09pm UTC](https://discourse.julialang.org/t/rasterizing-feature-data/59594/4 "2021-04-19T16:09:18Z")

</div>

Isn’t [this](https://github.com/esa-esdl/ESDL.jl/blob/master/src/Shapes.jl#L104) method for AbstractMultiPolygons?

---

<div class="post-metadata">

**Author:** ![mihalybaci](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mihalybaci/32/13528_2.png) [@mihalybaci](https://discourse.julialang.org/u/mihalybaci)\
**Post date:** [April 19, 2021, 4:14pm UTC](https://discourse.julialang.org/t/rasterizing-feature-data/59594/5 "2021-04-19T16:14:46Z")

</div>

It is. I’m just not sure how to go from WKT strings (or `ArchGAL.IGeometry`) to the `Shapefile.AbstractMultiPolygon` type that the `rasterize!` method uses. Though, I might be missing something simple.

---

<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:** [April 19, 2021, 4:28pm UTC](https://discourse.julialang.org/t/rasterizing-feature-data/59594/6 "2021-04-19T16:28:55Z")

</div>

We really need to make this kind of thing easier… I’ve been putting off doing this in GeoData.jl forever.  
But it could be better if we make a simple shared package that does it well.

---

<div class="post-metadata">

**Author:** ![mihalybaci](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mihalybaci/32/13528_2.png) [@mihalybaci](https://discourse.julialang.org/u/mihalybaci)\
**Post date:** [April 19, 2021, 4:34pm UTC](https://discourse.julialang.org/t/rasterizing-feature-data/59594/7 "2021-04-19T16:34:21Z")

</div>

Yeah, I saw [this thread](https://discourse.julialang.org/t/how-can-i-rasterize-the-polygons-of-a-shapefile-in-julia/15089/17) from 2018 and was going to try it until I found the `AG.unsafe_rasterize` function. But now that I look at that old thread, I would still need to get my shapes into the `Dataset` type somehow.

It seems the natural place would be within `Shapefile.jl` since (unless I’m missing something) WKT is just a string version of a shapefile.

---

<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:** [April 19, 2021, 4:43pm UTC](https://discourse.julialang.org/t/rasterizing-feature-data/59594/8 "2021-04-19T16:43:16Z")

</div>

The issue there is reading WKT into e.g. `AbstractMultiPolygon`, which you would probably do with ArchGDAL unless someone writes a Julia WKT parser.

So it seems that conversion the `IGeometry` to polygon is the bottleneck. Probably @yeesian or @visr can direct you here. I would like to know too!

---

<div class="post-metadata">

**Author:** ![mihalybaci](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mihalybaci/32/13528_2.png) [@mihalybaci](https://discourse.julialang.org/u/mihalybaci)\
**Post date:** [April 19, 2021, 5:17pm UTC](https://discourse.julialang.org/t/rasterizing-feature-data/59594/9 "2021-04-19T17:17:36Z")

</div>

Yeah, the WKT reading part was easy, it’s just converting into other formats that is not straightforward

---

<div class="post-metadata">

**Author:** ![visr](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/visr/32/17204_2.png) [@visr](https://discourse.julialang.org/u/visr)\
**Post date:** [April 19, 2021, 6:59pm UTC](https://discourse.julialang.org/t/rasterizing-feature-data/59594/10 "2021-04-19T18:59:55Z")

</div>

I understand that you want to rasterize a WKT geometry, using `ArchGDAL.gdalrasterize`, but this expects a Dataset, whereas `ArchGDAL.fromWKT` returns a Geometry, right?

The GDAL docs explain the [Vector Data Model](https://gdal.org/user/vector_data_model.html), which is useful to understand the steps required. In short we need to create an empty Dataset, add a Layer to it, and for each Geometry put a Feature in this layer, and add the Geometry to the Feature.

For a single geometry this looks something like:

```julia
import ArchGDAL as AG

geom = AG.fromWKT("POINT (1 2)")
geomtype = AG.getgeomtype(geom)
dataset = AG.create(AG.getdriver("Memory"))
layer = AG.createlayer(;dataset, geom=geomtype, spatialref=AG.importEPSG(4326))

AG.createfeature(layer) do feature
    AG.setgeom!(feature, geom)
end

AG.gdalrasterize(AG.Dataset(dataset.ptr), ["-of","MEM","-tr","0.05","0.05"]) do ds_raster
    # e.g. save the raster
end

```

This isn’t quite as simple as it probably could be. I looked a bit into the [GeoDataFrames.jl IO code](https://github.com/evetion/GeoDataFrames.jl/blob/de2ffde3807e6fc6fb01b7843afe7f71f100a6f0/src/io.jl#L46-L68) for how to do this.

I tried to do it using the interactive IDataset, but had to work around them a bit in the last two calls, but that would probably be better if [RFC: Allow to use interactive Datasets in unsafe\_gdal functions in utilities.jl by felixcremer · Pull Request #167 · yeesian/ArchGDAL.jl · GitHub](https://github.com/yeesian/ArchGDAL.jl/pull/167) is resolved.

---

<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:** [April 19, 2021, 7:05pm UTC](https://discourse.julialang.org/t/rasterizing-feature-data/59594/11 "2021-04-19T19:05:37Z")

</div>

Hm, I feel I’ll get into troubles because this is WIP work and many bugs floating around (making this example already revealed a couple) but can’t resist, so here it goes. With current GMT.jl you can do

```julia
# For sake of example, create a GMT dataset with a square line
D = mat2ds([1. 1; 1 2; 2 2; 2 1; 1 1]);

# Save it as a shapefile
ogr2ogr(D, save="square.shp");

# Read it back as a GDAL dataset
ds = readgd("square.shp");

# Rasterize it (just picked something similar to gdal_rasterize man page)
ds_ras = gdalrasterize(ds, ["-burn", "255", "-burn", "0", "-burn", "0", "-ts", "512", "512", "-te", "0", "0", "3", "3", "-ot", "Byte"]);

# Save it as a PNG file
gdaltranslate(ds_ras, save="square.png")

# See it (one GMT bug here. One should be able to do directly imshow(ds_ras))
imshow("square.png", fmt=:png)

```

 ![GMTjl_tmp](https://global.discourse-cdn.com/julialang/original/3X/a/a/aad968602de2353d95abe88ffe30ae73de6fa0a5.png)

---

<div class="post-metadata">

**Author:** ![mihalybaci](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mihalybaci/32/13528_2.png) [@mihalybaci](https://discourse.julialang.org/u/mihalybaci)\
**Post date:** [April 20, 2021, 4:32pm UTC](https://discourse.julialang.org/t/rasterizing-feature-data/59594/12 "2021-04-20T16:32:48Z")

</div>

Thanks for the replies!

@joa-quim - This seems pretty easy, though I don’t have GMT installed because I couldn’t get past versioning issues with GMT and GDAL (mentioned in [this thread](https://discourse.julialang.org/t/gmt-jl-errors-when-trying-to-plot-problem-with-png/47860)).

@visr - This method seems easy enough as well, I just don’t have enough experience with GDAL (thought I could fumble around a bunch). I tried to save out the resulting raster with

```julia
AG.gdalrasterize(AG.Dataset(dataset.ptr), ["-of","MEM","-tr","0.05","0.05"]) do ds_raster
    # e.g. save the raster
    println(ds_raster)
    AG.write(ds_raster, "test.tif")
end

```

but no file gets saved. Is this the right way to go about it? The result of `println` tells me the resulting file is a 1x1 dataset, even when I use my real WKT, which seems odd.

---

<div class="post-metadata">

**Author:** ![visr](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/visr/32/17204_2.png) [@visr](https://discourse.julialang.org/u/visr)\
**Post date:** [April 20, 2021, 4:49pm UTC](https://discourse.julialang.org/t/rasterizing-feature-data/59594/13 "2021-04-20T16:49:34Z")

</div>

Ah yes, by passing the arguments `"-of","MEM"` to `gdalrasterize`, we are telling it to create an in memory raster, that cannot be written to disk. The `println(ds_raster)` also confirms this (`GDAL Dataset (Driver: MEM/In Memory Raster)`. To get a GeoTIFF, use `"-of","GTiff"`. Then your `AG.write` call will do something.

Rasterizing a single point will result in a 1x1 raster, unless you give both explicit resolution and extent with the `-tr` and `-te` flags described in the [gdal\_rasterize docs](https://gdal.org/programs/gdal_rasterize.html). I don’t know your real geometry, but if the extent is below the 0.05 resolution you specify, that would be expected as well.

---

<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:** [April 20, 2021, 4:56pm UTC](https://discourse.julialang.org/t/rasterizing-feature-data/59594/14 "2021-04-20T16:56:06Z")

</div>

> [@mihalybaci](#):
>
> @joa-quim - This seems pretty easy, though I don’t have GMT installed because I couldn’t get past versioning issues with GMT and GDAL (mentioned in [this thread](https://discourse.julialang.org/t/gmt-jl-errors-when-trying-to-plot-problem-with-png/47860)).

Hmm, there was nothing GDAL related in that thread and the issue, related to a bad PS closing for older GMT versions, was fixed months ago (you even said that all tests run fine).

---

<div class="post-metadata">

**Author:** ![mihalybaci](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mihalybaci/32/13528_2.png) [@mihalybaci](https://discourse.julialang.org/u/mihalybaci)\
**Post date:** [April 20, 2021, 5:13pm UTC](https://discourse.julialang.org/t/rasterizing-feature-data/59594/15 "2021-04-20T17:13:23Z")

</div>

@visr - Ah, great, that works now. After looking at the GDAL docs, I tried putting the filename in place of “MEM” thinking it would recognize it. I am trying to make image masks, so I’ll need to add the size flags.

@joa-quim - Not to get too far off topic, but I think that’s because I had resolved the GDAL before finding that other issue. If I fully remember the original problem, the system version of GDAL (currently v2.4) for my OS (MX Linux) was too old to work with GMT. I tried installing the latest GDAL from the website, but then I didn’t have compatible system libraries so it still didn’t work. In the end I had to completely change OS’s to Ubuntu (which I had wanted to do anyway), and that was the only way I got both GDAL and GMT to install properly. Now, after not-so-great Ubuntu experience, I am back on MX Linux, so I am hesitant to manually install GDAL again.

---

<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:** [April 20, 2021, 5:43pm UTC](https://discourse.julialang.org/t/rasterizing-feature-data/59594/16 "2021-04-20T17:43:03Z")

</div>

OK, but just a note. GMT works with any GDAL you link it with. Only that many features won’t be available when GDAL version is old (and may even error if care was not taken to shelter those cases).

---

<div class="post-metadata">

**Author:** ![mihalybaci](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mihalybaci/32/13528_2.png) [@mihalybaci](https://discourse.julialang.org/u/mihalybaci)\
**Post date:** [April 20, 2021, 5:45pm UTC](https://discourse.julialang.org/t/rasterizing-feature-data/59594/17 "2021-04-20T17:45:18Z")

</div>

That’s good to know. I’ll try to install again sometime.

---

<div class="post-metadata">

**Author:** ![mihalybaci](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mihalybaci/32/13528_2.png) [@mihalybaci](https://discourse.julialang.org/u/mihalybaci)\
**Post date:** [April 22, 2021, 12:27pm UTC](https://discourse.julialang.org/t/rasterizing-feature-data/59594/18 "2021-04-22T12:27:18Z")

</div>

I think I have working version now. My output .tif matches WKT representation in QGIS after switching to GTiff output. I did also have to add the “burn” keyword like @joa-quim suggested, otherwise the files were just all zeros. Thanks for your help everyone!

Edit: Here is the final result:

```julia
"""
Arguments for the rasterizeWKT function
"""
Base.@kwdef struct RasterArgs
    saveas = "geom.tif" # Output file(path)
    burn = "1" # Feature fill value
    ext = ["0", "0", "1", "1"] # Image extent
    size = ["256", "256"]
end

"""
Rasterize a WKT string
"""
function rasterizeWKT(wkt::String; kwargs...)
    args = RasterArgs(; kwargs...)
    geom = AG.fromWKT(wkt)
    
    geomtype = AG.getgeomtype(geom)

    dataset = AG.create(AG.getdriver("Memory"))
    layer = AG.createlayer(;dataset, geom=geomtype, spatialref=AG.importEPSG(4326))
    
    AG.createfeature(layer) do feature
        AG.setgeom!(feature, geom)
    end
    
    gdal_kws = ["-at", # Set "all touched" to true
                "-burn", args.burn, # Set burn value 
                "-of", "GTiff", # Save as geotiff 
                "-te", args.ext..., # Define image extent 
                "-ts", args.size...] # Define image size (pixels)
    AG.gdalrasterize(AG.Dataset(dataset.ptr), gdal_kws) do ds_raster
        AG.write(ds_raster, args.saveas)
    end

    return nothing
end

```