Geod. Az?

[+bicknell]

I want to make sure I'm interpreting this correctly, because I'm about to write a script to process some files for me....

 

I'll pick on FE2483 so we have a concrete example.

 

We have:

 

FE2483 *CURRENT SURVEY CONTROL

FE2483 ___________________________________________________________________

FE2483* NAD 83(1995)- 35 05 53.99253(N) 089 14 53.51921(W) ADJUSTED

FE2483* NAVD 88 - 163. (meters) 535. (feet) SCALED

FE2483 ___________________________________________________________________

and

 

FE2483|---------------------------------------------------------------------|

FE2483| PID Reference Object Distance Geod. Az |

FE2483| dddmmss.s |

FE2483| FE2484 GRAND JCT MUNICIPAL TANK APPROX. 7.8 KM 1331930.9 |

FE2483| CG3191 BEA RM 1 30.810 METERS 15704 |

FE2483| CG3190 BEA AZ MK 1635620.8 |

FE2483| CG3192 BEA RM 2 25.257 METERS 29953 |

FE2483|---------------------------------------------------------------------|

If we focus on reference mark #1 it's 30.81 Meters in the direction of 157 degrees 4 minutes from the mark. So, if I take the coordinates and I do some sort of calculation I can come up with the coordinates of the reference mark.

 

As you might guess, I'm interested in putting these in my GPS. After all, I'm lazy, and if I can walk over and look down for an easy recovery, that's good for me.

 

It would also allow me to project out a line for the azimuth mark. The lack of a distance is annoying, but if I project out a point 0.5 miles away or something, I think I can make my GPS draw a line, so at least I have a line to look along for the mark.

 

Now, it appears the NGS program "FORWARD" does this computation. If I put in the two marks above using the forward web page http://www.ngs.noaa.gov/cgi-bin/Inv_Fwd/forward2.prl, I get LAT = 35 5 54.99229 North LON = 89 14 53.51921 West.

 

Here's part that makes no sense to me. I would have expected this to be South East of the mark (157 degrees), however, plotting the output of forward on Google Earth places it North East of the mark. So I clearly don't understand something about Geodetic Azimuth.

 

So really, three questions. Do I have this right so far? Can someone explain why 157 degrees is North East? And most importantly, where do I find the forumla, or if it's really complicated source code to do the conversion? The NGS took is ok for one at a time, but I want to write up some perl scripts or something to generate points for all of the reference marks for all the trangulation stations in a county and output a kml file to load into google earth, and a gpx file to load into Mapsource and thus on to my GPS....

 

Lastly, a suggestion for the NGS. Can you just put the lat-long of the reference marks on the box score? That would be SO much easier.

[+Klemmer]

Google the word "Vincenty". You will find many sites that have the full sets of equations for geodetic calculation of point to point, or a point projection via azimuth & distance. Pretty serious math. One of them has a set of excel spreadsheets I have used that seem VERY accurate. Helped me find an RM.

 

Hope this helps. Good luck!

[+Kewaneh & Shark]

Grid azimuth is simply using true North as a basis of bearing. An azimuth bearing of 157 degrees should be, as you suspect, south and east of the original position. I think you may have put the data into the FORWARD incorrectly. When I look at your numbers, your results look to me to be due North of your starting position.

 

From the datasheet as shown on your post:

Lat: 35 05 53.99253 (N) Long: 089 14 53.51921 (W)

 

Your results as shown on your post:

Lat: 35 05 54.99229 (N) Long: 089 14 53.51921 (W)

 

Deltas:

Lat: 00 00 00.99976 (N) Long: 000 00 00.00000 (W)

 

Obviously not the intended position from the setting survey crew, or the intended results from you.

 

When I run your data through FORWARD, these are my results:

 

Ellipsoid : GRS80 / WGS84 (NAD83)

Equatorial axis, a = 6378137.0000

Polar axis, b = 6356752.3141

Inverse flattening, 1/f = 298.25722210088

 

First Station : FE2483

----------------

LAT = 35 5 53.99253 North

LON = 89 14 53.51921 West

 

Second Station : CG3191

----------------

