Repository navigation
transform() function rotates shapes incorrectly -> incorrect extend #353
Description
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
Activity
LucaMarconato commentedon Oct 9, 2023
Quick question @timtreis. Running the code produces this:
I think I need to use one of the new branches of
spatialdata-plotthat will make it into the next release. Which one?LucaMarconato commentedon Oct 9, 2023
I guess scverse/spatialdata-plot#162, gonna try now.
timtreis commentedon Oct 9, 2023
Yes, this one. Can you reproduce?
LucaMarconato commentedon Oct 9, 2023
Yes thanks. This is the cause of the bug:
timtreis commentedon Oct 9, 2023
Ahhh yeah, makes sense. That's why it scales with
rLucaMarconato commentedon Oct 9, 2023
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 commentedon Oct 9, 2023
Fixed
Gonna make a PR to the
get_extent()PR, I will also fix the failing tests.LucaMarconato commentedon Oct 9, 2023
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:
spatialdata-plotspatialdata-plotto add tests for the cases above; in fact in the tests inspatialdataI will just check that the coordinates are correct; I am making the plots now just to see what's going onLucaMarconato commentedon Oct 9, 2023
All correctly computed. I have opened the aforementioned issues in
spatialdata-plotand I have made a PR to merge against theget_extent()PR.LucaMarconato commentedon Oct 10, 2023
Fixed by #373 and merged via #318