Calculate bearing of lines in QGIS

 Bearing refers to the direction or orientation of a line or object relative to a fixed reference point, typically north. It is often measured in degrees, with 0° representing north, 90° representing east, 180° representing south, and 270° representing west. In simpler terms, bearing is like a compass pointing in a specific direction.

Bearings are commonly used in the field of Navigation, Surveying, Astronomy and Engineering. In some instances, different discipline use different terms such as Azimuth, Heading, Direction, Course etc to refer to bearing.

As at the time of writing, calculating the bearing of a line in QGIS doesn't work well for the context of surveying in Nigeria. Nigerian land surveyors work with Whole Circle Bearing (WCB) which is a method of expressing direction using degrees measured clockwise from north. WCB values range from 0° to 360° where North is 0°, East is 90°, South is 180° and West is 270°.

In QGIS under the "Measure" tool/icon, you will find "Measure Bearing" tool. Currently, the bearing is in Quadrant format.

In this post, lets explore how to calculate the Whole Circle Bearing of lines and use it to label the lines.

Electronic Nautical Charts in QGIS

 Different organizations are responsible for the production of Electronic Navigational Charts a.k.a Electronic Nautical Charts in different countries. In Nigeria, the Nigerian Navy Hydrographic Office (NNHO) commenced the production of Nautical Charts in 2018.

In the US, charts are produced by National Oceanic and Atmospheric Administration (NOAA). In the UK, it is UK Hydrographic Office (UKHO), they called their charts 'Admiralty Charts'. Other countries have similar agencies and organizations that makes these charts.

Most of these organizations have provided access to the charts in WMS/WMTS layers in QGIS software. In this post, I will be using the service provided by ArgoGIS, a marine geographic information based on nautical charts as web map service.

The ArgoGIS WMS url is: https://enc.argogis.com/wms/wms_enc

To add it to QGIS, navigate to the Layer menu and add WMS/WMTS layer, then provide the url above. The result around 'Bonny Town' looks like the image below.

These charts are used by oceanographers, mariners, kayakers, navigators, boat captains, commercial fishermen, nautical tourists, coastal engineers, marine scientists, hydrographic surveyors (hydrographers) etc to plan voyages and navigate ships safely by providing detailed information about water depths, locations of navigational hazards like rocks and shoals, positions of aids to navigation, anchorages, coastlines, and other features essential for safe navigation on waterways (in streams, rivers, lakes, oceans, seas, and coastal areas), particularly in coastal areas and open seas.

Happy navigaton!

AutoLISP Program that iterates through selected lines, calculates their angles and lengths, and then annotates them

 In other words, the autoLISP program will select a line, get its bearing and distance and write the values on beside the line. The distance is in meters and the bearing format is in Whole Circle Bearing and DD°MM'SS

