I got some comments on my comments about calculating azimuth yesterday that are worth sharing. Bill Huber provided me with some code to do the calculation with fewer nested If...then statements (and if you can be sure that your line end-point's x-coordinates are not the same, you can even eliminate the remaining If statement.) The code, translated into MapBasic appears below. Note also that I've added a floating point version of the Mod function, since MapInfo's Mod statement supports only integers.
Mark Neuhaus also pointed out that projected coordinate systems do not align to north anywhere except along their central meridian, and if you want your azimuths measured relative to north you will find them to be less accurate the further from the center of the projection they are. The best answer to that is to use projections that are as local to your lines as possible, and remember that accuracy degrades with distance. Again, using Pythagorus' method for making measurements on a round surface is at best a compromise. Really, if you want to do it right, use spherical trig. Working examples of great circle calculations in MapBasic are available from Jacques Paris' MB Resources site at http://www.paris-pc-gis.com/mb_r/mbr_start.htm (yours truly wrote the examples.) The following is a MapBasic program to calculate azimuths assuming that the world is flat and that the coordinates are cartesian. The best ideas in here come from Bill Huber. 'LineAzimuth.mb '============== 'Demonstration program illustrating how to calculate azimuth from 'line coordinates. Define PI2 1.570796327 '(PI / 2) Define ERR_DIVIDE_BY_ZERO 309 Declare Sub Main Declare Function Azimuth (ByVal dx As Float, ByVal dy As Float) As Float Declare Function fMod (ByVal fNum As Float, ByVal fModifier As Float) As Float Sub Main Dim x0, y0, x1, y1 As Float 'Clear message window Print Chr$(12) 'Set coordinates of starting point x0 = 0 y0 = 0 'Set coordinates of ending point x1 = 0.9 y1 = -0.1 'Print the azimuth angle that ending point bears to 'starting point measured clockwise from north. Print Azimuth (x1-x0, y1-y0) End Sub Function Azimuth ( ByVal dx As Float, 'difference in x coordinates of line ByVal dy As Float) 'difference in y coordinates of line As Float 'Azimuth in degrees (0=North, 90=East, etc.) 'Calculates the azimuth between two coordinates given their 'Differences in x and y. Dim a As Float If dx = 0 Then a = Sgn(dy) * PI2 Else a = Atn(dy/dx) + (1-Sgn(dx)) * PI2 End If Azimuth = fmod(450 - a * 90 / PI2, 360) End Function Function fMod ( ByVal fNum As Float, 'the r in r mod s ByVal fModifier As Float) 'the s in r mod s As Float 'Remainder of r/s 'This is a floating point version of the Mod statement, implemented 'as a function. (MapInfo's Mod statement is integer-based). If fModifier = 0 Then 'Error. Raise divide by zero error. Error ERR_DIVIDE_BY_ZERO Else fMod = fNum - Int (fNum / fModifier) * fModifier End If End Function -- - Bill Thoen ------------------------------------------------------------ GISnet, 1401 Walnut St., Suite C, Boulder, CO 80302 tel: 303-786-9961, fax: 303-443-4856 mailto:[EMAIL PROTECTED], http://www.gisnet.com/ ------------------------------------------------------------ _______________________________________________________________________ List hosting provided by Directions Magazine | www.directionsmag.com | To unsubscribe, send e-mail to [EMAIL PROTECTED] and put "unsubscribe MapInfo-L" in the message body.
