Add geofileops to speed test?

[theroggy]

Hello!

I encountered this repo by coincidence: really interesting excercise.

You might also be interested in having a look at geofileops. At the bottom of the readme of the repo there are also some benchmarks. As you can see, intersection isn't the operation where the most spectacular speed improvements are to be expected because it isn't the "heaviest" spatial operation but...

It works different compared to the tools you already tested because it works file-based rather than in-memory... GPKG has some specific properties that make this possible, so for spatial operations the files should be in GPKG. I've never really looked optimized to get memory usage explicitly low, but because files are not taken in memory even really large files (larger than memory available) should work fine.

For illustration, I wrote a rough "geofileops" version of your test and ran it on my ~desktop: a windows virtual machine with 8 virtual CPU's and 32 GB RAM:

• Because your inputs are in many files seperate files, they are first merged in a single gpkg file: this took 190 seconds.
• Then to actually calculate the intersection the load/intersection/save are in one go, and this took 105 seconds.

"""Test the speed of intersection with geofileops."""

import logging
import time
import warnings
from pathlib import Path

import geofileops as gfo

if __name__ == "__main__":
    logging.basicConfig(level=logging.INFO)
    warnings.filterwarnings("ignore")

    start = time.time()

    building_cols = [
        "oid",
        "aktualit",
        "gebnutzbez",
        "funktion",
        "anzahlgs",
        "gmdschl",
        "lagebeztxt",
        # "geometry",
    ]
    parcels_cols = [
        "oid", "aktualit", "nutzart", "bez", "flstkennz", # "geometry"
    ]

    alkis_dir = Path("H:/temp/ALKIS")
    buildings_paths = list(alkis_dir.glob("./*/GebauedeBauwerk.shp"))
    parcels_paths = list(alkis_dir.glob("./*/NutzungFlurstueck.shp"))

    buildings_path = alkis_dir / "GebauedeBauwerk.gpkg"
    if not buildings_path.exists():
        for path in buildings_paths:
            gfo.copy_layer(src=path, dst=buildings_path, dst_layer=buildings_path.stem, append=True, create_spatial_index=False)
        gfo.create_spatial_index(buildings_path)

    parcels_path = alkis_dir / "NutzungFlurstueck.gpkg"
    if not parcels_path.exists():
        for path in parcels_paths:
            gfo.copy_layer(src=path, dst=parcels_path, dst_layer=parcels_path.stem, append=True, create_spatial_index=False)
        gfo.create_spatial_index(parcels_path)

    # buildings_gdf = buildings_gdf.drop_duplicates(subset="oid", keep="first")
    # parcels_gdf = parcels_gdf.drop_duplicates(subset="oid", keep="first")
    print(f"geofileops: Prepare data duration: {(time.time() - start):.0f} s.")

    start_intersection = time.time()
    buildings_with_parcels_path = alkis_dir / "buildings_with_parcels.gpkg"
    gfo.intersection(buildings_path, parcels_path, buildings_with_parcels_path, input1_columns=building_cols, input2_columns=parcels_cols)
    print(f"geofileops: Load, intersection, save takes: {(time.time() - start_intersection):.0f} s.")

    print(f"geofileops: Total duration: {(time.time() - start):.0f} s.")

[sehHeiden]

Thanks. I already Rad about geofileops. When I remember correctly it uses multiprocessing on sliced geodata frames?

I will leave it open as an reminder, to try it out when I find the time to execute on both machines and append it as chapter to the corresponding Blog Post.

[theroggy]

Thanks. I already Rad about geofileops. When I remember correctly it uses multiprocessing on sliced geodata frames?

No, for most spatial operations it doesn't use dataframes. Typically it uses SQL statements to filter the data directly using the spatial index available in GPKG files to avoid having to read/keep all data in memory.

[sehHeiden]

@theroggy I added a ci that execute the benchmarks. Currently, it's still a part of the dataset on Ubuntu latest 2 core. But I can't get this: https://github.com/sehHeiden/geospeed/blob/master/geospeed/geofileops.py

