Tweaked the intersects call to set the geometry srid to the srid of t… #34
Conversation
|
would be nice if we could add tests for tables with geom not in EPSG:4326 |
|
How about https://gisdata.mn.gov/dataset/bdry-mn-county-open-data-status? I've been searching for a while and came across it, the shapefile is in 26915 and it's free to use without any restrictions. |
|
@RemcoMeeuwissen we can also create something like https://github.com/developmentseed/tipg/blob/main/tests/fixtures/my_data.sql but with geometries in other CRS |
|
Yeah that's possible as well, but I figured that having a existing dataset makes things a little easier since it's not hard to convert the shapefile to SQL so I can do that as well. Just wanted to check that adding that specific one as a fixture is fine 😃 |
Hah, I was part of the group that dataset together! |
|
Small world haha. |
@RemcoMeeuwissen that dataset is quite enormous to include as a fixture. We could just take a small sample of rows from that, but I would not include the whole parcels dataset. |
tipg/filter/filters.py
Outdated
| "INTERSECTS": lambda f, a: Func( | ||
| "st_intersects", | ||
| f, | ||
| Func("st_setsrid", cast(a, "geometry"), Func("st_srid", f)), |
There was a problem hiding this comment.
st_setsrid does not actually transform the geometry to the new srid, to actually change the geometry srid you need to use st_transform(geometry, srid)
|
Yeah, so there's basically two options that I see, A) The user sends a polygon in the coordinate system of the table and then we use setsrid to let postgres know that this is the case (the approach I've used here) B) The user sends a polygon in a fixed coordinate system (e.g. 4326) and we transform it into the local coordinate system. I've chosen for A since there's less overhead/performance costs, but some documentation somewhere on this edge case would be needed to avoid confusion. |
Unless we add a parameter specifying the SRID, we always assume that arguments are in geographic coordinates. I think this follows the OGC spec, but I would need to double check - whatever we do here, we should make sure that it aligns with the OGC spec. There should be negligible overhead/performance impacts to transforming the incoming geometry since it is only transforming a single geometry - we do want to avoid running the transform on the table geometries as that would be incurred for every geometry in the table. |
|
@RemcoMeeuwissen so right now we do not support OGC Features API Part 2 (https://ogcapi.ogc.org/features/) so output for features is always in geographic coordinates. You can see here that we hard code the return to always be in 4326 for features. For vector tiles we do allow different tile matrix sets. For filters, there is a filter-crs parameter in the specification that we would need to add in order to implement the behavior from your option B. Note that in the spec, filter-crs takes a full url definition for the crs and not just a simple srid number. |
|
I see that that would also require the implementation of https://docs.ogc.org/is/18-058r1/18-058r1.html, which would take a decent amount of work to implement. Would you say that the transform of the geometry is a acceptable solution then (taking the spec into account)? Because that would be more feasible to implement in the short run 😅 |
|
So, sticking to the spec without implementing "Part 2", the only path that I see right now would be to still require that filter geometries must be in geographic coordinates, but transform them to the dataset coordinates as per your option B. I've added a ticket #35 for implementing "Part 2" so that we can make sure to include filter-crs support when we add that. |
|
I'll start work on option B with a fixed CRS then, including adding tests and so forth. Should hopefully have a updated PR ready by the end of the week. |
…and renamed the geom field in the minnesota fixture
…tem of the column
…ed the unneeded cast
|
Didn't get around to finishing it up on Friday but this should hopefully do the trick |
|
@bitner Poke 😅 |
Back with another pull request 😄
What I am changing
For our data we have a lot of tables which are in a projection other than 4326. This is working fine for the tiles, but we ran into issues trying to do a intersect filter with CQL. Because we can send the polygon in the correct projection but since there is no way to pass along a coordinate system the resulting query still assumes it's 4326, thus causing the database to throw a error since it can't intersect from two different coordinate systems.
How I did it
The intersect query has been changed from
st_intersects(geom, 'polygon')tost_intersects(geom, st_setsrid('polygon'::geometry, st_srid(geom))). More specifically:cast(a, "geometry")is to typecast the polygon being passed to a geometry type to keep setsrid happy.Func("st_srid", f)is to get the srid of the geom columnst_setsridcall, allowing intersect to use the same coordinate system on both sides when it's doing it's thing.Related Issues
Not so much a issue, but I assume other queries such as CONTAINS and WITHIN have the same issues so if this tweak is good I can look into implementing it there as well.
I've also not written any tests for this yet, but for this we first need a dataset in a different coordinate system available.