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. The easiest way to get pyproj is as part of the Matplotlib Basemap package. On Ubuntu Linux, this can be installed with the command: sudo apt-get install python-mpltoolkits.basemap 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. 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 can be installed with a single command on recent Linux systems (sudo apt-get install python-pyproj). 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 rest of the Matplotlib Basemap module. 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. Below is a Python module that carries out the second step for you. Save the code as ‘BNG.py’ in a directory where Python can find it. Note: The script below was updated to fix a rounding error on 4 October 2013. If you downloaded the code before then, you should replace it with the updated version below. The function converts between a BNG grid reference as a string, and a tuple of OSGB36 (x,y) coordinates. When converting to BNG coordinates, there is an opportunity to specify how many figures to use. Note: the coordinate transform between WGS84 and OSGB36 is complicated by some distortions in the British National Grid. This is mainly important for high-precision work such as surveying or construction, and is described by the National Grid Transformation, OSTN02. To get the most accurate results, give Proj a datum shift grid file that describes the offsets. 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 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. 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”). Source.