running. From my current understanding, it is a problem with the gdal version vs geofileops version.
See the install process here:
https://github.com/sehHeiden/geospeed/blob/master/.github/workflows/benchmark.yml
Do have any tips?

Btw.: If it is running, I would like to use a larger machine.

[theroggy]

Not sure, I actually always use conda, but it looks to me that the wrong version of the GDAL python bindings are being installed:

This line seems to resolve to GDAL 3.8.5, while it should be 3.8.4 as that's the version of the binaries that are installed.

uv add 'geofileops>=0.8.0' "GDAL>=$GDAL_MAJOR_MINOR.0

So I guess changing the GDAL version like this might solve it:

uv add 'geofileops>=0.8.0' "GDAL==$GDAL_MAJOR_MINOR.4

[sehHeiden]

Hi, I just merged the request.

Not so familiar with the code base anymore, but looks like geofileops runs successfully.
Just the speedup is not so high, as other tools. Also the Ram usage is higher.
Perhaps these are factors that help you with improving geofileops.

If you don't have any other comments, I gonna close this issue after New Year. To give you some time to answer.

[sehHeiden]

@theroggy just want to make sure you get a message and will read my last message.

[theroggy]

Hi, I just merged the request.

Great!

Not so familiar with the code base anymore, but looks like geofileops runs successfully. Just the speedup is not so high, as other tools.

I updated to the newest version of geofileops in #19

Now the test is faster, but I think that the small test file in the CI run is difficult to draw conclusions on in general. I see relatively large differences between runs in all benchmarks where the differences between runs are sometimes larger than the differences between the different libraries... So, the CI test is great to see if everything works fine, but drawing conclusions about performance seems difficult.

Also the Ram usage is higher. Perhaps these are factors that help you with improving geofileops.

Comparing RAM usage isn't super easy. I added an additional memory counter, and tried to explain the difference in this PR: #21

To make it easier to compare results in a PR, I added a log of the result in the "benchmark" action in this PR: #20

[sehHeiden]

@theroggy yes using peak_memory_less_available_mb variable sounds sound...
Perhaps that should be done in general and also used for the final table.
yes counting Ram usage is tricky, as you can easily miss a peak. Otherwise don't want to do so often AND on these tiny 2 core runners this easily creates new problems when tasks or processes do run in parallel.
Also as you said the timing can not be stable on a tiny virtual machine. You never know how much actual load there is due to other users and cores are clocking down etc.
I started it as a test. But in best cases it could show effects of updated libraries. But for that you would need repetition of the tasks and probably a dedicated server for that time of the run.

As the geometries are really simple, this also a best case where tree indixes can shine.
With long or complex geometries the use of smaller files should even be more beneficial.

[sehHeiden]

@theroggy what was worrying about the Ram use is that it can be a hint of not optimized operations. Loading and writing from ram is so much slower than pure compute.
The other fear I did have, that time is spend with writing files and therefore the code is slower.

[theroggy]

As the geometries are really simple, this also a best case where tree indixes can shine. With long or complex geometries the use of smaller files should even be more beneficial.

Yes, an intersection between "simple" geometries is indeed not a super demanding/complex use case. In the performance benchmarking I do myself, between geopandas and geofileops, this is indeed the use case with the least gain:
https://github.com/geofileops/geobenchmark/blob/main/results_vector_ops/GeoBenchmark.png

what was worrying about the Ram use is that it can be a hint of not optimized operations.

OK. As shown with the difference in free memory it is not the case that geofileops uses that much memory, it's mainly double-counting. When you run on large files you will even see that geofileops most likely will use only a fraction of the memory of most other options.

Loading and writing from ram is so much slower than pure compute.

That depends. In GIS there are a lot of CPU-intensive operations where the processing takes a lot longer than reading RAM or even reading from SSD/Disk. See the benchmarks I linked to above, e.g. for dissolve or overlays with complex geometries even the I/O to disk is negligable compared to CPU time...

The other fear I did have, that time is spend with writing files and therefore the code is slower.

There are always tradeoffs to be made... but the big advantage of working file-based is that memory usage for large files is a fraction of solutions that pull everything in memory... and parallellization reduces the impact of this significantly. But, this should become clear with a test on larger files...