LAT = 35 5 53.07179 North

LON = 89 14 53.04520 West

 

Forward azimuth FAZ = 157 4 0.0000 From North

Back azimuth BAZ = 337 4 0.2725 From North

Ellipsoidal distance S = 30.8100 m

 

Deltas:

Lat: 00 00 00.92074 (S) Long: 000 00 00.47401 (W)

 

You may need to be sure that you completely 'spell out' all of the bearings when you use the FORWARD program. For example, the azimuth as shown on the datasheet is 157 04. When that is put into the FORWARD program, it should be input as 157 04 00.

 

Hope this helps.

 

- Kewaneh

Edited May 5, 2006 by Kewaneh & Shark

[holograph]

And most importantly, where do I find the forumla, or if it's really complicated source code to do the conversion? The NGS took is ok for one at a time, but I want to write up some perl scripts or something to generate points for all of the reference marks for all the trangulation stations in a county and output a kml file to load into google earth, and a gpx file to load into Mapsource and thus on to my GPS....

 

Over the short distances from stations to reference marks, you don't need Vincenty. A simple "flat earth" calculation based on trigonometry will be just as good.

 

1. Calculate the DX and DY distance to the reference mark using DX = d * cos(90° - az), DY = d * sin(90° - az), where d = the distance in the box score, and az is the azimuth.

 

2. Convert the DY value to delta degrees latitude based on deltaLAT = (180 * DY) / (r * PI), where r is the radius of the earth (6389000 meters is a good enough value for this calculation).

 

3. Convert the DX value to delta degrees longitude based on deltaLON = (180 * DX) / (r * cos(station lat)* PI).

 

4. Compute the latitude of the reference mark = LAT + deltaLAT, where LAT is the latitude of the station.

 

5. Compute the longitude of the reference mark = LON + deltaLON, where LON is the longitude of the station.

 

Your example:

 

