You may be able to use Rasters.coverage for each land use separately and then combine them with maximum?
I guess we could build that into rasterize…
Is there another GIS function somewhere that does this we could copy?
Yes, geometries are not squares, so I may have a destination raster pixel which is “covered” by 40% forest (id=1), 30% agriculture (id=2) and 30% artificial area (id=3), and in that case the pixel should have the value 1.
I’ll try the Rasters.coverage…
Thank you, this seems to achieve what I was looking for:
clc2018df = GeoDataFrames.read(clc2018_path)
DataFrames.rename!(clc2018df,"Shape" => "geometry")
clc_311 = clc2018df[clc2018df.Code_18 .== "311",:]
clc_312 = clc2018df[clc2018df.Code_18 .== "312",:]
clc_313 = clc2018df[clc2018df.Code_18 .== "313",:]
clc_311_share = Rasters.coverage(clc_311; to=mask) # `mask` is my raster grid mask
clc_312_share = Rasters.coverage(clc_312; to=mask)
clc_313_share = Rasters.coverage(clc_313; to=mask)
Just a nuisance, the coverage command took forever because lots of info were printed :
[ Info: A geometry was ignored at Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}((X = (3.4718729800000004e6, 3.47309808e6), Y = (1.991912e6, 1.9929848200000003e6))) as it was outside of the supplied extent Extents.Extent{(:X, :Y), Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}((X = (2.984e6, 4.322e6), Y = (1.993e6, 3.211e6)))
In my case this is ok, as the vector is way larger than the area under study. I think the info message should be printed only if verbose is set to true.
1 Like
You can use Oh you mean that still doesn’t stop it. Yes we should fix that: Dont print coverage info with `verbose=true` · Issue #680 · rafaqz/Rasters.jl · GitHubverbose=false to turn off the info printing (check the docs for all these keywords)
coverage will always be kinda slow compared to rasterize because it actually has to rasterize 100 points for each pixel. But still its very likely that Rasters.coverage is the fastest coverage function available in any language.
It would be good to somehow put your code into a single function in Rasters.jl, like rasterize(last, geoms; by_coverage=maximum, kw...) could do it?
(I think we need both last and maximum in case there are multiple with the same coverage)