`transform()` function rotates shapes incorrectly -> incorrect extend

[timtreis]

MRE

import geopandas as gpd
import spatialdata
import spatialdata_plot
from shapely import MultiPolygon, Point, Polygon

points = []
for p in [[0, 0], [0, 1], [1, 1], [1, 0]]:
    points.append(Point(p))

gdf = gpd.GeoDataFrame(geometry=points)
gdf["radius"] = 0.05
circles_df = ShapesModel.parse(gdf)
sdata = spatialdata.SpatialData(shapes={"circles": circles_df})

theta = math.pi / 4
rotation = Affine(
[
    [math.cos(theta), -math.sin(theta), 0],
    [math.sin(theta), math.cos(theta), 0],
    [0, 0, 1],
],
input_axes=("x", "y"),
output_axes=("x", "y"),
)
set_transformation(element=sdata["circles"], transformation=rotation, to_coordinate_system="transformed")

sdata.pl.render_shapes().pl.show()

for cs in ["global", "transformed"]:
    print(cs, get_extent(sdata, coordinate_system=cs))

Produces the following plot:

We expect the axis of the transformed plot to also touch the circles. This is not the case since the extent is calculated to transformed {'y': (-0.07071067811865475, 1.48492424049175), 'x': (-0.7778174593052023, 0.7778174593052024)}.

However, the error doesn't come from the extent calculation, but from L305 in _compute_extent_in_coordinate_system:

Bounding box before and after transformation:

Since we're rotating by 45°, we'd expect them to be

[ 0,      0]
[−0.7071, 0.7071]
[ 0,      1.4142]
[ 0.7071, 0.7071]

(+/- the radius of the circle 0.05) which would result in a different bounding box. With larger radii, the inaccuracy becomes bigger

radius = 0.5

[LucaMarconato]

Quick question @timtreis. Running the code produces this:

I think I need to use one of the new branches of spatialdata-plot that will make it into the next release. Which one?

[LucaMarconato]

I guess https://github.com/scverse/spatialdata-plot/pull/162, gonna try now.

[timtreis]

Yes, this one. Can you reproduce?

[LucaMarconato]

Yes thanks. This is the cause of the bug:

[timtreis]

Ahhh yeah, makes sense. That's why it scales with r

[LucaMarconato]

The current approach works for raster data, the affected part is points, shapes and polygons. I will implement the exact calculation for the vector data, enabled by default, but I will also make it optional with a parameter of get_extent.

I suspect that the current implementation can be a lot faster since it doesn't need to transform the data, and for most cases actually one doesn't need the exact extent.

[LucaMarconato]

Fixed Gonna make a PR to the get_extent() PR, I will also fix the failing tests.

[LucaMarconato]

I am adding tests, not just for circles but also for points, polygons and multipolygons.

These are the plots that one would get with the old code (note that I have changed the positions of the circles in the "global" coordinate system). 4 circles 4 polygons 4 multipolygons 4 points

This is with the new code. 4 circles 4 polygons 4 multipolygons 4 points

Things work, except for points. I will:

• [x] check if get_extent() works for points, and if it does report the bug in spatialdata-plot
• [x] open an issue in spatialdata-plot to add tests for the cases above; in fact in the tests in spatialdata I will just check that the coordinates are correct; I am making the plots now just to see what's going on

[LucaMarconato]

All correctly computed. I have opened the aforementioned issues in spatialdata-plot and I have made a PR to merge against the get_extent() PR.

[LucaMarconato]

Fixed by https://github.com/scverse/spatialdata/pull/373 and merged via https://github.com/scverse/spatialdata/pull/318