1. DX = 30.81 * cos(90 - 157°04') = 12.01 meters.

DY = 30.81 * sin(90 - 157°04') = -28.37 meters.

 

2. deltaLAT = (180 * -28.37) / (6389000 * 3.142) = -0.0002544 degrees = -0.916 seconds.

 

3. deltaLON = (180 * 12.01) / (6389000 * cos(35°05'54") * 3.142) = 0.0001316 degrees = 0.474 seconds.

 

4. Lat RM1 = 35°05'53.993" - 0.916" = 35°05'53.077"

 

5. Lon RM1 = -89°14'53.518" + 0.474 = -89°14'53.044"

 

This position is within 1 mm of the position calculated by the much more complicated Vincenty algorithm. Compared to the 4.5 mm uncertainty due to the azimuth being given to only the nearest minute, and the several meter uncertainty in your GPS receiver, this is more than enough accuracy for your purposes.

[+bicknell]

Holograph to the rescue! That I can deal with!

 

See the script below. I'm still testing, it may be broken, but initial testing looks good. Run "ref2gpx.pl ngs-file.dat > reference-gpx-file.gpx". The result is a GPX file with "PID-x" for the various PIDs. Note! The x in -1 does NOT necessarily match the RM #. I'll work on that.

 

If anyone else can use this feel free, if some others can check the output for accuracy that would be good too.

 

$ cat ref2gpx.pl

#!/usr/bin/perl

 

use constant PI => 4 * atan2 1, 1;

 

# This program assumes north america.

sub deg2rad { PI * $_[0] / 180 }

 

print <<EOM;

<?xml version="1.0" encoding="ISO-8859-1" standalone="yes"?>

<gpx version="1.0" creator="ref2gpx.pl"

xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance"

xmlns="http://www.topografix.com/GPX/1/0"'>http://www.topografix.com/GPX/1/0"

xsi:schemaLocation="http://www.topografix.com/GPX/1/0

http://www.topografix.com/GPX/1/0/gpx.xsd">

EOM

 

while (<>) {

if (m/^ (\S+)\s+DESIGNATION -\s+(\S+)/) {

$des = $1;

}

if (m/^ (\S+)\* NAD 83\(1995\)- (\d+) (\d+) (\d+\.\d+)\(N\)\s+(\d+) (\d+) (\d+\.\d+)\(W\)\s+ADJUSTED/) {

$count = 1;

$mark = $1;

$lat = $2 + ($3 / 60) + ($4 / 3600);

$lon = $5 + ($6 / 60) + ($7 / 3600);

$lon *= -1;

# print "Found adjusted mark $mark $lat $lon\n";

}

if (m/^ (\S+)\|\s+(\S+) (.*RM\s+\d+)\s+(\d+\.\d+)\s+METERS\s+(\d+)\s+\|/) {

$name = $1;

$reference = $3;

$distance = $4;

$azimuth = $5;

 

if (length($azimuth) == 9) {

($d, $m, $s) = unpack("A3A2A4", $azimuth);

} elsif (length($azimuth) == 5) {

($d, $m) = unpack("A3A2", $azimuth);

$s = 0;

} elsif (length($azimuth) == 3) {

$d = $azimuth;

$m = $s = 0;

} else {

print "Can't parse azimuth $azimuth.\n";

}

 

if ($name eq $mark) {

# print "Compute distance from $name $distance meters bearing $azimuth to $reference.\n";

$rotation = $d + ($m / 60) + ($s / 3600);

$dx = $distance * cos(deg2rad(90 - $rotation));

$dy = $distance * sin(deg2rad(90 - $rotation));

# print "dx = $dx, dy = $dy\n";

$dlat = (180 * $dy) / (6389000 * 3.142);

$dlon = (180 * $dx) / (6389000 * 3.142 * cos(deg2rad($lat)));

$newlat = $lat + $dlat;

$newlon = $lon + $dlon;

# print "New lat $newlat ($dlat) new lon $newlon ($dlon)\n";

printf("<wpt lat=\"%.8f\" lon=\"%.8f\">\n", $newlat, $newlon);

print <<EOM;

<name>$name-$count</name>

<desc>$reference</desc>

<sym>Geocache</sym>

</wpt>

EOM

$count++;

}

}

}

 

print "</gpx>\n";

 

Edited May 6, 2006 by bicknell

[+PFF]

Hi, Bicknell:

 

Yes, the FORWARD program will calculate the precise coordinates for you. In this case, the bearing must be entered as 157 04 00, since that is the degree of accuracy represented in the box. I have pasted the results for Ref Mark 1, below. It does, indeed, plot to the southeast, just as you thought it should.

 

For the Az Mark, you can use the same proceedure, but with various distances. This will produce a line along which you can walk, using GPS waypoints as a guide. I typically calcuate coordates at the following distances from the station: 0.4 mile, 0.45 mile, 0.5 mile, and 0.55 mile. Once in the field, I start at the coordinates closest to the station, and then walk along the bearing--looking to the right and to the left, so I don't miss anything.

 

Best regards,

Paul

 

Output from FORWARD

 

 

First Station : Main Station

----------------

LAT = 35 4 53.99253 North

LON = 89 14 53.31921 West

 

Second Station : REF MARK 1

----------------

LAT = 35 4 53.07179 North

LON = 89 14 52.84530 West

 

Forward azimuth FAZ = 157 4 0.0000 From North

Back azimuth BAZ = 337 4 0.2724 From North

Ellipsoidal distance S = 30.8100 m

[+GEO*Trailblazer 1]

Well I am not that good at programming.

 

But there is another way I do it.

 

I use the GPS Trackmaker..now functions and is interactive with GOOGLE EARTH.

 

You can set all your coordinates up in the datum(s) DGTV TERRAIN VIEWER TRACK MAKER you prefer or the ones listed for the benchmarks.

 

Any how you can add the lines of azimuths and distaces then create a waypoint at that point and you have the reading (coordinates) of the reference marks and azimuth.

 

You can also change the icons to Benchmark icons and the other to the many other things related.

 

Here is the close up of GF0914 EAGLE TRIANGULATION STATION and the 2 reference marks.

All the tools are available within the freeware program.

 

This is the 3 marks and the azimuth.

 

This it the tying the triangulation stations other stations and USGS TOPO BM's and reference marks close to home.

 

It is the easy way for me.