I think this has been discussed before, but I just wanted to double-check... am I correct in thinking that operations like "distance" do not produce useful results when working in a lat/lon coordinate system?
I.E. they produce results in "degrees" but since a degree is not a constant distance along the surface of the earth, it can't be compared usefully with any other distance. I know PostGIS has a "DistanceSphere" function that will calculate distance (in meters) across the surface of a sphere between two lat/lon points, is there anything similar in JTS? Or another library I might use? Thanks, Jeff _______________________________________________ jts-devel mailing list [hidden email] http://lists.refractions.net/mailman/listinfo/jts-devel |
You're correct. Operations in JTS assume a uniform Cartesian grid.
This makes metric operations such as distance problematic. There's nothing currently in JTS to address this. I have had some thoughts about starting work on supporting ellipsoidal models for distance calculations. Haven't got very far with this, though. If you want to help out it would be helpful to know what operations you want to perform in the ellipsoidal space (e.g. distance between points, length/area, buffering, etc). It's going to be a big job to retrofit all existing operations, but the first place to start is to get some use cases identified, and then perhaps to provide some simple prototype implementations which can be used to test out the design. Martin Jeff Adams wrote: > I think this has been discussed before, but I just wanted to > double-check... am I correct in thinking that operations like > "distance" do not produce useful results when working in a lat/lon > coordinate system? > > I.E. they produce results in "degrees" but since a degree is not a > constant distance along the surface of the earth, it can't be compared > usefully with any other distance. > > I know PostGIS has a "DistanceSphere" function that will calculate > distance (in meters) across the surface of a sphere between two > lat/lon points, is there anything similar in JTS? Or another library > I might use? > > Thanks, > Jeff > ------------------------------------------------------------------------ > > _______________________________________________ > jts-devel mailing list > [hidden email] > http://lists.refractions.net/mailman/listinfo/jts-devel > -- Martin Davis Senior Technical Architect Refractions Research, Inc. (250) 383-3022 _______________________________________________ jts-devel mailing list [hidden email] http://lists.refractions.net/mailman/listinfo/jts-devel |
I would consider projecting the data first onto some kind of Cartesian grid, chosen depending on what you need to do (e.g. to measure distances, choose an equidistant projection, like an azimuthal). Any reason this is unfeasible?
-rory
On Thu, Feb 5, 2009 at 8:51 AM, Martin Davis <[hidden email]> wrote: You're correct. Operations in JTS assume a uniform Cartesian grid. This makes metric operations such as distance problematic. _______________________________________________ jts-devel mailing list [hidden email] http://lists.refractions.net/mailman/listinfo/jts-devel |
In reply to this post by Martin Davis
More thoughts on this...
I can imagine a JTS enhancement which would offer the following: - ability to define a Datum (Cartesian, spherical, ellipsoidal) - ability to define a mapping from the Geometry coordinate system onto the Datum coordinate system (Cartesian or spherical). This would support using existing Map Projection libraries to perform this transformation. - this information would be carried as parameters on the GeometryFactory - functions to compute distance, area, and azimuth on the Datum - a way of injecting this information into operations which required it - updates to operations to use the Datum-based metrics in internal computation I'm not sure about the interface to Map Projection libraries. That would imply coordinate conversion on the fly, which is perhaps too expensive. Maybe it's better to leave this as an external concern, and just support 3 datums (Cartesian, spherical, ellipsoidal) using standard spherical coordinates (lon/lat). (Although if this idea is taken to the extreme, you wind up with Rory's suggestion of simply transforming in to a Cartesian space and computing there. I'm interested in see what people think about this.) Martin Davis wrote: > You're correct. Operations in JTS assume a uniform Cartesian grid. > This makes metric operations such as distance problematic. > > There's nothing currently in JTS to address this. I have had some > thoughts about starting work on supporting ellipsoidal models for > distance calculations. Haven't got very far with this, though. If > you want to help out it would be helpful to know what operations you > want to perform in the ellipsoidal space (e.g. distance between > points, length/area, buffering, etc). It's going to be a big job to > retrofit all existing operations, but the first place to start is to > get some use cases identified, and then perhaps to provide some simple > prototype implementations which can be used to test out the design. > > Martin > > Jeff Adams wrote: >> I think this has been discussed before, but I just wanted to >> double-check... am I correct in thinking that operations like >> "distance" do not produce useful results when working in a lat/lon >> coordinate system? >> >> I.E. they produce results in "degrees" but since a degree is not a >> constant distance along the surface of the earth, it can't be >> compared usefully with any other distance. >> >> I know PostGIS has a "DistanceSphere" function that will calculate >> distance (in meters) across the surface of a sphere between two >> lat/lon points, is there anything similar in JTS? Or another library >> I might use? >> >> Thanks, >> Jeff >> ------------------------------------------------------------------------ >> >> _______________________________________________ >> jts-devel mailing list >> [hidden email] >> http://lists.refractions.net/mailman/listinfo/jts-devel >> > -- Martin Davis Senior Technical Architect Refractions Research, Inc. (250) 383-3022 _______________________________________________ jts-devel mailing list [hidden email] http://lists.refractions.net/mailman/listinfo/jts-devel |
In reply to this post by Rory Plaire
Is the issue that for large geometries there is no suitable projection
which reduces distortion to an acceptable level? E.g. if you want to measure the coastline of Canada there's no projection which will give you the "correct" answer. I know that people often complain about this with reference to buffering on the PostGIS list. Buffers of large objects require a distance metric which varies depending on location. Rory Plaire wrote: > I would consider projecting the data first onto some kind of Cartesian > grid, chosen depending on what you need to do (e.g. to measure > distances, choose an equidistant projection, like an azimuthal). Any > reason this is unfeasible? > > -rory > > On Thu, Feb 5, 2009 at 8:51 AM, Martin Davis <[hidden email] > <mailto:[hidden email]>> wrote: > > You're correct. Operations in JTS assume a uniform Cartesian > grid. This makes metric operations such as distance problematic. > > There's nothing currently in JTS to address this. I have had some > thoughts about starting work on supporting ellipsoidal models for > distance calculations. Haven't got very far with this, though. > If you want to help out it would be helpful to know what > operations you want to perform in the ellipsoidal space (e.g. > distance between points, length/area, buffering, etc). It's going > to be a big job to retrofit all existing operations, but the first > place to start is to get some use cases identified, and then > perhaps to provide some simple prototype implementations which can > be used to test out the design. > > Martin > > Jeff Adams wrote: > > I think this has been discussed before, but I just wanted to > double-check... am I correct in thinking that operations like > "distance" do not produce useful results when working in a > lat/lon coordinate system? > > I.E. they produce results in "degrees" but since a degree is > not a constant distance along the surface of the earth, it > can't be compared usefully with any other distance. > > I know PostGIS has a "DistanceSphere" function that will > calculate distance (in meters) across the surface of a sphere > between two lat/lon points, is there anything similar in JTS? > Or another library I might use? > > Thanks, > Jeff > ------------------------------------------------------------------------ > > _______________________________________________ > jts-devel mailing list > [hidden email] > <mailto:[hidden email]> > http://lists.refractions.net/mailman/listinfo/jts-devel > > > > -- > Martin Davis > Senior Technical Architect > Refractions Research, Inc. > (250) 383-3022 > > _______________________________________________ > jts-devel mailing list > [hidden email] > <mailto:[hidden email]> > http://lists.refractions.net/mailman/listinfo/jts-devel > > > ------------------------------------------------------------------------ > > _______________________________________________ > jts-devel mailing list > [hidden email] > http://lists.refractions.net/mailman/listinfo/jts-devel > -- Martin Davis Senior Technical Architect Refractions Research, Inc. (250) 383-3022 _______________________________________________ jts-devel mailing list [hidden email] http://lists.refractions.net/mailman/listinfo/jts-devel |
In reply to this post by Martin Davis
In SharpMap, we use on-the-fly reprojection a lot, and the computational cost of doing this is a small percentage of application CPU consumption, even with millions of coordinates. I'd suspect it would be a fairly efficient route, computationally, if the error which can result from projection can be tolerated.
On the java side of things, I've used JTS with GeoTools' OGC-compliant reprojection code quite successfully for similar tasks as this - projecting from lat/lon to mercator, performing some measurements, and unprojecting back to lat/lon.
-r
On Thu, Feb 5, 2009 at 9:21 AM, Martin Davis <[hidden email]> wrote: More thoughts on this... _______________________________________________ jts-devel mailing list [hidden email] http://lists.refractions.net/mailman/listinfo/jts-devel |
Thanks, Rory.
This adds to my impression is that there isn't a big driver out there to add the ability to compute on the ellipsoid to JTS. It would be nifty to do - but perhaps not worth the time and effort? Any comments, Jeff? Rory Plaire wrote: > In SharpMap, we use on-the-fly reprojection a lot, and the > computational cost of doing this is a small percentage of application > CPU consumption, even with millions of coordinates. I'd suspect it > would be a fairly efficient route, computationally, if the error which > can result from projection can be tolerated. > On the java side of things, I've used JTS with GeoTools' > OGC-compliant reprojection code quite successfully for similar tasks > as this - projecting from lat/lon to mercator, performing some > measurements, and unprojecting back to lat/lon. > > -r > On Thu, Feb 5, 2009 at 9:21 AM, Martin Davis <[hidden email] > <mailto:[hidden email]>> wrote: > > More thoughts on this... > > I can imagine a JTS enhancement which would offer the following: > - ability to define a Datum (Cartesian, spherical, ellipsoidal) > - ability to define a mapping from the Geometry coordinate system > onto the Datum coordinate system (Cartesian or spherical). This > would support using existing Map Projection libraries to perform > this transformation. > - this information would be carried as parameters on the > GeometryFactory > - functions to compute distance, area, and azimuth on the Datum > - a way of injecting this information into operations which > required it > - updates to operations to use the Datum-based metrics in internal > computation > > I'm not sure about the interface to Map Projection libraries. > That would imply coordinate conversion on the fly, which is > perhaps too expensive. Maybe it's better to leave this as an > external concern, and just support 3 datums (Cartesian, spherical, > ellipsoidal) using standard spherical coordinates (lon/lat). > (Although if this idea is taken to the extreme, you wind up with > Rory's suggestion of simply transforming in to a Cartesian space > and computing there. I'm interested in see what people think > about this.) > > > > Martin Davis wrote: > > You're correct. Operations in JTS assume a uniform Cartesian > grid. This makes metric operations such as distance problematic. > > There's nothing currently in JTS to address this. I have had > some thoughts about starting work on supporting ellipsoidal > models for distance calculations. Haven't got very far with > this, though. If you want to help out it would be helpful to > know what operations you want to perform in the ellipsoidal > space (e.g. distance between points, length/area, buffering, > etc). It's going to be a big job to retrofit all existing > operations, but the first place to start is to get some use > cases identified, and then perhaps to provide some simple > prototype implementations which can be used to test out the > design. > > Martin > > Jeff Adams wrote: > > I think this has been discussed before, but I just wanted > to double-check... am I correct in thinking that > operations like "distance" do not produce useful results > when working in a lat/lon coordinate system? > > I.E. they produce results in "degrees" but since a degree > is not a constant distance along the surface of the earth, > it can't be compared usefully with any other distance. > > I know PostGIS has a "DistanceSphere" function that will > calculate distance (in meters) across the surface of a > sphere between two lat/lon points, is there anything > similar in JTS? Or another library I might use? > > Thanks, > Jeff > ------------------------------------------------------------------------ > > _______________________________________________ > jts-devel mailing list > [hidden email] > <mailto:[hidden email]> > http://lists.refractions.net/mailman/listinfo/jts-devel > > > > > -- > Martin Davis > Senior Technical Architect > Refractions Research, Inc. > (250) 383-3022 > > _______________________________________________ > jts-devel mailing list > [hidden email] > <mailto:[hidden email]> > http://lists.refractions.net/mailman/listinfo/jts-devel > > > ------------------------------------------------------------------------ > > _______________________________________________ > jts-devel mailing list > [hidden email] > http://lists.refractions.net/mailman/listinfo/jts-devel > -- Martin Davis Senior Technical Architect Refractions Research, Inc. (250) 383-3022 _______________________________________________ jts-devel mailing list [hidden email] http://lists.refractions.net/mailman/listinfo/jts-devel |
In reply to this post by Jeff Adams-4
Hi.
I've discussed this issue before, and have found a couple of approximate ways of doing distance calculations:- 1. Calculate a Vincenty Distance (equivalent to the DistanceSphere) - I can let you have the code. This one is pretty accurate. 2. I have a reprojection algorythm that's more approximate, but will work for any size of geometry. You will just have to accept that you can't draw an sphere on a piece of paper and live with the approximation. I use X = 111230 * cosine(radians(latitude)) m /degree Longitude Y = 111230 m /degree Latitude to reproject based on the lon /lat of the geometry's centroid. The driver for this requirement is when working with GPS systems internationally. You can get into a real mess with regional datums. Best to stick to lat /long. _______________________________________________ jts-devel mailing list [hidden email] http://lists.refractions.net/mailman/listinfo/jts-devel |
In reply to this post by Martin Davis
I'd point out that there is some existing code in Geotools that can
help out with this problem. They've got a class called GeodeticCalculator or something like that. It can return distances on the surface of the Ellipsoid between two points. Martin is the King of JTS, but I'll throw in my 2 cents: I think JTS is popular and successful because it does one thing (represent 2D goemetries on a plane) and it does it well. JTS doesn't mess with curves, it doesn't mess with 3D, and I don't think it should mess with geometries on the ellipsoid either. Package that stuff in a separate library that could interface with JTS at some level. Think about Coordinate objects in JTS for a minute. The z ordinate is there, but it isn't really used. In some sense that makes JTS a 2.5 D library. I think any good library that deals with geometry calculations on the ellipsoid should be set up for true 3D calculations. I'd also point out that map projections don't solve all problems, as mentioned above. I think geometry calculations on the ellipsoid is a unique problem domain that stands all on its own. It deserves its own library. (Or clearly segregated packages at least.) :] The Sunburned Surveyor On Thu, Feb 5, 2009 at 9:37 AM, Martin Davis <[hidden email]> wrote: > Thanks, Rory. > This adds to my impression is that there isn't a big driver out there to add > the ability to compute on the ellipsoid to JTS. It would be nifty to do - > but perhaps not worth the time and effort? > > Any comments, Jeff? > > Rory Plaire wrote: >> >> In SharpMap, we use on-the-fly reprojection a lot, and the computational >> cost of doing this is a small percentage of application CPU consumption, >> even with millions of coordinates. I'd suspect it would be a fairly >> efficient route, computationally, if the error which can result from >> projection can be tolerated. >> On the java side of things, I've used JTS with GeoTools' OGC-compliant >> reprojection code quite successfully for similar tasks as this - projecting >> from lat/lon to mercator, performing some measurements, and unprojecting >> back to lat/lon. >> -r >> On Thu, Feb 5, 2009 at 9:21 AM, Martin Davis <[hidden email] >> <mailto:[hidden email]>> wrote: >> >> More thoughts on this... >> >> I can imagine a JTS enhancement which would offer the following: >> - ability to define a Datum (Cartesian, spherical, ellipsoidal) >> - ability to define a mapping from the Geometry coordinate system >> onto the Datum coordinate system (Cartesian or spherical). This >> would support using existing Map Projection libraries to perform >> this transformation. >> - this information would be carried as parameters on the >> GeometryFactory >> - functions to compute distance, area, and azimuth on the Datum >> - a way of injecting this information into operations which >> required it >> - updates to operations to use the Datum-based metrics in internal >> computation >> >> I'm not sure about the interface to Map Projection libraries. >> That would imply coordinate conversion on the fly, which is >> perhaps too expensive. Maybe it's better to leave this as an >> external concern, and just support 3 datums (Cartesian, spherical, >> ellipsoidal) using standard spherical coordinates (lon/lat). >> (Although if this idea is taken to the extreme, you wind up with >> Rory's suggestion of simply transforming in to a Cartesian space >> and computing there. I'm interested in see what people think >> about this.) >> >> >> >> Martin Davis wrote: >> >> You're correct. Operations in JTS assume a uniform Cartesian >> grid. This makes metric operations such as distance problematic. >> >> There's nothing currently in JTS to address this. I have had >> some thoughts about starting work on supporting ellipsoidal >> models for distance calculations. Haven't got very far with >> this, though. If you want to help out it would be helpful to >> know what operations you want to perform in the ellipsoidal >> space (e.g. distance between points, length/area, buffering, >> etc). It's going to be a big job to retrofit all existing >> operations, but the first place to start is to get some use >> cases identified, and then perhaps to provide some simple >> prototype implementations which can be used to test out the >> design. >> >> Martin >> >> Jeff Adams wrote: >> >> I think this has been discussed before, but I just wanted >> to double-check... am I correct in thinking that >> operations like "distance" do not produce useful results >> when working in a lat/lon coordinate system? >> >> I.E. they produce results in "degrees" but since a degree >> is not a constant distance along the surface of the earth, >> it can't be compared usefully with any other distance. >> >> I know PostGIS has a "DistanceSphere" function that will >> calculate distance (in meters) across the surface of a >> sphere between two lat/lon points, is there anything >> similar in JTS? Or another library I might use? >> >> Thanks, >> Jeff >> >> ------------------------------------------------------------------------ >> >> _______________________________________________ >> jts-devel mailing list >> [hidden email] >> <mailto:[hidden email]> >> http://lists.refractions.net/mailman/listinfo/jts-devel >> >> >> >> -- Martin Davis >> Senior Technical Architect >> Refractions Research, Inc. >> (250) 383-3022 >> >> _______________________________________________ >> jts-devel mailing list >> [hidden email] >> <mailto:[hidden email]> >> http://lists.refractions.net/mailman/listinfo/jts-devel >> >> >> ------------------------------------------------------------------------ >> >> _______________________________________________ >> jts-devel mailing list >> [hidden email] >> http://lists.refractions.net/mailman/listinfo/jts-devel >> > > -- > Martin Davis > Senior Technical Architect > Refractions Research, Inc. > (250) 383-3022 > > _______________________________________________ > jts-devel mailing list > [hidden email] > http://lists.refractions.net/mailman/listinfo/jts-devel > jts-devel mailing list [hidden email] http://lists.refractions.net/mailman/listinfo/jts-devel |
In reply to this post by Paul Uszak
JGS = Java Geodesy Suite :]
Now all we need to do is get Martin funding for the development of the new library. Is anyone on good terms with a US senator? If so, we might have a shot at getting the money put in the 900 billion dollar US stimulus package. :] SS On Thu, Feb 5, 2009 at 9:41 AM, Paul Uszak <[hidden email]> wrote: > Hi. > > I've discussed this issue before, and have found a couple of approximate ways > of doing distance calculations:- > > 1. Calculate a Vincenty Distance (equivalent to the DistanceSphere) - I can > let you have the code. This one is pretty accurate. > > 2. I have a reprojection algorythm that's more approximate, but will work for > any size of geometry. You will just have to accept that you can't draw an > sphere on a piece of paper and live with the approximation. I use > > X = 111230 * cosine(radians(latitude)) m /degree Longitude > Y = 111230 m /degree Latitude > > to reproject based on the lon /lat of the geometry's centroid. > > The driver for this requirement is when working with GPS systems > internationally. You can get into a real mess with regional datums. Best to > stick to lat /long. > _______________________________________________ > jts-devel mailing list > [hidden email] > http://lists.refractions.net/mailman/listinfo/jts-devel > jts-devel mailing list [hidden email] http://lists.refractions.net/mailman/listinfo/jts-devel |
Woohoo! Bring it on... (as a famous former US politician once said)
I'd only need about 10% of it... 8^) I mostly agree with you, SS. The only nagging doubt is sparked by the buffering issue I mentioned before. I don't see any way of handling that without having JTS itself use a distance/azimuth calculation which is ellipsoidal-aware. I'd be very happy to not have the actual distance code in JTS (although it might not hurt to bundle at least a sample version) - but I don't see any way around providing an interface which can be used internally. Sunburned Surveyor wrote: > JGS = Java Geodesy Suite :] > > Now all we need to do is get Martin funding for the development of the > new library. Is anyone on good terms with a US senator? If so, we > might have a shot at getting the money put in the 900 billion dollar > US stimulus package. :] > > SS > > On Thu, Feb 5, 2009 at 9:41 AM, Paul Uszak <[hidden email]> wrote: > >> Hi. >> >> I've discussed this issue before, and have found a couple of approximate ways >> of doing distance calculations:- >> >> 1. Calculate a Vincenty Distance (equivalent to the DistanceSphere) - I can >> let you have the code. This one is pretty accurate. >> >> 2. I have a reprojection algorythm that's more approximate, but will work for >> any size of geometry. You will just have to accept that you can't draw an >> sphere on a piece of paper and live with the approximation. I use >> >> X = 111230 * cosine(radians(latitude)) m /degree Longitude >> Y = 111230 m /degree Latitude >> >> to reproject based on the lon /lat of the geometry's centroid. >> >> The driver for this requirement is when working with GPS systems >> internationally. You can get into a real mess with regional datums. Best to >> stick to lat /long. >> _______________________________________________ >> jts-devel mailing list >> [hidden email] >> http://lists.refractions.net/mailman/listinfo/jts-devel >> >> > _______________________________________________ > jts-devel mailing list > [hidden email] > http://lists.refractions.net/mailman/listinfo/jts-devel > > -- Martin Davis Senior Technical Architect Refractions Research, Inc. (250) 383-3022 _______________________________________________ jts-devel mailing list [hidden email] http://lists.refractions.net/mailman/listinfo/jts-devel |
Another issue with ellipsoidal (geodetic) calculation is the whole date
line/poles wrapping issue. This is a major pain to deal with. I guess local projections can be used to deal with this as well, though. Seems like it would be nice to have a library that would automatically choose the best projection, do the transformation, compute the result, and untransform... -- Martin Davis Senior Technical Architect Refractions Research, Inc. (250) 383-3022 _______________________________________________ jts-devel mailing list [hidden email] http://lists.refractions.net/mailman/listinfo/jts-devel |
Wow, what a conversation I started!
So here's my use case: I have a web app that is displaying polygons and points (different layers) on a "North America" scale (although of course you can zoom in to any particular area of interest). All the data is in "Web Mercator" projection. 1) I need to allow the user to choose a location and "search by radius (in miles)" for features in one of the point or polygon layers. 2) I need to find the nearest polygon to a point (although this one is less critical, as it's a hack to get around bad data and there's a chance we can get better data). Any reprojecting I will be doing will be on the fly, as the underlying data has to remain in web mercator. Since we're talking a web app, I'm a little concerned about reprojecting a whole lot of polygons on the fly just to check distance, I'm concerned that performance will be a problem (I.E. it won't be very responsive). Some of the polygons follow county boundaries and stuff so they can have a lot of points (although I could simplify them first if it were faster). Paul, I would love to get your code for "Vicinity Distance". In the interests of full disclosure however I must admit this is paid work for a commercial customer, so I can understand if you do not wish to give me a copy of the code. Thanks, Jeff On Thu, Feb 5, 2009 at 1:08 PM, Martin Davis <[hidden email]> wrote: Another issue with ellipsoidal (geodetic) calculation is the whole date line/poles wrapping issue. This is a major pain to deal with. I guess local projections can be used to deal with this as well, though. _______________________________________________ jts-devel mailing list [hidden email] http://lists.refractions.net/mailman/listinfo/jts-devel |
Hmmm... isn't your underlying data stored in a some kind of datastore?
So you will need to perform your search using the capabilities of that datastore? But I can see how you might be retrieving a set of candidate polygons (based on an intersection with a candidate bounding box for the search area), and then wanting to filter the results based on true distance. I think you're going to have to compute a suitable bounding box for your initial search in web mercator. This will require some ellipsoidal calculation, I think. To implement the secondary true distance filter you have two options - reproject first, or use the (non-existant) ellipsoidal JTS distance computation. It's not clear to me which would be more performant. Reprojection costs cycles, but computing ellipsoidal distance takes (a lot?) more cycles than Euclidean distance. It might be a toss-up about which is faster. You can easily try it - it shouldn't be too hard to code up a Distance class using the Vincenty distance. Jeff Adams wrote: > Wow, what a conversation I started! > > So here's my use case: > > I have a web app that is displaying polygons and points (different > layers) on a "North America" scale (although of course you can zoom in > to any particular area of interest). All the data is in "Web > Mercator" projection. > > 1) I need to allow the user to choose a location and "search by radius > (in miles)" for features in one of the point or polygon layers. > 2) I need to find the nearest polygon to a point (although this one is > less critical, as it's a hack to get around bad data and there's a > chance we can get better data). > > Any reprojecting I will be doing will be on the fly, as the underlying > data has to remain in web mercator. Since we're talking a web app, > I'm a little concerned about reprojecting a whole lot of polygons on > the fly just to check distance, I'm concerned that performance will be > a problem (I.E. it won't be very responsive). Some of the polygons > follow county boundaries and stuff so they can have a lot of points > (although I could simplify them first if it were faster). > > Paul, I would love to get your code for "Vicinity Distance". In the > interests of full disclosure however I must admit this is paid work > for a commercial customer, so I can understand if you do not wish to > give me a copy of the code. > > Thanks, > Jeff > > > On Thu, Feb 5, 2009 at 1:08 PM, Martin Davis <[hidden email] > <mailto:[hidden email]>> wrote: > > Another issue with ellipsoidal (geodetic) calculation is the whole > date line/poles wrapping issue. This is a major pain to deal > with. I guess local projections can be used to deal with this as > well, though. > > Seems like it would be nice to have a library that would > automatically choose the best projection, do the transformation, > compute the result, and untransform... > > > > -- > Martin Davis > Senior Technical Architect > Refractions Research, Inc. > (250) 383-3022 > > _______________________________________________ > jts-devel mailing list > [hidden email] > <mailto:[hidden email]> > http://lists.refractions.net/mailman/listinfo/jts-devel > > > ------------------------------------------------------------------------ > > _______________________________________________ > jts-devel mailing list > [hidden email] > http://lists.refractions.net/mailman/listinfo/jts-devel > -- Martin Davis Senior Technical Architect Refractions Research, Inc. (250) 383-3022 _______________________________________________ jts-devel mailing list [hidden email] http://lists.refractions.net/mailman/listinfo/jts-devel |
In reply to this post by Jeff Adams-4
Check this out:
http://www.movable-type.co.uk/scripts/latlong-vincenty.html Should be easy to port the code. If you do this, can you report out on your experience? E.g. does this run fast, does it give you the results you're looking for, etc. Jeff Adams wrote: > Wow, what a conversation I started! > > So here's my use case: > > I have a web app that is displaying polygons and points (different > layers) on a "North America" scale (although of course you can zoom in > to any particular area of interest). All the data is in "Web > Mercator" projection. > > 1) I need to allow the user to choose a location and "search by radius > (in miles)" for features in one of the point or polygon layers. > 2) I need to find the nearest polygon to a point (although this one is > less critical, as it's a hack to get around bad data and there's a > chance we can get better data). > > Any reprojecting I will be doing will be on the fly, as the underlying > data has to remain in web mercator. Since we're talking a web app, > I'm a little concerned about reprojecting a whole lot of polygons on > the fly just to check distance, I'm concerned that performance will be > a problem (I.E. it won't be very responsive). Some of the polygons > follow county boundaries and stuff so they can have a lot of points > (although I could simplify them first if it were faster). > > Paul, I would love to get your code for "Vicinity Distance". In the > interests of full disclosure however I must admit this is paid work > for a commercial customer, so I can understand if you do not wish to > give me a copy of the code. > > Thanks, > Jeff > > > On Thu, Feb 5, 2009 at 1:08 PM, Martin Davis <[hidden email] > <mailto:[hidden email]>> wrote: > > Another issue with ellipsoidal (geodetic) calculation is the whole > date line/poles wrapping issue. This is a major pain to deal > with. I guess local projections can be used to deal with this as > well, though. > > Seems like it would be nice to have a library that would > automatically choose the best projection, do the transformation, > compute the result, and untransform... > > > > -- > Martin Davis > Senior Technical Architect > Refractions Research, Inc. > (250) 383-3022 > > _______________________________________________ > jts-devel mailing list > [hidden email] > <mailto:[hidden email]> > http://lists.refractions.net/mailman/listinfo/jts-devel > > > ------------------------------------------------------------------------ > > _______________________________________________ > jts-devel mailing list > [hidden email] > http://lists.refractions.net/mailman/listinfo/jts-devel > -- Martin Davis Senior Technical Architect Refractions Research, Inc. (250) 383-3022 _______________________________________________ jts-devel mailing list [hidden email] http://lists.refractions.net/mailman/listinfo/jts-devel |
Yes, my underlying data is being served up by ArcGIS Server, which lets me search by intersection with a bounding box, then I was planning on filtering. If I try porting the distance algorithm (probably not until next week) I'll let you know how it goes.
The other suggestion I got from a coworker was to reproject the center point of the circle into a Cartesian system, use it there to create an approximate circle (100 points around the edge, or maybe even 50 or 25 would probably be plenty circular enough for this type of search), reproject that circle back to the original coordinate system, and do an intersection search using that instead of using a bounding box. I think I'd have to figure out the correct projection to use so my circle didn't get messed up going back to web mercator, but that may be the most performant option since I'm reprojecting a relatively small number of points. That does only solve case 1 however, for case 2 I still need to calculate the shortest distance among various polygons. Jeff On Thu, Feb 5, 2009 at 2:40 PM, Martin Davis <[hidden email]> wrote: Check this out: _______________________________________________ jts-devel mailing list [hidden email] http://lists.refractions.net/mailman/listinfo/jts-devel |
In reply to this post by Jeff Adams-4
The following is code that I've carefully ripped off from the Interweb and
messed with. You should be able to understand /extract the bit you need. It's a scary amount of trig. and I have absolutely no idea how it works, but seems to be some sort of iterative process. If you feel guilty using it for free in your project - donate some money to charity ;-) I think that it works, but I do not have the resources to check it on a larger scale than the Scottish tests (about 20 km distance). Perhaps others in this list can help verify the code. I would appreciate that a lot... ======================================================= package com.klickngo.geometry; import org.apache.commons.lang.builder.*; public class LatLongCoordinate { /* These are directly accessible */ public double altitude; public double latitudeY; public double longitudeX; private LatLongCoordinate() { this(0.0, 0.0, 0.0); } public LatLongCoordinate(double longitudeX, double latitudeY) { this(longitudeX, latitudeY, 0.0); } public LatLongCoordinate(double longitudeX, double latitudeY, double altitude) throws IllegalArgumentException { if (longitudeX > 180.0 | longitudeX < -180.0) throw new IllegalArgumentException("Longitude out of +/- 180 bounds"); if (latitudeY > 90.0 | latitudeY < -90.0) throw new IllegalArgumentException("Latitude out of +/- 90 bounds"); if (altitude < 0.0) throw new IllegalArgumentException("Altitude below sea level"); this.longitudeX = longitudeX; this.latitudeY = latitudeY; this.altitude = altitude; } public LatLongCoordinate(LatLongCoordinate coordinate) { this(coordinate.longitudeX, coordinate.latitudeY, coordinate.altitude); } public void setCoordinate(LatLongCoordinate otherCoord) { this.longitudeX = otherCoord.longitudeX; this.latitudeY = otherCoord.latitudeY; this.altitude = otherCoord.altitude; } /********************************************************** * * Did some testing in Scotland (via Memory Map). * Two points 18 by 9 km grid squares apart, results were: * VD = 20.125 km, Pythagarus = 20.116 km, MM = 20 km * * Two points 1 by 1 km grid square:- * VD = 1.413 km, Pythagarus = 1.414 km, MM = 1.41 km * **********************************************************/ public double getVincentyDistance(LatLongCoordinate otherCoord) { double retVal = 0.0D; double a = 6378137, b = 6356752.3142, f = 1 / 298.257223563; double l = Math.toRadians(otherCoord.longitudeX) - Math.toRadians(this.longitudeX); double u1 = Math.atan((1 - f) * Math.tan(Math.toRadians(this.latitudeY))); double u2 = Math.atan((1 - f) * Math.tan(Math.toRadians(otherCoord.latitudeY))); double sinU1 = Math.sin(u1), cosU1 = Math.cos(u1); double sinU2 = Math.sin(u2), cosU2 = Math.cos(u2); double lambda = l, lambdaP = 2 * Math.PI; int iterLimit = 20; double sinLambda, cosLambda; double sinSigma = 0.0D; double cosSigma = 0.0D; double sigma = 0.0D; double alpha; double cosSqAlpha = 0.0D; double cosToSigmaM = 0.0D; double c; while (Math.abs(lambda - lambdaP) > 1e-12 && --iterLimit > 0) { sinLambda = Math.sin(lambda); cosLambda = Math.cos(lambda); sinSigma = Math.sqrt( (cosU2 * sinLambda) * (cosU2 * sinLambda) + (cosU1 * sinU2 - sinU1 * cosU2 * cosLambda) * (cosU1 * sinU2 - sinU1 * cosU2 * cosLambda)); cosSigma = sinU1 * sinU2 + cosU1 * cosU2 * cosLambda; sigma = Math.atan2(sinSigma, cosSigma); alpha = Math.asin(cosU1 * cosU2 * sinLambda / sinSigma); cosSqAlpha = Math.cos(alpha) * Math.cos(alpha); cosToSigmaM = cosSigma - 2 * sinU1 * sinU2 / cosSqAlpha; c = f / 16 * cosSqAlpha * (4 + f * (4 - 3 * cosSqAlpha)); lambdaP = lambda; lambda = l + (1 - c) * f * Math.sin(alpha) * (sigma + c * sinSigma * (cosToSigmaM + c * cosSigma * (- 1 + 2 * cosToSigmaM * cosToSigmaM))); } // IF no convergence if (iterLimit == 0) { retVal = Double.NaN; } else { double uSq = cosSqAlpha * (a * a - b * b) / (b * b); double a2 = 1 + uSq / 16384 * (4096 + uSq * (-768 + uSq * (320 - 175 * uSq))); double b2 = uSq / 1024 * (256 + uSq * (-128 + uSq * (74 - 47 * uSq))); double deltaSigma = b2 * sinSigma * (cosToSigmaM + b2 / 4 * (cosSigma * (-1 + 2 * cosToSigmaM * cosToSigmaM) - b2 / 6 * cosToSigmaM * (-3 + 4 * sinSigma * sinSigma) * (-3 + 4 * cosToSigmaM * cosToSigmaM))); retVal = b * a2 * (sigma - deltaSigma); } return retVal; } public String toString() { return new ToStringBuilder(this). append("Long.", this.longitudeX). append("Lat.", this.latitudeY). append("Alt.", this.altitude). toString(); } public boolean equals(Object obj) { if (obj == null) { return false; } if (obj == this) { return true; } if (obj.getClass() != getClass()) { return false; } LatLongCoordinate rhs = (LatLongCoordinate) obj; return new EqualsBuilder() .append(this.altitude, rhs.altitude) .append(this.longitudeX, rhs.longitudeX) .append(this.latitudeY, rhs.latitudeY) .isEquals(); } public int hashCode() { // you pick a hard-coded, randomly chosen, non-zero, odd number // ideally different for each class return new HashCodeBuilder(307, 31). append(this.altitude). append(this.longitudeX). append(this.latitudeY). toHashCode(); } /* These private accessors are only for Hibernate to use */ private double getAltitude() { return altitude; } private void setAltitude(double altitude) { this.altitude = altitude; } private double getLongitudeX() { return longitudeX; } private void setLongitudeX(double longitudeX) { this.longitudeX = longitudeX; } private double getLatitudeY() { return latitudeY; } private void setLatitudeY(double latitudeY) { this.latitudeY = latitudeY; } } ============================================ _______________________________________________ jts-devel mailing list [hidden email] http://lists.refractions.net/mailman/listinfo/jts-devel |
In reply to this post by Martin Davis
2009/2/6 Martin Davis <[hidden email]>:
> Check this out: > > http://www.movable-type.co.uk/scripts/latlong-vincenty.html > > Should be easy to port the code. > Just for general info, the Geoscience Australia web site also has useful technical info and online calculators for geodetic distances. E.g. http://www.ga.gov.au/geodesy/datums/distance.jsp I must admit though, my gut feeling is to side with SS's view that JTS would do best to stick to being a great 2D cartesian library. Michael _______________________________________________ jts-devel mailing list [hidden email] http://lists.refractions.net/mailman/listinfo/jts-devel |
In reply to this post by Martin Davis
(Off topic)
Hey Martin, On Thu, 2009-02-05 at 09:24 -0800, Martin Davis wrote: > Is the issue that for large geometries there is no suitable projection > which reduces distortion to an acceptable level? E.g. if you want to > measure the coastline of Canada there's no projection which will give > you the "correct" answer. This is a particularly unfortunate example since this very question (for England, not Canada) was the title of a landmark scientific paper which showed the question to be meaningless. Coastlines are irregular down to the plank length. Any linear 'length' we assign will merely reflect the resolution with which we choose to describe the coast. It seems that the best way to tackle such shapes is to consider them to exist in a different dimension from linear space. That dimension to use turns out to be fractional (it lies between 1 and 2) which gives rise to, and names, 'fractals'. But that's got nothing to do with the discussion, and the proposed work which is of interest for other reasons. cheers, --adrian _______________________________________________ jts-devel mailing list [hidden email] http://lists.refractions.net/mailman/listinfo/jts-devel |
2009/2/7 Adrian Custer wrote:
> (Off topic) > > This is a particularly unfortunate example since this very question (for > England, not Canada) was the title of a landmark scientific paper which > showed the question to be meaningless. > Even more off topic.. but it's interesting to see that over 40 years after its publication, this paper is still regularly cited. It hasn't yet undergone the final transition of a classic paper - to become so much part of our worldview that we take it for granted and don't even think to cite it. Apologies to the OP for this little highjack :) Michael _______________________________________________ jts-devel mailing list [hidden email] http://lists.refractions.net/mailman/listinfo/jts-devel |
Free forum by Nabble | Edit this page |