pyproj¶
PROJ is a generic coordinate transformation software, that transforms coordinates from one coordinate reference system (CRS) to another. This includes cartographic projections as well as geodetic transformations.
PROJ.4 was started in the 1980s, and Version 4 has been around since 1995. Now there is a version 5 it is referred to as PROJ (the history of the project can be found on Wikipedia). In version 4 everything was converted by passing through WGS84. If you are interested in projections, have a read of https://proj4.org/development/migration.html for more details. PROJ4 could only guarantee metre accuracy for many transformations, so PROJ5 provides much higher accuracy. A fund raising effort was started to improve the way projections are handled at https://gdalbarn.com/ - as projections are critical to the entire open source stack.
It is used throughout GIS software, for example in ArcGIS, QGIS, and MapServer.
pyproj
is a Cython wrapper to provide python interfaces to PROJ. Access to PROJ is also available through many other languages (e.g. C#, R).
Links:
- proj sourcecode at Github: https://github.com/OSGeo/proj.4/
- proj documentation at https://proj4.org/
- pyproj documentation at https://jswhit.github.io/pyproj/
- pyproj sourcode at https://github.com/jswhit/pyproj
Projection Definitions¶
To install pyproj you can run the following command from the command line (Windows or Linux). It is however preinstalled on OGGeoLive.
pip install pyproj
http://epsg.io/ is a very useful site to find project definitions. For example to see projections relating to Ireland you can run the following search: http://epsg.io/?q=ireland The Irish Transverse Mercator projection is found at http://epsg.io/2157 - note the structure URLs which make these easy to access programatically.
Each projection on the site hasPROJ.4 definitions
# Irish Transverse Mercator http://epsg.io/2157.proj4
+proj=tmerc +lat_0=53.5 +lon_0=-8 +k=0.99982 +x_0=600000 +y_0=750000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs
# Irish National Grid https://epsg.io/29902.proj4
+proj=tmerc +lat_0=53.5 +lon_0=-8 +k=1.000035 +x_0=200000 +y_0=250000 +ellps=mod_airy +towgs84=482.5,-130.6,564.6,-1.042,-0.214,-0.631,8.15 +units=m +no_defs
Session Code¶
Note - when saving the script, don’t name it pyproj.py
as this will clash with the name of the pyproj
library that needs
to be imported and cause the script to fail!
"""
This script does the following:
+ Downloads a list of Irish towns and cities from https://data.gov.ie in a CSV format
+ Takes the ITM coordinates for each settlement and converts them to Lat/Lng
The dataset details can be found at https://data.gov.ie/dataset/settlements-ungeneralised-osi-national-statistical-boundaries/resource/c7d0e002-8fcf-4e31-820e-5d628f41eea2
Note there is no metadata indicating the coordinate projection!
As the Irish settlement names include special characters such as a fada the text needs to be encoded in UTF (Universal Text Format).
Otherwise errors such as the following will occur:
UnicodeEncodeError: 'ascii' codec can't encode character u'\ufeff' in position 0: ordinal not in range(128)
"""
import csv # a standard Python library
import pyproj
import requests # a third party library also installed on OSGeoLive
# show the location of the folder where all the projection definitions are stored
# On Windows: C:\VirtualEnvs\pytraining\lib\site-packages\pyproj\data
# On Linux /usr/share/proj/epsg
print("Projection definitions folder location : {}".format(pyproj.pyproj_datadir))
url = "http://data-osi.opendata.arcgis.com/datasets/059f9de9770147a2a18c5e5d710839ad_3.csv"
r = requests.get(url)
# encode the text in UTF (see above)
txt = r.text.encode(r.encoding) # r.encoding is "utf-8"
# strip any whitespace at the start and end of the text which could cause extra blank records
txt = txt.strip()
# now split the text by new line characters \n, to return a list of lines
rows = txt.split("\n")
# get an "iterable" from the text, with a record for each row in the list
records = csv.reader(rows)
next(records, None) # skip the headers at the start of the CSV file
# we'll create some projection objects here, outside of the loop
proj1 = pyproj.Proj(init='epsg:2157') # ITM
proj2 = pyproj.Proj(init='epsg:4326') # WGS84
# loop through the records
for r in records:
name = r[2] # access values by list index
x = r[3]
y = r[4]
# print out the ITM coords
print("{} (X: {} Y: {})".format(name, x, y))
# now let's reproject to Lat Lng
x2, y2 = pyproj.transform(proj1, proj2, x, y)
# print out the new coords
print("{} (Lng: {} Lat: {})".format(name, x2, y2))
