[R-sig-Geo] OS Grid conversion

[Barry Rowlingson]

On Sat, Nov 27, 2010 at 2:12 AM, Scott Abramson <sabramso at princeton.edu> wrote:
> Hi all,
> I have a vector of UK OS National Grid Reference points (ST7475 for
> example). I would like to convert them to WSG84  latitude/longitude
> but cannot figure out a simple way to do this.
>
> I apologize if there is an immediate answer that I have overlooked.

 No time for a full working answer right now, but the function you
need is spTransform from rgdal (or is it in sp or maptools, I forget -
get all of them). You then need to know the offsets for the letter
codes (see the wikipedia page on OSGB national grid) and then you need
to add offsets to your eastings and northings corresponding to the
letter code. Extract parts of your character coordinate "ST7475" using
the substr function.

 Having done all that, you then need to know that OSGB is EPSG:27700 -
put your points in a SpatialPointsDataFrame with that CRS and then
convert to EPSG:4326 (which is WGS84).

 So for example:

 ST7475

ST is 3 east and 1 north from the origin:
http://en.wikipedia.org/wiki/Ordnance_Survey_National_Grid

 So you have (374, 175) (assuming I've got the order right)

 Now let's see... the letter square is 100km, so you've got 3 digits,
thats 1km references, so multiply by 1000 to get metres:

 (374000,175000)

 those should be epsg:27700 coordinates - here's how we transform:

# make a point:

 > require(rgdal);require(sp)
 > pt = data.frame(x=374000,y=175000)
 > coordinates(pt)=~x+y
 > pt
SpatialPoints:
          x      y
[1,] 374000 175000
Coordinate Reference System (CRS) arguments: NA

 okay, now say its osgb grid coords:

 > proj4string(pt)=CRS("+init=epsg:27700")
 > pt
SpatialPoints:
          x      y
[1,] 374000 175000
Coordinate Reference System (CRS) arguments: +init=epsg:27700
+proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000
+y_0=-100000 +ellps=airy +datum=OSGB36 +units=m +no_defs
+towgs84=446.448,-125.157,542.060,0.1502,0.2470,0.8421,-20.4894

 now we have a non NA CRS. So we can spTransform:

 > spTransform(pt,CRS("+init=epsg:4326"))
SpatialPoints:
             x        y
[1,] -2.375734 51.47336
Coordinate Reference System (CRS) arguments: +init=epsg:4326
+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs +towgs84=0,0,0

 How's that? I have to go chop a tree down now before it gets dark...

Barry