Tweaked the intersects call to set the geometry srid to the srid of t…

[RemcoMeeuwissen]

Back with another pull request :smile:

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') to st_intersects(geom, st_setsrid('polygon'::geometry, st_srid(geom))). More specifically:

• The 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 column
• So that we can then update the srid of the polygon to this srid using the st_setsrid call, 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.

[vincentsarago]

would be nice if we could add tests for tables with geom not in EPSG:4326

[RemcoMeeuwissen]

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.

[vincentsarago]

@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

[RemcoMeeuwissen]

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 😃

[bitner]

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.

Hah, I was part of the group that dataset together!

[RemcoMeeuwissen]

Small world haha.

[bitner]

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 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.

[RemcoMeeuwissen]

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.

[bitner]

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.

[bitner]

@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.

[RemcoMeeuwissen]

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 😅

[bitner]

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.

[RemcoMeeuwissen]

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.

[RemcoMeeuwissen]

Didn't get around to finishing it up on Friday but this should hopefully do the trick

[RemcoMeeuwissen]

@bitner Poke 😅