To clarify, Rasters.jl polygon ops are also pure julia via PolygonInbounds.jl and other internal code. It’s pretty fast for vectors of points, but not so much for single points. Its very fast for rasterize because the points of a raster are already sorted.
(Rasters.jl only uses GDAL for loading some file types, because GDAL has so many and is reliable, and for warping and reprojecting, because projection formats are hard)