How to properly benchmark/profile to disambiguate IO/function call overhead

[AndrewAnnex]

Intro

I am producing web-optimized COGs with rio-cogeo and requesting specific tiles using the same TMS used to create the cogs with rio-tiler but I am getting pretty poor performance when using the standard rasterio reader methods in rio-tiler of about half a second per tile. This is with using a single cog, the Reader class using the tms, and the .tile method to look up a single tile at the native zoom level (highest zoom level defined in the COG), so these tiles should be essentially the same as the blocks within the COG. However I ran this test with a cog on a local NVME disk so I don't think that IO is the issue.

Main Question:

As I am still trying to understand this, my main question is: What are the best practices to profile/benchmark rio-tiler? Any tools/techniques to share for how to differentiate between IO/compute issues would be great.

Now, it could be that the performance of rio-tiler is already good enough for what it aims to do, and that my specific use case is so specific that rio-tiler isn't the right library to use. This is why I haven't posted this as an issue, but I hope this question could lead to some documentation additions/best practices that others could benefit from.

What I've looked into so far:

I found some prior issues/discussions that mentioned activating various GDAL and Proj debug environment variables CPL_DEBUG=ON, PROJ_DEBUG=3, (I even used CPL_LOG to capture messages) externally before starting the ipython repl or running any scripts I made, but I didn't see anything notable from them (for a single tile read, just two quick calls to pj_open_lib and then gdalClose). Disabling GDAL_DISABLE_READDIR_ON_OPEN didn't matter much/at all.

As I initially thought this was an IO issue, but I get the same performance when I moved the file over to a NVME disk. Doing some quick checks, I attempted to use rasterio directly, converting the TMS tiles to the correct Window object and with boundless=False enabled I got to about 10000 iterations per second. After enabling boundless request this drops by about an order of magnitude to 1000 iterations per second. This compared to the initial ~2 iterations per second tells me something performance bound is more likely the cause.

I then wrote a little side script and used cProfile to profile a single tile lookup and attached the log below. The top of the file with calls that actually gave perceptible amounts of time is:

         228 function calls (223 primitive calls) in 0.440 seconds

   Ordered by: cumulative time

   ncalls  tottime  percall  cumtime  percall filename:lineno(function)
        1    0.000    0.000    0.440    0.440 /home/aannex/micromamba/envs/marspre/lib/python3.12/site-packages/rio_tiler/io/rasterio.py:317(tile)
        1    0.000    0.000    0.326    0.326 /home/aannex/micromamba/envs/marspre/lib/python3.12/site-packages/rio_tiler/io/rasterio.py:365(part)
        1    0.217    0.217    0.326    0.326 /home/aannex/micromamba/envs/marspre/lib/python3.12/site-packages/rio_tiler/reader.py:296(part)
        1    0.113    0.113    0.114    0.114 /home/aannex/micromamba/envs/marspre/lib/python3.12/site-packages/rio_tiler/io/base.py:88(tile_exists)
        1    0.108    0.108    0.109    0.109 /home/aannex/micromamba/envs/marspre/lib/python3.12/site-packages/rio_tiler/reader.py:85(read)
        1    0.000    0.000    0.001    0.001 /home/aannex/micromamba/envs/marspre/lib/python3.12/functools.py:977(__get__)
        1    0.000    0.000    0.001    0.001 /home/aannex/micromamba/envs/marspre/lib/python3.12/site-packages/morecantile/models.py:560(rasterio_crs)
        1    0.001    0.001    0.001    0.001 

Looking at this it looks like the half second time is accurate, and that a big chunk of the delay could be in the tile_exists function. The rest seems to be in the part/read calls in reader.py. But I'm not sure yet how to interpret past that point.

Using snakeviz I get a similar view, but I am not sure why I don't see more depth of function calls here.

At this point I am not sure what more to look at, but I wonder if there isn't some opportunity to inspect the metadata for a given cog, see if it is web optimized, and go down some other code path to more optimally query it before performing whatever post processing is needed or requested at the higher level calls. But maybe this is already happening?

