Hello, being only marginally on the domain, I am a bit lost with the various GIS packages…
I need to rasterize a geopackage (gpkg) vector file (Corine Land Cover).
Visualizing the gpkg file in QGis I can see there are various layers in a single file (in this case the main one is the European land use cover, the one I am interested in, and the others refer to overseas French territories).
My final objective is to take the rasterization at high resolution and count the various pixels of a given class within a low resolution pixel, so I have as final product a set of rasters, each one indicating the percentual of land use for a specific class within each pixel.
I think I can do the second step, but I need first to do the first rasterization. I have tried in qgis to trsnform the geopackage layer as Shapefile, for which there is an example of Rasterization in Rasters.jl, but QGis complain that the file is too big for the Shapefile format.
You should be able to open the geopackage using ArchGDAL and then it will be the same as the shapefile example for rasterize. There is no need to convert the data beforehand, because both packages adhere to the GeoInterface and therefore these geo types are interchangeable.
Yes ArchGDAL or GeoDataframes.jl will open that fine, or also GeoIO.jl.
Then Rasters.jl should just rasterize any of those objects directly in rasterize without modification, because GeoInterface.jl handles the conversions for you.
There are no comprehensive benchmarks at this stage but it looks like Rasters.rasterize its one of the fastest implementations that exists. If you find anything else as fast, make an issue
yes, it is a multi-vector layer.
The read call of ArchGDAL just get the whole package, but then one has to use, as I discovered, use ArchGDAL.getlayer() to get the specific layer.
I will try 2 more attempts now:
use GeoDataframes.jl to load the layer
trying to use the Rasterize transform from GeoStats.jl on the ArchGDAL.getlayer output…
Edit:
First attempt succeeded in creating a df (got a warning " Warning: This file has multiple layers, you only get the first layer by default now." but that’s fine for me).
I obtained a df that looks like:
However passing this df to Rasters.rasterize I got the same error as the gdal object (ERROR: ArgumentError: Object is not a GeoInterface.jl compatible geometry: ())
Ok, prepared a hotfix to handle this unexpected geometry column name:
# helper function to find the
# geometry column of a table
function geomcolumn(names)
snames = string.(names)
gnames = ["geometry","geom","shape"]
gnames = [gnames; uppercasefirst.(gnames)]
gnames = [gnames; uppercase.(gnames)]
gnames = [gnames; [""]]
select = findfirst(∈(snames), gnames)
if isnothing(select)
throw(ErrorException("geometry column not found"))
else
Symbol(gnames[select])
end
end
Managed to load the file without problems with the master branch of GeoIO.jl:
@sylvaticus if you decide to try the Rasterize transform on the loaded geotable with GeoIO.jl, please feel free to report further issues. These reports really help improve robustness of the packages.
If people want to use another package that’s fine, but honestly I find these threads with two competing methods from different packages are pretty distracting.
The question was for Rasters.jl and the solution is very simple. Let’s not do this all the time. We can wait to see if there is no easy solution to the direct question.
In this case there is a keyword for exactly this problem.
I wouldn’t join the thread with an alternative solution if the title was requesting Rasters.jl specifically. I thought that @sylvaticus just wanted to accomplish a generic task as a first-time user of GIS software.
If this is a Rasters.jl question, sorry for the confusion.
I am very sorry, I didn’t want to spark an issue between GIS packages.
Just GIS is not my main domain, but I needed to do “some GIS” within a larger model.
I have manually renamed the df outputed by GeoDataFrames and this then indeed worked with Rasters.Rasterize(), I will also try the GeoIO → GeoStats.Rasterize method and report if it works…
I have a related issue but I am now afraid to ask… :-/ How do I rasterize by setting the pixel value the one with the majority land use in the pixel (land use is in the vector format in the Code_18 column).
