Level Extreme platform
Subscription
Corporate profile
Products & Services
Support
Legal
Français
Map points within a circle or distance -- Pythagoras --
Message
From
24/10/2011 01:25:32
 
 
To
23/10/2011 22:58:30
General information
Forum:
Visual FoxPro
Category:
Coding, syntax & commands
Environment versions
Visual FoxPro:
VFP 9 SP1
OS:
Windows XP SP2
Network:
Windows 2003 Server
Database:
MS SQL Server
Application:
Desktop
Miscellaneous
Thread ID:
01527151
Message ID:
01527154
Views:
64
>Or maybe the Haversine Formula? Which? In VFP, If I have two coordinates on a map with x/y for 0/0 like the examples below and a radius to draw a circle around that position (let's say 500 ft), how can I tell if any other x/y coordinates are within the radius? I know I need to check and see if the distance to the 2nd x/y coordinates is less than the radius, but it escapes me -- math in general tonight. I've tried checking within the radius, just returning the distance in feet, but just cannot get it right....
>
>None of these are accurate--the distance should be about 392 ft (a mish-mash of attempts from me, the web, etc):
>
>*circumference is 360 degrees times 60 minutes/degree times 1.852 km/minute = 40003.2 km. 
>*The implied radius is the circumference divided by 2 pi: R = 6367 km = 3956 mi 
>lon1=  1698101.25   &&1713881.3750
> lat1=   813831.69 && 833601.8125
>lon2 =    1698215.13 &&1697432.8750
>lat2 =    813458.25 && 810161.0625
> R = 3956 && 3963 && 3956
> dlon = lon2 - lon1
>  dlat = lat2 - lat1
>  a = (sin(dlat/2))^2 + cos(lat1) * cos(lat2) * (sin(dlon/2))^2
>  ? a
>  c = 2 * atn2(sqrt(a), sqrt(1-a)) 
>  ? c
>  d = R * c
>  ?[Example1: ]
>  ??d
>
>   
>*--is point within a radius
>radius = 500  &&(what would 500 ft be?)
>centerx =  1698101.25   &&1713881.3750
>centery =   813831.69 && 833601.8125
>x =    1698215.13 &&1697432.8750
>y =    813458.25 && 810161.0625
>
>?[Example2: ] 
>?? (x-centerx)^2 + (y - centery)^2 
>?[     Inside radius: ] 
>?? (x-centerx)^2 + (y - centery)^2 < radius^2
>
>
>
>*--Get distance between two points
>#DEFINE KMRAD 6377
>#DEFINE MIRAD 3956 && 3963  
>#DEFINE RADIANS DTOR
>
>
>tlKM = .F.
>
>lon1=  1698101.25   &&1713881.3750
>lat1=   813831.69 && 833601.8125
>lon2 =    1698215.13 &&1697432.8750
>lat2 =    813458.25 && 810161.0625
>
>    lnDistance=IIF(tlKM, KMRAD, MIRAD)*;
>     ACOS(COS(DTOR(90-lat1))*COS(DTOR(90-lat2))+;
>     SIN(DTOR(90-lat1))*SIN(DTOR(90-lat2))*;
>     COS(DTOR(lon1 - lon2)))
>
>?[Example3: ]
>??lnDistance  && *5280 && 5280 feet in a mile
>
>
>
>*--Select all within the radius  &&Courtesy Carl Karsten 
>#DEFINE cdbl
>#DEFINE _PI 3.1415926535893
>#DEFINE _StatuteMile 5280 &&' 5280 feet = 1 statute mile
>#DEFINE _NauticalMile 6076.11549 &&' 6076.11549 feet = 1 nautical mile
>#DEFINE _Seconds 60 && ' 60 seconds = 1 nautical mile
>lon1=  1698101.25   &&1713881.3750
>lat1=   813831.69 && 833601.8125
>lon2 =    1698215.13 &&1697432.8750
>lat2 =    813458.25 && 810161.0625
>
>? InRadius(lat1,lon1,500,lat2,lon2)
>return
>
>FUNCTION InRadius( tlat, tlong, tnRadius, tlat2, tlong2 )
>    LOCAL lnResults, lnSeconds, ;
>      lnLatRange, lnLonRange,;
>      lnLowLat, lnLowLon, ;
>      lnHighLat, lnHighLon
>      
>   *' 1 degree of latitude = 60 nautical miles or 1 minute = 1 nautical mile
>    lnLatRange = tnRadius / ((_NauticalMile / _StatuteMile) * _Seconds)
>
>    *  Longitude is a bit more complicated; 
>    * 1 degree of longitude at the equator = 60 nautical miles, 
>    * At the poles 1 degree = 0 nautical miles.
>    lnLonRange = tnRadius / (((COS(cdbl(tlat * _PI / 180)) * ;
>				_NauticalMile) / _StatuteMile) * _Seconds)
>
>    lnLowLat = tlat - lnLatRange
>    lnHighLat = tlat + lnLatRange
>    lnLowLon = tlong - lnLonRange
>    lnHighLon = tlong+ lnLonRange
>
>?[Example4: ]
>?? tlat2 <= lnHighLat;
>      AND tlat2 >= lnLowLat;
>      AND tlong2 >= lnLowLon;
>      AND tlong2 <= lnHighLon
>?[Example5: ] 
>?? getdistance2(lon1,lat1,lon2,lat2)
>
>RETURN
>
>      
>FUNCTION getdistance2  &&Courtesy Carl Karsten
>LPARAMETERS lon1,lat1,lon2,lat2
>
>*--Get distance between two points
>#DEFINE KMRAD 6377
>#DEFINE MIRAD 3956 && 3963  
>#DEFINE RADIANS DTOR
>
>
>tlKM = .F.
>
>lon1=  1698101.25   &&1713881.3750
>lat1=   813831.69 && 833601.8125
>lon2 =    1698215.13 &&1697432.8750
>lat2 =    813458.25 && 810161.0625
>
>    lnDistance=IIF(tlKM, KMRAD, MIRAD)*;
>     ACOS(COS(DTOR(90-lat1))*COS(DTOR(90-lat2))+;
>     SIN(DTOR(90-lat1))*SIN(DTOR(90-lat2))*;
>     COS(DTOR(lon1 - lon2)))
>
>?[Example6: ]
>??lnDistance  
>
>return lndistance
>
>
>
>I'm brain dead! I've mucked it all up.....

Long time no see! Great to see you back!

It looks like for sure the Haversine function is what you need. Here's some code I converted from a Java implementation on Wikipedia:
* FUNCTION GreatCircleDistance
LPARAMETERS ;
	tnLat1 ;
	, tnLon1 ;
	, tnLat2 ;
	, tnLon2

LOCAL ;
	lnDeltaLat ;
	, lnDeltaLon ;
	, lnA

* All 4 parameters are in degrees
* Return value is in kilometers

* Java implementation from Wikipedia:
*!*	   double haversine(double lat1, double lng1, double lat2, double lng2) {
*!*	        double dLat = Math.toRadians(lat2 - lat1);
*!*	        double dLon = Math.toRadians(lng2 - lng1);
*!*	        double a = Math.sin(dLat / 2) * Math.sin(dLat / 2) + Math.cos(Math.toRadians(lat1)) * 
*!*	           Math.cos(Math.toRadians(lat2)) * Math.sin(dLon / 2) * Math.sin(dLon / 2);
*!*	        return 2 * Math.atan2(Math.sqrt(a), Math.sqrt(1 - a)) * 6371; // 6371 = average earth radius

* Calculate in radians and kilometers
* There are 2 * PI() radians in a circle, one radian is 180 / PI() degrees
lnDeltaLat = ( tnLat2 - tnLat1 ) * ( PI( ) / 180 )
lnDeltaLon = ( tnLon2 - tnLon1 ) * ( PI( ) / 180 )

lnA = ( SIN( lnDeltaLat / 2 ) ) ^ 2 ;
	+ COS( tnLat1 * PI( ) / 180 ) * COS( tnLat2 * PI( ) / 180 ) * ( SIN( lnDeltaLon / 2 ) ) ^ 2

RETURN 2 * Atan2( SQRT ( lnA ), SQRT( 1 - lnA ) ) * 6371
*****************************
FUNCTION Atan2

LPARAMETERS ;
	tnY ;
	, tnX

LOCAL ;
	lnRetVal

* Parameters must be passed in radians

* From Wikipedia definition of ATan2() function:
IF tnX > 0
	lnRetVal = ATAN( tnY / tnX )

ELSE
	IF tnX < 0
		IF tnY < 0
			lnRetVal = ATAN( tnY / tnX ) - PI( )
		
		ELSE
			lnRetVal = ATAN( tnY / tnX ) + PI( )
		
		ENDIF
	
	ELSE
		IF tnY > 0
			lnRetVal = PI( ) / 2
		
		ELSE
			IF tnY < 0
				lnRetVal = - ( PI( ) / 2 )
			
			ELSE
				* Return value where both X and Y are zero is undefined, we'll return 0:
				lnRetVal = 0
			
			ENDIF
		
		ENDIF
	
	ENDIF

ENDIF

RETURN lnRetVal
*****************************
I tested this exactly once, with a known distance near my home. I got exact lat/lons from Google Maps using a cool site http://www.getlatlon.com/ . I moved the crosshairs to two locations that I know are 350m apart on a curved track, and copy/pasted the exact lat/lons the site provided. The functions above returned 0.29km ( 290m ) distance between them, which is quite believable.

I have absolutely no idea what the lat/lons are in your code samples. If they're some sort of military grid system you'll need to convert them to conventional lat/lon.
Regards. Al

"Violence is the last refuge of the incompetent." -- Isaac Asimov
"Never let your sense of morals prevent you from doing what is right." -- Isaac Asimov

Neither a despot, nor a doormat, be

Every app wants to be a database app when it grows up
Previous
Next
Reply
Map
View

Click here to load this message in the networking platform