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?
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 Maybe it works for your use case…
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.
Isn’t this method for AbstractMultiPolygons?
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.
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.
Yeah, I saw this thread 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.
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!
Yeah, the WKT reading part was easy, it’s just converting into other formats that is not straightforward
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, 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:
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 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 is resolved.
2 Likes
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
# 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)
1 Like
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).
@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
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.
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. 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.
1 Like
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).
@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.
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).
That’s good to know. I’ll try to install again sometime.
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:
"""
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
1 Like