GPX and projection ok

README.md

@@ -79,6 +79,26 @@ $ screen -X -S sesddboat quit
 where sesddboat is the session name the has been used at the start.
 
 
+## working with GPS
+
+The following packets must have been installed :
+```
+$ sudo apt install python3-gpxpy python3-pyproj
+```
+
+The GPS is acquired and the GPGLL message is recorded in GPX format (so you can plot it in GeoPortail for example) in the **tst.gpx** file with the command:
+```
+$ python3 tst_gpx.py
+```
+
+In the NMEA message GPGLL, the latitude and longitude are given in DDmm.mmmm with DD in degrees and mm.mmmm in decimal minutes. To get the value in decimal degrees we do DD + mm.mmmm/60.0
+
+Using a projection from degrees (longitude,latitude in WGS84) to meters (UTM zone 30N) we can compute the distance in meter to a reference point. The following command executes a test:
+```
+$ python3 tst_proj.py
+```
+
+
 
 
 

tst_gpx.py

@@ -3,6 +3,7 @@ import sys
 import gps_driver_v2 as gpsdrv
 import gpxpy.gpx
 
+# record GPS (GPGLL ) in a GPX file
 # requires install of gpxpy
 # sudo apt install python3-gpxpy
 
@@ -32,7 +33,7 @@ gpx.tracks.append(gpx_track)
 gpx_segment = gpxpy.gpx.GPXTrackSegment()
 gpx_track.segments.append(gpx_segment)
 
-tmax = 10*60.0 # 10 minutes max
+tmax = 10*60.0 # 10 minutes max (can be stopped before the end with Ctrl+C)
 t0 = time.time()
 while True:
     print ("---------------------------------------------------")

tst_proj.py

@@ -7,55 +7,41 @@ import gpxpy
 from pyproj import Proj, transform
 import numpy as np
 
-#projMeter = Proj("+init=EPSG:2154",preserve_units=False) # Lambert-93 - France 
-#projMeter = Proj("+init=EPSG:5627",preserve_units=False) # ED 50
-projMeter = Proj("+init=EPSG:32620",preserve_units=False) # UTM Zone 30
-projDegre = Proj("+init=EPSG:4326") # WGS84 in degrees and not EPSG:3857 in meters)
-
-# reference point "the flag"
-reference_lat = 48.4188
-reference_lon = -4.4725
-reference_x,reference_y = transform(projDegre,projMeter,reference_lon,reference_lat)
-print (reference_x,reference_y,reference_lat,reference_lon)
-x = reference_x
-y = reference_y
-# reference point center of the rugby field
-reference_lat = 48.4185
-reference_lon = -4.4739
-reference_x,reference_y = transform(projDegre,projMeter,reference_lon,reference_lat)
-print (reference_x,reference_y,reference_lat,reference_lon)
-dx = x - reference_x
-dy = y - reference_y
-distance = np.sqrt(dx*dx+dy*dy)
-print (dx,dy,distance)
-exit()
-
+# define the projection from WGS84 (lon,lat in degrees) to UTM (meters)
+# UTM zone for Brittany is 30N 
+projDegree2Meter = Proj("+proj=utm +zone=30, +north +ellps=WGS84 +datum=WGS84 +units=m +no_defs")
+
+iref = 0 # select reference point
+if iref==0:
+    # reference point "the flag"
+    reference_lat = 48.4188
+    reference_lon = -4.4725
+elif iref==1:
+    # reference point center of the rugby field
+    reference_lat = 48.4185
+    reference_lon = -4.4739
+
+reference_x,reference_y = projDegree2Meter(reference_lon,reference_lat)
+print ("ref",reference_x,reference_y,reference_lon,reference_lat)
+
+# open the ddboat gps log file (gpx format) 
 gpx_file = open('tst_ensta_stadium_20210929.gpx', 'r')
 
-
-init = True
 gpx = gpxpy.parse(gpx_file)
+i = 0
 for track in gpx.tracks:
     for segment in track.segments:
         for point in segment.points:
             #print('Point at ({0},{1}) -> {2}'.format(point.latitude, point.longitude, point.elevation))
-            x,y = transform(projDegre,projMeter,point.longitude, point.latitude)
-            lat,lon = transform(projMeter,projDegre,x,y)
-            print (x,y,lat,lon,point.latitude,point.longitude)
+            x,y = projDegree2Meter(point.longitude, point.latitude)
+            lat,lon = projDegree2Meter(x,y,inverse=True)
+            #print (x,y,lat,lon,point.latitude,point.longitude)
             dx = x - reference_x
             dy = y - reference_y
             distance = np.sqrt(dx*dx+dy*dy)
-            print (dx,dy,distance)
-
-
-
-# test project lat,lon (deg) to meters
-#from pyproj import Proj, transform
-#inProj  = Proj("+init=EPSG:2263",preserve_units=True)
-#outProj = Proj("+init=EPSG:4326") # WGS84 in degrees and not EPSG:3857 in meters)
-#x2,y2 = -74.2700000001129, 40.46999999990434
-#print (x2,y2)
-#x1,y1 = transform(outProj,inProj,x2,y2)
-#print (x1,y1)
-#print(transform(inProj,outProj,x1,y1))
-
+            # print  first distance
+            if i == 0:
+                print (dx,dy,distance)
+            i += 1
+# print last distance
+print (dx,dy,distance)