Easily change coordinate projection systems in Python with pyproj

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.

Categories: Uncategorized

13 Comments

  1. Ria says:

    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

    • Asphota says:

      Hi Ria!
      I think I could use your help. I live in Nepal too and I’m stuck with the same problem.

  2. Dan Harasty says:

    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 in TQ30100. The “100” portion of “TQ30100” is clearly wrong for nDigits = 4. I suggest the correct answer — using floor() — is TQ2999. However, if you insist on rounding, note that “bumps” you to the next tile north, an the correct BNG would be TL3000 (note “TL” instead of “TQ”).

  3. Pingback: A robot’s guide to fieldwork | Volcan01010

  4. Pingback: All the software a geoscientist needs. For free! | Volcan01010

  5. Pingback: Reprojecting elevation dataset for "tilted Earth" | DL-UAT

  6. Pingback: All the software a geoscientist needs. For free! | The Spatial Blog

  7. Pingback: change letter case python - Ais.gencook.com

  8. Mike says:

    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

  9. Krishna says:

    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