[postgis-users] creating a spatial index for finding objects withing a certaindistance in miles based on lat/lon
Lowther, David W
dlowther at ou.edu
Fri Sep 12 12:09:27 PDT 2003
Daniel,
Here's a bit of code that may help...
Given a lat/ lon (dGNISLat and dGNISLon) and 2 dir / distance (miles) pairs
(szDir1 nDist1, szDir2, nDist2) this function will calc a new lat/lon. If
you do this once for each corner (or modify in some other way) you can build
your box... Its in VB, but nothing too complex to port.
-------------------
Global Const FtPerDegLat = 365224.46
Global Const PI = 3.14159265359
Public Function CalcLocation(dGNISLat As Double, dGNISLon As Double,
szDir1$, szDir2$, dOutLat As Double, dOutLon As Double, Optional nDist1 As
Double, Optional nDist2 As Double) As Integer
Dim dCurrentLon As Double, dCurrentLat As Double
Dim dLonChange As Double, dLatChange As Double
Dim dLatPart As Double, dLonPart As Double
If nDist1 = 0 And nDist2 = 0 Then
dOutLat = dGNISLat
dOutLon = dGNISLon
Exit Function
End If
Select Case UCase(szDir1)
Case "N"
dLatChange = (5280 * nDist1) / FtPerDegLat
dCurrentLat = dGNISLat + dLatChange
dCurrentLon = dGNISLon
Case "S"
dLatChange = (5280 * nDist1) / FtPerDegLat
dCurrentLat = dGNISLat - dLatChange
dCurrentLon = dGNISLon
Case "E"
dLonChange = (5280 * nDist1) / (FtPerDegLat * Cos((dGNISLat * PI
/ 180)))
dCurrentLon = dGNISLon + dLonChange
dCurrentLat = dGNISLat
Case "W"
dLonChange = (5280 * nDist1) / (FtPerDegLat * Cos((dGNISLat * PI
/ 180)))
dCurrentLon = dGNISLon - dLonChange
dCurrentLat = dGNISLat
Case "NW"
dLonChange = ((5280 * nDist1) / 1.414) / (FtPerDegLat *
Cos((dGNISLat * PI / 180)))
dLatChange = ((5280 * nDist1) / 1.414) / FtPerDegLat
dCurrentLat = dGNISLat + dLatChange
dCurrentLon = dGNISLon - dLonChange
Case "NE"
dLonChange = ((5280 * nDist1) / 1.414) / (FtPerDegLat *
Cos((dGNISLat * PI / 180)))
dLatChange = ((5280 * nDist1) / 1.414) / FtPerDegLat
dCurrentLat = dGNISLat + dLatChange
dCurrentLon = dGNISLon + dLonChange
Case "SW"
dLonChange = ((5280 * nDist1) / 1.414) / (FtPerDegLat *
Cos((dGNISLat * PI / 180)))
dLatChange = ((5280 * nDist1) / 1.414) / FtPerDegLat
dCurrentLat = dGNISLat - dLatChange
dCurrentLon = dGNISLon - dLonChange
Case "SE"
dLonChange = ((5280 * nDist1) / 1.414) / (FtPerDegLat *
Cos((dGNISLat * PI / 180)))
dLatChange = ((5280 * nDist1) / 1.414) / FtPerDegLat
dCurrentLat = dGNISLat - dLatChange
dCurrentLon = dGNISLon + dLonChange
Case Else
'//gonna handle this?
End Select
Select Case UCase(szDir2)
Case "N"
dLatChange = (5280 * nDist2) / FtPerDegLat
dCurrentLat = dCurrentLat + dLatChange
dCurrentLon = dCurrentLon
Case "S"
dLatChange = (5280 * nDist2) / FtPerDegLat
dCurrentLat = dCurrentLat - dLatChange
dCurrentLon = dCurrentLon
Case "E"
dLonChange = (5280 * nDist2) / (FtPerDegLat * Cos((dGNISLat * PI
/ 180)))
dCurrentLon = dCurrentLon + dLonChange
dCurrentLat = dCurrentLat
Case "W"
dLonChange = (5280 * nDist2) / (FtPerDegLat * Cos((dGNISLat * PI
/ 180)))
dCurrentLon = dCurrentLon - dLonChange
dCurrentLat = dCurrentLat
Case "NW"
dLonChange = ((5280 * nDist2) / 1.414) / (FtPerDegLat *
Cos((dGNISLat * PI / 180)))
dLatChange = ((5280 * nDist2) / 1.414) / FtPerDegLat
dCurrentLat = dCurrentLat + dLatChange
dCurrentLon = dCurrentLon - dLonChange
Case "NE"
dLonChange = ((5280 * nDist2) / 1.414) / (FtPerDegLat *
Cos((dGNISLat * PI / 180)))
dLatChange = ((5280 * nDist2) / 1.414) / FtPerDegLat
dCurrentLat = dCurrentLat + dLatChange
dCurrentLon = dCurrentLon + dLonChange
Case "SW"
dLonChange = ((5280 * nDist2) / 1.414) / (FtPerDegLat *
Cos((dGNISLat * PI / 180)))
dLatChange = ((5280 * nDist2) / 1.414) / FtPerDegLat
dCurrentLat = dCurrentLat - dLatChange
dCurrentLon = dCurrentLon - dLonChange
Case "SE"
dLonChange = ((5280 * nDist2) / 1.414) / (FtPerDegLat *
Cos((dGNISLat * PI / 180)))
dLatChange = ((5280 * nDist2) / 1.414) / FtPerDegLat
dCurrentLat = dCurrentLat - dLatChange
dCurrentLon = dCurrentLon + dLonChange
Case Else
'//gonna handle this?
End Select
dOutLat = dCurrentLat
dOutLon = dCurrentLon
End Function
David Lowther
Software Engineer
GEO Information Systems
University of Oklahoma
dlowther at ou.edu
(405) 325-3131
http://www.geo.ou.edu
> -----Original Message-----
> From: Daniel Ceregatti [mailto:vi at sh.nu]
> Sent: Friday, September 12, 2003 2:06 PM
> To: Hugh W. O'Brien
> Cc: postgis-users at postgis.refractions.net
> Subject: Re: [postgis-users] creating a spatial index for
> finding objects withing a certaindistance in miles based on lat/lon
>
>
> Hugh W. O'Brien wrote:
>
> > Since I only had to support North America, I picked 3
> postal codes in
> > Houston that formed the opposite and adjacent sides of a right
> > triangle. Then, using the geocoded postal code data my company
> > purchased, I identified the 3 coordinates that stood for my postal
> > codes. Then using the distance_spheriod(...) ( I think )
> function of
> > postgis, I determined a degrees/mile conversion for latitude and
> > longitude. Note, at the time I did this the distance_spheroid(...)
> > function of postgis had a bug ( it couldn't handle either
> the verticle
> > or horizontal calc because of a bug ). So I tweak the
> geocoded postal
> > code value ( decreased the last decimal point value of one of the
> > points by 1 to get around the bug ).
>
> Unfortunately I can't narrow my search area down to such a
> small area. Ultimately I want to have the entire world be
> searchable via this facility.
>
> It seems that I simply need to incorporate the "degrees per
> distance unit" in the query that creates the box3d
> dynamically. Perhaps all that's needed now is to work out the
> math and test it.
>
> Oh, and figure out how to get around boxes that traverse 90
> degree latitudes and 0/180 degree longitudes.
>
> >
> >
> > I wish I could give you more help like sample code and my actual
> > conversions and sample query but I no longer work for the company I
> > did this for ( been unemployed for 3+ month atm ).
>
> Sorry to hear that. :(
>
> >
> >
> > I hope this helps.
>
> I believe this will help me understand this whole thing
> better. I'm just surprised that this is so difficult to
> implement in postgres, seeing
> how easy it was to implement in oracle. You'd think that this
> type of application would be standard fare in such a package.
>
> Daniel Ceregatti
>
> >
> >
> > -- Hugh
> >
> > Daniel Ceregatti wrote:
> >
> >> Sigh...I'm still having some trouble wrapping my mind around all
> >> this. I followed the thread referenced in this post:
> >>
> >>
> http://postgis.refractions.net/pipermail/postgis-users/2002-December/
> >> 001905.html
> >>
> >>
> >> This post shows a "bastardized" query for doing what I need. Thing
> >> is, it seems that there are some assumptions made. For
> example, the
> >> "Miles per lat, Miles per lon" deal. Obviously this will vary over
> >> great distances. As I want to be able to query distances
> around the
> >> entire world, and across 90 lon/180 lat, I don't see how
> this will do
> >> what I need.
> >>
> >> Hugh, perhaps you can tell me how you ultimately did this?
> My thread
> >> starts here:
> >>
> >>
> http://postgis.refractions.net/pipermail/postg> is-users/2003-September
> >> /003087.html
> >>
> >>
> >> Thanks,
> >>
> >> Daniel Ceregatti
> >>
> >> Paul Ramsey wrote:
> >>
> >>
> >>
> >>> Daniel Ceregatti wrote:
> >>>
> >>>
> >>>
> >>>>> Daniel, lat/lon *are* degrees. Degrees are the units of
> a lat/lon
> >>>>> data set. You can pretend that lat/lon coordinates are
> planar for
> >>>>> the
> >>>>>
> >>>>
> >>>> I'm not too "geometrically inclined", hence my
> confusion. So what
> >>>> you're saying is that my existing index is already
> appropriate for
> >>>> the query?
> >>>>
> >>>
> >>> Yes, that is what I am saying. You just have to question
> your index
> >>> in the right way. Instead of asking "what are all the
> points within
> >>> 3 miles of point X?" you have to ask "what are all the
> points within
> >>> a bounding box of (lon1 lat1,lon2 lat2) and also within 3
> miles of
> >>> point X?" So the only "tricky" bit is constructing the
> bounding box
> >>> so that it is about 6 miles on a side but expressed in degrees of
> >>> latitude and longitude.
> >>>
> >>> P
> >>>
> >>>
> >>
> >
> >
>
>
> _______________________________________________
> postgis-users mailing list postgis-users at postgis.refractions.net
> http://postgis.refractions.net/mailman/listinfo/postgis-users
>
>
More information about the postgis-users
mailing list