Projections - Convert your Survey Data to a Georeferenced GIS file
Posted by Andre Kruger in Articles
My first blog post will tackle the basis of working with GIS data. It is also a topic that is normally not covered in tertiary courses training for civil engineers. This topic is projections. All I know about projections was self taught. It could be a very complicated subject if you want it to be. But to just understand enough of the basics should be sufficient for engineers.
To learn about projections we are going to take a point text file from a surveyor and convert it into a georeferenced GIS file. We will create a shapefile, a GeoPackage and a KML which can be loaded into Google Earth.
Normally I would not use this methodology to load the survey point file into GIS. My normal workflow would be as follows:
- Load the survey into a spreadsheet. (I prefer LibreOffice Calc but MS Excel does the same)
- Covert the survey points to Transverse Mercator as discussed in Point 4 below.
- Aplly column headings with an "X" and "Y" for the coordinates.
- Save the file as a CSV
- Exit the spreadsheet software
- Load the CSV file into QGIS
- It automatically recognises the "X" and "Y" as coordinates
- Load and view the data.
- Under the "Properties" of the layer set the coordinate reference system.
- For a survey in South Africa just search for "HBK" and choose the applicable one as discussed in point 5 below.
This excercise below automates the above workflow and show cases what is possible with Python and Jupyter.
1. Input and Output¶
1.1 Input¶
- The file name of the surveyor data. Type is a string.
- The central meridian. Discussed inpoint 5. Type is integer.
1.2 Output¶
All files are georeferenced.
- An ESRI Shapefile
- A Geopacakage (GPKG) file.
- A KML file. To view your survey's location in Google Earth.
- An embedded OpenStreetMap map that zooms to the area of your survey.
# Input
survey_file = 'BANDELIER8.TXT'
central_meridian = 29
2. Inspect the Survey File¶
The easiest way to call shell from Jupyter is to prepend an exclamation point ! to a shell command.
To list the files in your current working directory I would normall use !ls which is the list command on Unix like systems. The equivalent on Windows will be !dir. But to keep it cross platform and to teach you Python at the same time we will use Python to list the directory contents
import os
os.listdir('./')
The Unix command for looking at the first rows of a text file is !head. But sticking to Python.
with open(survey_file) as source:
for i, line in enumerate(source):
if i > 7:
break
print(line, end='')
From the inspection of the file we can see that each record has got 4 fields. Something like an "ID" an X and a Y or an Y and X, we will determine that now. The last two fields is an Elevation and a Description. No Header.
3. The Transverse Mercator Projection¶
When I work with South African data I almost always use the Transverse Mercator projection. The reason for using this projection is as follows:
- It is easy to understand.
- All geographical and CAD packages supports it.
- It is easy to quickly convert survey data to Transverse Mercator as I will show now.
- If you receive projected aeriel photography it will almost certainly be in Transverse Mercator.
- The South African LiDAR data I have seen also uses Transverse Mercator.
How I understand it:
In South Africa your Y coordinate will always be negative and its absolute value will always be greater than two million. It indicates the amount of metres South of the equator. The X coordinate will be postive or negative depending on whether the coordinate is (+) East or (-) West of the central meridian.
The survey coordinates traditionally used by South African surveyors is known as South Oriented Transverse Mercator or TMSO for short. It is a crazy projection. Here is a excellent write-up or rant. It is the equivalent of the Transverse Mercator except that the X and Y axis are swapped and the signs on both are changed.
4. Loading the Survey Data and Changing it from TMSO to Transverse Mercator¶
So changing the coordinates is only a matter of taking the large coordinate of 2 million + making it the Y coordinate and making it negative. The other value will be your X coordinate. Remember to change the sign of this value as well. From negative to positive or vice versa.
import csv
header = False
with open(survey_file) as source:
survey = []
reader = csv.reader(source, delimiter='\t')
for i, line in enumerate(reader):
if header and i==0:
continue
if i < 5:
print(line)
x = float(line[1])*-1
y = float(line[2])*-1
elevation = float(line[3])
if line[4] == '':
description = line[0]
else:
description = line[4]
survey.append([x, y, elevation, description])
print('\nThere is', len(survey), 'points in the survey.\n')
unique = set()
for point in survey:
unique.add(point[3])
print('There is', len(unique), 'unique descriptions in the survey.')
print('The descriptions are:')
unique = list(unique)
unique.sort()
for i, item in enumerate(unique):
if i%5 == 0:
print()
print(item, end='\t\t')
print()
5. Determine the Central Meridian¶
Usually refered to in South Africa as Lo27 or Lo29 for example. Lo is short for Longitude of origin. I prefer to use the term Central meridian. Depending on where you are in the country it can be anything from 17 to 31 degrees in uneven increments.
If you know where the site is get a Easting value of the site on Google Earth or use a GPS coordinate reading from the site. This site above is near Bandelierkop in the Limpopo province. I read a coordinate of about 29°50' East or this site. This gives us a central meridian 29.
print('The chosen Central Meridian is', central_meridian, 'degrees.')
6. Write the Survey points to various file formats¶
File formats supported by this Fiona installation.
- r = read
- a = append
- w = write
import fiona
print(fiona.supported_drivers)
Write the survey data to a GeoPackage:
%%time
from collections import OrderedDict
crs = {'ellps': 'WGS84',
'k': 1,
'lat_0': 0,
'lon_0': central_meridian,
'no_defs': True,
'proj': 'tmerc',
'towgs84': '0,0,0,0,0,0,0',
'units': 'm',
'x_0': 0,
'y_0': 0}
schema = {'geometry': 'Point',
'properties': OrderedDict([('X', 'float'),
('Y', 'float'),
('Elevation', 'float'),
('Description', 'str')])}
with fiona.open('survey.gpkg', 'w', driver='GPKG', crs=crs, schema=schema) as sink:
for point in survey:
rec = {}
rec['geometry'] = {}
rec['geometry']['type'] = 'Point'
rec['geometry']['coordinates'] = (point[0], point[1])
rec['properties'] = {}
rec['properties']['X'] = point[0]
rec['properties']['Y'] = point[1]
rec['properties']['Elevation'] = point[2]
rec['properties']['Description'] = point[3]
sink.write(rec)
Convert the GeoPackage to an ESRI Shapefile and a KML file for Google Earth.
Bind the ogr2ogr command line tool to the binary in your Ananconda installation. It may be installed at other places as well.
%alias ogr2ogr /Users/andre/anaconda3/envs/eng/bin/ogr2ogr
%%time
%ogr2ogr -f "ESRI Shapefile" -overwrite "shp_survey.shp" "survey.gpkg"
%ogr2ogr -f "KML" -overwrite "kml_survey.kml" "survey.gpkg"
if 'map_survey.json' in os.listdir('./'):
os.remove('map_survey.json')
%ogr2ogr -f "GeoJSON" -t_srs EPSG:4326 "map_survey.json" "survey.gpkg"
with fiona.open('map_survey.json') as source:
print(source.bounds)
xmin, ymin, xmax, ymax = source.bounds
import folium
map_osm = folium.Map(location=[(ymin+ymax)/2, (xmin+xmax)/2])
map_osm.fit_bounds([[ymin, xmin], [ymax, xmax]])
#map_osm.geo_json(geo_path='map_survey.json')
map_osm