Using the GDAL library to convert spatial coordinates from Northings and Eastings to WGS84

I have written before about geospatial data types and the open source libraries which can be used to manipulate them.

Today, I needed to take some geospatial data which was in a database with northing and easting fields (basically ordinance survey gridpoints), and convert it to a WGS-84 coordinate system (latitude, longitude used by most GPS systems these days). This is so that I could spatially match it to another data-set.

The conversion took a few steps and had a couple of ‘gotchas’, so I’ve documented it here for prosperity. Although I’m converting from SRID 27700 to 4326 here, you can of course use the same steps to convert between other types.

The steps below assume you have already installed the GDAL suite.

• Get your source data into to a plain old csv including at least a primary key and the easting and northing data fields. Include a your field names on the top row, and let’s assume they’re called ID, Easting andNorthing.
• Place the data in (C:\Example\Location\) as SourceData.csv
• Create a new text file in the same directory, called SourceData.vrt, and enter the following:

  <OGRVRTDataSource>
      <OGRVRTLayer name="EastingNorthing">
          <SrcDataSource>C:\Example\Location\SourceData.csv</SrcDataSource>
          <GeometryType>wkbPoint</GeometryType>
          <LayerSRS>EPSG:27700</LayerSRS>
          <GeometryField encoding="PointFromColumns" x="Easting" y="Northing"/>
      </OGRVRTLayer>
  </OGRVRTDataSource>

• Open the ‘GDAL command prompt’
  1. Enter the following to ensure your vrt is working. You should see information about the file come back, including ‘Number of features’, which should be equal to the number of rows in your source file, and then each of the features.

     ogrinfo -ro -al "C:\Example\Location\SourceData.vrt"

  2. If the number of features is 0, or you get some other error, check everything.
     • Are there stray commas in the file?
     • Are the field names correct, and match the vrt specification?
     • Are the file addresses right?
  3. You should also see, for each of the features, the easting and northing that has been picked up. Make surethat these numbers haven’t been truncated. In my case, my source data had commas and GDAL had truncated the numbers after the comma (e.g. 234,567 -> 234). I didn’t notice at first, and was confused when all my converted coordinates were clustered in the sea to the South West of the UK… I resolved this by simply removing the extra commas from the source csv, although I’m sure there are ways of doing this with the vrt specification as well…
  4. Assuming the above works, you can go ahead and run the conversion. It’s easiest to drop your output file into a new folder, as otherwise it causes some name clashes with the source file. Create a new folder (C:\Example\Location\Output\), and run the following:

     ogr2ogr -s_srs EPSG:27700 -t_srs EPSG:4326 -f CSV -lco "GEOMETRY=AS_WKT" "C:\Example\Location\Output\output.csv" "C:\Example\Location\SourceData.vrt"

     (s_srs and t_srs set the source and target geography projection; lco instructs GDAL to include a field in the output with the new coordinates; and finally, f sets the target filetype to csv.)

  5. Check your output field. You should see a new field containing data like 'POINT (0.000428633570006 52.700838815264168)'. This is your shiny new WGS-84 output, in the 'well known text' format.

  6. If you read this into your database as a varchar field called “WKT”, you will be able to convert it to a spatial data type with something like geography::STPointFromText(“WKT”, 4326), depending on your tool of choice! Note that you also have to specify the ‘4326’ to specify the data is in the WGS-84 format!

     Select "ID", geography::STPointFromText(replace("WKT",'"',''),4326) AS "Coordinate"
     into clean.LoadedData
     from raw.LoadedData

  7. Do some spatial joining. In my case, I had two sets of spatial data, it was acceptable to brute force the distance calculation for all combinations of records, and save any results where the distance was less than 50m. That way, I got all the expensive spatial calculations done up front, and could then work with the results table:

     select a.ID AS "A_ID",
          b.ID AS "B_ID",
          a.Coordinate.STDistance(b.Coordinate) as "Distance"
     into working.CloseNodes
     from clean.LoadedData a
     inner join clean.ReferenceData b
     on a.Coordinate.STDistance(b.Coordinate) < 50