I have some spatial data that has all of its coordinates as lat/lon pairs (with about 10 digits decimal precision), it's stored in a database as WGS84 data. Some of the data is represented as polygons which are the resulting union of some smaller polygons whose boundaries are stored. I also have a number of points from which I build line segments (just 2 points in each segment) which I use later for intersection tests with the polygons.
I'm using a SpatialIndex to improve my queries so I insert the envelopes of all polygons in a tree (tested with both QuadTree and STRtree). Then, I connect two points into a line segment and I'm using its envelope to query the tree for possible intersections. The problem is that I get pretty much all the polygons as a result which is clearly wrong. To give you some idea about the real scale of my data, I have about 100 polygons that cover the whole North America, each line covers a very very small part of a single polygon. Ideally, I would expect no more than 2 polygons as a result.
I'm using JTS to do this calculation and I'm aware that it's not really suited for spherical data so can you suggest me another library/tool to achieve the desired behaviour or possible a workaround (for example, projecting before using JTS)?
EDIT: I'm using postgresql + postgis extension, storing the all spatial data as geometry in 4326.
Best Answer
Fixed it.
The problem was a weird bug that was reversing latitude/longitude at some place.
Anyway, in general JTS can be used for calculations on latitude/longitude data. Although JTS doesn't try to project the coordinates before applying any spatial functions, all of them work fine (including the index trees) for smaller areas... I haven't really tested this with large ones but according to my findings - the larger the area, the less accurate the measures will be. This makes sense taking into account the distortion caused by the usage of spherical coordinates as if they were Cartesian.
Note that PostGIS has a slightly different behaviour. Whenever the coordinates are stored as WG84 and you use a function that is not yet natively implemented for spherical data the following happens - a suitable projection is automatically chosen and applied by PostGIS. Then the function is executed on the resulting projection rather than performing calculations directly on the sphere.