(defun c:BRGDIST (/ ss i entData ptStart ptEnd dx dy angleRad angleDeg wcb bearingStr dist midpoint)
  ; Function to convert decimal degrees to DMS (degrees, minutes, seconds)
  (defun DEGtoDMS (deg / d m s)
    (setq d (fix deg)) ; Degrees
    (setq m (fix (* (- deg d) 60))) ; Minutes
    (setq s (* (- (* (- deg d) 60) m) 60)) ; Seconds
    (strcat (itoa d) "°" (itoa m) "'" (rtos s 2 2) "\"") ; Format as DD°MM'SS"
  )

  ; Select lines
  (setq ss (ssget '((0 . "LINE"))))
  (if ss
    (progn
      (setq i 0)
      (repeat (sslength ss)
        ; Extract line data
        (setq entData (entget (ssname ss i)))
        (setq ptStart (cdr (assoc 10 entData))) ; Start point
        (setq ptEnd (cdr (assoc 11 entData)))   ; End point
        
        ; Calculate bearing
        (setq angleRad (angle ptStart ptEnd)) ; Angle in radians (from X-axis, counterclockwise)
        (setq angleDeg (/ (* angleRad 180.0) pi)) ; Convert to degrees
        
        ; Convert to Whole Circle Bearing (clockwise from north, 0°-360°)
        (setq wcb (rem (+ (- 90 angleDeg) 360) 360)) ; Adjust for WCB
        (setq bearingStr (DEGtoDMS wcb)) ; Convert WCB to DMS format
        
        ; Calculate distance
        (setq dist (rtos (distance ptStart ptEnd) 2 2)) ; Distance with 2 decimal places
        
        ; Calculate midpoint
        (setq midpoint (list
          (/ (+ (car ptStart) (car ptEnd)) 2.0)
          (/ (+ (cadr ptStart) (cadr ptEnd)) 2.0)
        ))
        
        ; Create text at midpoint
        (entmake
          (list
            '(0 . "TEXT")
            (cons 10 midpoint)          ; Insertion point
            (cons 40 5)               ; Text height
            (cons 1 (strcat bearingStr "      " dist"m")) ; Bearing (DMS) + distance
            (cons 50 angleRad)          ; Rotation angle (aligned with line)
            (cons 7 "Standard")         ; Text style
          )
        )
        (setq i (1+ i))
      )
    )
  )
  (princ)
)

This will create a new layer for the annotation

(defun c:BRGDIST (/ ss i entData ptStart ptEnd dx dy angleRad angleDeg wcb bearingStr dist midpoint)
  ; Function to convert decimal degrees to DMS (degrees, minutes, seconds)
  (defun DEGtoDMS (deg / d m s)
    (setq d (fix deg)) ; Degrees
    (setq m (fix (* (- deg d) 60))) ; Minutes
    (strcat (itoa d) "°" (itoa m) "'") ; Format as DD° MM'
    ;; (setq s (* (- (* (- deg d) 60) m) 60)) ; Seconds
    ;; (strcat (itoa d) "°" (itoa m) "'" (rtos s 2 2) "\"") ; Format as DD° MM' SS"
  )

  ; Create or use the annotation layer
  (if (not (tblsearch "LAYER" "BRG_DIST"))
    (command "._LAYER" "_M" "BRG_DIST" "_C" "1" "" "") ; Create layer and set color to red
  )
  (setvar "CLAYER" "BRG_DIST") ; Set current layer to BRG_DIST

  ; Select lines
  (setq ss (ssget '((0 . "LINE"))))
  (if ss
    (progn
      (setq i 0)
      (repeat (sslength ss)
        ; Extract line data
        (setq entData (entget (ssname ss i)))
        (setq ptStart (cdr (assoc 10 entData))) ; Start point
        (setq ptEnd (cdr (assoc 11 entData)))   ; End point
        
        ; Calculate bearing
        (setq angleRad (angle ptStart ptEnd)) ; Angle in radians (from X-axis, counterclockwise)
        (setq angleDeg (/ (* angleRad 180.0) pi)) ; Convert to degrees
        
        ; Convert to Whole Circle Bearing (clockwise from north, 0°-360°)
        (setq wcb (rem (+ (- 90 angleDeg) 360) 360)) ; Adjust for WCB
        (setq bearingStr (DEGtoDMS wcb)) ; Convert WCB to DMS format
        
        ; Calculate distance
        (setq dist (rtos (distance ptStart ptEnd) 2 2)) ; Distance with 2 decimal places
        
        ; Calculate midpoint
        (setq midpoint (list
          (/ (+ (car ptStart) (car ptEnd)) 2.0)
          (/ (+ (cadr ptStart) (cadr ptEnd)) 2.0)
        ))
        
        ; Create text at midpoint
        (entmake
          (list
            '(0 . "TEXT")
            (cons 10 midpoint)          ; Insertion point
            (cons 40 5)               ; Text height
            (cons 1 (strcat bearingStr "   " dist "m")) ; Bearing (DMS) + distance
            (cons 50 angleRad)          ; Rotation angle (aligned with line)
            (cons 7 "Standard")         ; Text style
            (cons 8 "BRG_DIST")         ; Layer
          )
        )
        (setq i (1+ i))
      )
    )
  )
  (princ)
)

Happy Coding!

Generating AutoCAD Script (.SCR) file from GPS/GNSS Data

 Recently, I was working on a project where I had to label some points in several AutoCAD files from data points obatined from GPS/GNSS reciever. So, I wrote this Python code to 'Generate AutoCAD Script (.SCR) file from GNSS Data'.

easting = '356662'
northing = '946418'

# construct text string eg: _TEXT easting,northing 10 90 PA10
script_str_txt = f'_TEXT {easting},{northing} 10 90 {name}\n'

# construct point string eg: _POINT easting,northing
script_str_pt = f'_POINT {easting},{northing}\n'

with open(f"{file_num}_{sch_name}.scr", "a") as f:
    f.write(script_str_txt)

The result of the script:

A more practical application code is this one below. It reads a large CSV file containing school names, then groups the data by the school names and generate .scr file for each group.

sch_dict = {'GSS TOTO WEST TOTO LGC':1, 'GJSS ASOPADA KARU LGC':2, 'GJSS KURMIN TAGWAYE AKWANG':3, 'GJSS LAFIA SOUTH LAFIA LGC':4, 'GJSS SAMBGOBARAU KOKONA LGC':5, 'GOV DAY SEC OBI LGC':6, 'GS LOKO NASARAWA LGC':7, 'GSS ALUSHI KEANA LGC':8, 'GSS AZUB CENTRAL LAFIA LGC':9, 'GSS EFUGBORINGO DOMA LGC':10, 'GSS GUDI AKWANGA LGC':11, 'GSS KEKURA AWE LGC':12, 'GSS MAMA WAMBA LGC':13, 'GSS SABONGARI KEFFI KEFFI LGC':14, 'ISLAMIYYA LGEA KEANA LGC':15, 'LGEA DHIZILA OBI LGC':16, 'LGEA EFURIGBORINGO DOMA LGC':17, 'LGEA EZZEN MADA STATION EGGON':18, 'LGEA GUDIGE TSOHO NASARAWA LGC':19, 'LGEA IBI AWE LGC':20, 'LGEA KIGBUNA LAFIA LGC':21, 'LGEA MADA HILL AKWANGA LGC':22, 'LGEA PILOT CEN UMAISHA TOTO LG':23, 'LGEA PRI ADOKASA KARU LGC':24, 'LGEA PRI KONZA WAMB LGC':25, 'LGEA PRI NASARAWA EAST NAS LGC':26, 'LGEA REFUGE CAMP LAFIA LGC':27, 'LGEA SAKWATO KOKONA LGC':28, 'LGEA TRANSFER UMAISHA TOT LGC':29, 'NADP PROJECT LAFIA LGC':30, 'RCM GSS NAS EGGON LGC':31}


# Read the GNSS created spread sheet....
df = pd.read_csv(r"sch-24.csv")

group_df = df.groupby('School Name')
group_keys = list(group_df.groups.keys())

for s in group_keys:
    temp_df = group_df.get_group(s)

    for idx, row in temp_df.iterrows():
        easting = row['X']
        northing = row['Y']
        name = row['Name']
        sch_name = row['School Name']
        file_num = sch_dict[sch_name]
    
        # construct text string eg: _TEXT easting,northing 10 90 PA10
        script_str_txt = f'_TEXT {easting},{northing} 10 90 {name}\n'
        # construct point string eg: _POINT easting,northing
        script_str_pt = f'_POINT {easting},{northing}\n'

        with open(f"{file_num}_{sch_name}.scr", "a") as f:
            f.write(script_str_txt)

        with open(f"{file_num}_{sch_name}.scr", "a") as f:
            f.write(script_str_pt)
        
    # break

print('Done...')

Happy coding!