rio_tiler_log.txt

[vincentsarago]

@AndrewAnnex could you share the files you're using?

I usually use https://github.com/developmentseed/tilebench?tab=readme-ov-file#api to benchmark

the timing of tile_exist is pretty weird, if the CRS of the COG and the TMS match the tile_exists method should just be in

rio-tiler/rio_tiler/io/base.py

Lines 199 to 204 in cd039d3

	return (
	(tile_bounds[0] < dst_bounds[2])
	and (tile_bounds[2] > dst_bounds[0])
	and (tile_bounds[3] > dst_bounds[1])
	and (tile_bounds[1] < dst_bounds[3])
	)

, I can't think about any reason why this would take 100ms

I attempted to use rasterio directly, converting the TMS tiles to the correct Window object and with boundless=False enabled I got to about 10000 iterations per second. After enabling boundless request this drops by about an order of magnitude to 1000 iterations per second.

Interesting, again is the COG and the TMS have the same CRS, we won't use WarpedVRT, and boundless would depends if the windows is crossing the dataset bounds

rio-tiler/rio_tiler/reader.py

Lines 176 to 189 in 09ed5c3

	boundless = False
	if window:
	if isinstance(window, tuple):
	window = windows.Window.from_slices(
	*window, height=dataset.height, width=dataset.width, boundless=True
	)
	(row_start, row_stop), (col_start, col_stop) = window.toranges()
	if (
	min(col_start, row_start) < 0
	or col_stop >= dataset.width
	or row_stop >= dataset.height
	):
	boundless = True

Again, if you can share both TMS and COG and also the code you're using it would be useful

[AndrewAnnex]

okay here's a faux cog zip file (~20 Mb uncompressed). I am going to look into some breakpoints now and see what might be going on

faux_mars_cog.zip

[AndrewAnnex]

Okay in the tile_exists function they are both rasterio CRSs but I wonder if this slowness is from the __eq__ in rasterio itself? If so it's exactly the same issue noted before in davenquinn/mars-tiler#2 but that issue never got moved up to rasterio itself.

in rasterio the eq is implemented in cython but does a lot of work:

https://github.com/rasterio/rasterio/blob/81ba49073cccef6384beb8d6c03eea7fb0b95e16/rasterio/crs.pyx#L923-L947

while writing this I just attempted this in ipython and sure enough:

In [3]: %timeit tms.rasterio_crs == tms.rasterio_crs
110 ms ± 829 μs per loop (mean ± std. dev. of 7 runs, 1 loop each)

I made rasterio/rasterio#3206 a few days ago that at least adds the IAU CRS to the to_authority method but maybe more work is need or there is something we could do to sidestep the equality and just compare strings...

[AndrewAnnex]

ahhh so comparing the wkt properties would be faster, the question I guess is does that make sense to switch over to within rio-tiler or should changes to rasterio be made.

tmsa = ptile.defaults.planetary_tms.get('MarsGeographicSphere')

In [36]: tmsa.rasterio_geographic_crs == tmsa.rasterio_crs
Out[36]: True

In [37]: %timeit tmsa.rasterio_geographic_crs == tmsa.rasterio_crs
21 ms ± 31.4 μs per loop (mean ± std. dev. of 7 runs, 10 loops each)

In [38]: %timeit tmsa.rasterio_geographic_crs == tmsa.rasterio_crs
15.5 ms ± 2.52 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)

In [39]: %timeit tmsa.rasterio_geographic_crs == tmsa.rasterio_crs
14 ms ± 7.86 μs per loop (mean ± std. dev. of 7 runs, 100 loops each)

In [40]: %timeit tmsa.rasterio_geographic_crs.to_wkt() == tmsa.rasterio_crs.to_wkt()
8.91 μs ± 997 ns per loop (mean ± std. dev. of 7 runs, 100,000 loops each)

In [41]: %timeit tmsa.rasterio_geographic_crs.to_wkt == tmsa.rasterio_crs.wkt
76.3 ns ± 8.41 ns per loop (mean ± std. dev. of 7 runs, 10,000,000 loops each)

