Python is an easy-to-use programming language which, thanks to a growing number of cool extension modules, is really taking off in the world of scientific data handling. The Proj4 libraries are a set of programs for performing coordinate system transformations. Both are open source, so you are free to install them on as many computers as you want and to share them with your friends. I had been using both for a while but only recently discovered the pyproj module that performs coordinate transformations inside Python itself. It makes life very easy.
Pyproj can be installed via the pip package manager (pip install pyproj).
1.) Setting up coordinate systems
The first step is to define pyproj ‘objects’ to represent the coordinate systems that you want to use. These can be defined using the Proj notation (see Proj4 documentation for details) but it is easier to set up commonly-used projections by referring to their standard code numbers. These are called EPSG codes and can be looked up on spatialreference.org.
import pyproj # Import the pyproj module
# Define a projection with Proj4 notation, in this case an Icelandic grid
isn2004=pyproj.CRS("+proj=lcc +lat_1=64.25 +lat_2=65.75 +lat_0=65 +lon_0=-19 +x_0=1700000 +y_0=300000 +no_defs +a=6378137 +rf=298.257222101 +to_meter=1")
# Define some common projections using EPSG codes
wgs84=pyproj.CRS("EPSG:4326") # LatLon with WGS84 datum used by GPS units and Google Earth
osgb36=pyproj.CRS("EPSG:27700") # UK Ordnance Survey, 1936 datum
UTM26N=pyproj.CRS("EPSG:32626") # UTM coords, zone 26N, WGS84 datum
UTM27N=pyproj.CRS("EPSG:32627") # UTM coords, zone 27N, WGS84 datum
UTM28N=pyproj.CRS("EPSG:32628") # ... you get the picture
Note: older versions of pyproj use pyproj.Proj("+init=EPSG:4326") syntax and return coordinates in lon, lat order.
2.) Changing between different coordinate systems
In most cases, you will want to change between coordinate systems. This is even the case with GPS or GoogleEarth data, which use the specific WGS84 datum. Coordinate system changes are done with the transform function.
pyproj.transform(wgs84, isn2004, 63.983, -19.700) # (1665725.2429655408, 186813.3884751596)
And when you have lots of data, transformations can be done with lists/tuples/arrays:
lon = [-19.5, -19.7, -19.9] lat = [63.183, 63.583, 63.983] xx, yy = pyproj.transform(wgs84, isn2004, lat, lon) xx # [1674812.314071126, 1665231.44553604, 1655933.043340466] yy # [97526.59264212535, 142220.3063699659, 186937.313488422]
It’s a simple as that.
For other systems, check out the pyproj website. If you are playing with coordinate transforms, then it is likely that at some point you are going to want plot stuff on a map. Python can do maps, too; check out the Cartopy, Fiona and Shapely libraries. Taking things even further, see the following of how to do GIS analysis using only Python with the Geopandas library: https://github.com/BritishGeologicalSurvey/geopandas-demo
British National Grid and the OSGB36 datum.
People working in the UK may have to go through another step to convert their data into the British National Grid format (BNG), which uses two-letter codes to define 100km-wide square regions instead of presenting the full 13-figure coordinates. Thus Arthur’s Seat, an extinct volcano in the centre of Edinburgh, has a BNG grid reference of NT2755072950.
The original version of this blog post contained a code snippet that converted between alphanumeric grid references and OSGB36 coordinates. It has since been converted into a standalone Python module. For instructions on how to install and use it, see https://pypi.org/project/bng/.
Happy mapping!
Update 2019-02-26 to use pyproj version 2.5.0 syntax.
Hi,
Recently I got stuck while transforming an everest transverse mercator Projected Coordinate System which has GCS_Everest as GCS, into WGS84(a vector file direcly digitized and obtained from Google Earth). The problem seemed really an
easy one to solve. I undefined one of the projections and redefined it and
later reprojected as the desired one. However, got no results. I have tried
numerous other ways to do it.
The problem here is I have two shapefiles in two different projections. If
I had raster images, it was not a problem.
I have been doing some reading on EPSG codes and trying to find an
application to run on Python.Any help would be appreciated
I found these projections on Spatial Reference.org.
ESRI:37006: GCS Everest Modified 1969
ESRI:37202: GCS Everest Bangladesh
ESRI:37203: GCS Everest India Nepal
You can click each to get the proj4 details, e.g.
+proj=longlat +a=6377276.345 +b=6356075.41314024 +no_defs
You can use GDAL (www.gdal.org) to reproject different data types. It contains the OGR package which is used on vector data. You can also control GDAL via Python, although I haven’t tried myself, yet.
https://pypi.python.org/pypi/GDAL/
Hope that helps
Thanks for the prompt reply. Trying to figure out things. Hope it gets fixed:-)
Hi Ria!
I think I could use your help. I live in Nepal too and I’m stuck with the same problem.
I’d like to propose you have a bug in your code (and a fix for it). I really think that
np.round()is the wrong operation at line 141 in the code. My justification is twofold:1) See the link provided, and watch the video. The BNG notation does not really identify points. Rather, it identifies squares. So you shouldn’t “round” a point to the “next nearest point”; rather, you should “floor” the value to find the grid square that contains the point — since that grid square is identified by it’s southwest (lower left) corner.
2) The code (as written as of this posting) produced non-sensical results for certain inputs, such as
from_osgb36((529900, 199900), nDigits=4). This results inTQ30100. The “100” portion of “TQ30100” is clearly wrong for nDigits = 4. I suggest the correct answer — using floor() — isTQ2999. However, if you insist on rounding, note that “bumps” you to the next tile north, an the correct BNG would beTL3000(note “TL” instead of “TQ”).Thanks Dan, well spotted. As of now, the code has been updated with your fix.
Pingback: A robot’s guide to fieldwork | Volcan01010
Pingback: All the software a geoscientist needs. For free! | Volcan01010
Pingback: Reprojecting elevation dataset for "tilted Earth" | DL-UAT
Pingback: All the software a geoscientist needs. For free! | The Spatial Blog
Pingback: change letter case python - Ais.gencook.com
Hi,
how would you go about converting lat, long, alt to x, y, z where you have an origin for x/y coordinate (specific place in London)? In other words, I have a lat long and alt data and would like to convert the to x,y,z with respect to a specific place, not the Earth center.
Thanks,
Mike
For anyone using this code on Python 2.* to find Grid References, note that this script will error “Coordinate location outside UK region” with valid OSGB36 coordinates due to lines 127 and 128.
Python 2.* outputs floats when evaluating these expressions, which then cannot be used to reference the correct letter-pair in the _regions array.
SIMPLE FIX is to pass through the inbuilt int function:
x_box = int(np.floor(x / 100000.0)) # Convert offset to index in ‘regions’
y_box = int(np.floor(y / 100000.0))
other options include using the asdatatype features in numpy.
Best,
Krishna