In [42]: tmsa.rasterio_crs
Out[42]: CRS.from_wkt('GEOGCS["Mars (2015) - Sphere XY",DATUM["Mars (2015) - Sphere",SPHEROID["Mars (2015) - Sphere",3396190,0]],PRIMEM["Reference Meridian",0],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AXIS["Longitude",EAST],AXIS["Latitude",NORTH]]')

In [43]: tmsa.rasterio_geographic_crs
Out[43]: CRS.from_wkt('GEOGCS["Mars (2015) - Sphere XY",DATUM["Mars (2015) - Sphere",SPHEROID["Mars (2015) - Sphere",3396190,0]],PRIMEM["Reference Meridian",0],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AXIS["Longitude",EAST],AXIS["Latitude",NORTH]]')

In [50]: e=rio.CRS.from_epsg(4326)

In [51]: %timeit e == e
71.4 ns ± 7.7 ns per loop (mean ± std. dev. of 7 runs, 10,000,000 loops each)

[vincentsarago]

🤯 I never though that the __eq__ method could be so costly

[AndrewAnnex]

yikes, well I'll post an issue (if not done already or discussed before) to rasterio and may need to go work on some other things for the rest of the day,

[AndrewAnnex]

replacing all crs equality or inequality checks with comparing the .wkt property in reader.py, base.py, rasterio.py results in a pretty significant speedup for a single tile read:

(marspre) aannex@daw:/scratch/aannex$ python profile_rio_tiler.py
         217 function calls (212 primitive calls) in 0.003 seconds

   Ordered by: cumulative time

   ncalls  tottime  percall  cumtime  percall filename:lineno(function)
        1    0.000    0.000    0.003    0.003 /home/aannex/projects/rio-tiler/rio_tiler/io/rasterio.py:251(tile)
        1    0.000    0.000    0.001    0.001 /home/aannex/projects/rio-tiler/rio_tiler/io/rasterio.py:299(part)
        1    0.000    0.000    0.001    0.001 /home/aannex/projects/rio-tiler/rio_tiler/reader.py:296(part)
        1    0.000    0.000    0.001    0.001 /home/aannex/projects/rio-tiler/rio_tiler/io/base.py:157(tile_exists)
        1    0.001    0.001    0.001    0.001 /home/aannex/projects/rio-tiler/rio_tiler/reader.py:86(read)
        1    0.000    0.000    0.001    0.001 /home/aannex/micromamba/envs/marspre/lib/python3.12/functools.py:977(__get__)
        1    0.000    0.000    0.001    0.001 /home/aannex/micromamba/envs/marspre/lib/python3.12/site-packages/morecantile/models.py:560(rasterio_crs)
        1    0.001    0.001    0.001    0.001 /home/aannex/micromamba/envs/marspre/lib/python3.12/site-packages/morecantile/utils.py:130(to_rasterio_crs)

What probably makes sense is to add a benchmark to rio-tiler that exhibits this issue to start, and then wait for a newer rasterio (which wouldn't need any changes here)

[vincentsarago]

we do have 2 benchmark in place

• https://github.com/cogeotiff/rio-tiler/blob/main/tests/benchmarks/benchmarks.py which test the tile() method for all the datetype/nodata-mask-alpha combination
• https://github.com/cogeotiff/rio-tiler/blob/main/tests/benchmarks/requests.py which tests the number of GET request for a Tile Request

I'm not quite sure we need another set of benchmark, I believe the changes in rasterio will be reflected in the current benchmark https://cogeotiff.github.io/rio-tiler/advanced/benchmark.html

[vincentsarago]

@AndrewAnnex It seems that the real improvement will be with the CRS.equals method but this will be in rasterio 1.5

At the same time, I'm about to release a new major version of rio-tiler and I don't want to wait on rasterio releases. As soon as rasterio 1.5 is released we will adapt the .equals method and release a new major version of rio-tiler

[AndrewAnnex]

@vincentsarago sounds great to me, thanks so much for working with me on this by the way!