[R-sig-Geo] Converting longitude and latitude to UK OS references

Roger Bivand Roger@B|v@nd @end|ng |rom nhh@no
Sun Oct 23 16:33:49 CEST 2022 

Sorry for top-posting, travelling:

Set the environment variable PROJ_NETWORK=ON, and start R, or Sys.setenv("PROJ_NETWORK"="ON") before loading sf:

library(sf)
station_pts<-read.csv(text="Station_ID, Lat, Lon
1, 55.8533, -2.38588
2, 55.8551, -2.39620
3, 55.8473, -2.11304")
stations<-st_as_sf(station_pts, coords=c("Lon", "Lat"), crs="OGC:CRS84")
sf_proj_network()
sf_proj_network(FALSE)
o <- sf_proj_pipelines("OGC:CRS84", "EPSG:27700")
o$definition
o$accuracy
stations_OS_no_cdn <- st_transform(stations, "EPSG:27700", pipeline=o$definition[1])
stations_OS_no_cdn
sf_proj_network(TRUE)
stations_OS <- st_transform(stations, "EPSG:27700")
stations_OS

There is about a 1m difference between the transformation using a grid and the best Helmert transformation. My guess is that even after other corrections, the online point-at-a-time transformation pages use a grid rather than a Helmert transformation. There is usually a reason for observed differences, and https://cdn.proj.org is a very useful resource where mapping agencies have contributed grids.

Roger

--
Roger Bivand
Emeritus Professor
Norwegian School of Economics
Postboks 3490 Ytre Sandviken, 5045 Bergen, Norway
Roger.Bivand using nhh.no

________________________________________
From: R-sig-Geo <r-sig-geo-bounces using r-project.org> on behalf of Micha Silver <tsvibar using gmail.com>
Sent: 23 October 2022 13:58
To: Nick Wray; r-sig-geo using r-project.org
Subject: Re: [R-sig-Geo]  Converting longitude and latitude to UK OS references

Just to combine Rob's and Roger's responses:


library(sf)


# Create a data.frame from text

station_pts<-read.csv(text="Station_ID, Lat, Lon
+ 1, 55.8533, -2.38588
+ 2, 55.8551, -2.39620
+ 3, 55.8473, -2.11304")


# Make an sf object from the DF

# Careful about the order of Lon/Lat

stations<-st_as_sf(station_pts, coords=c("Lon", "Lat"), crs=4326)


# CRS transform

stations_OS <- st_transform(stations, 27700)


str(stations_OS)
Classes ‘sf’ and 'data.frame':  3 obs. of  2 variables:
  $ Station_ID: int  1 2 3
  $ geometry  :sfc_POINT of length 3; first list element:  'XY' num
375940 662301
  - attr(*, "sf_column")= chr "geometry"
  - attr(*, "agr")= Factor w/ 3 levels "constant","aggregate",..: NA
   ..- attr(*, "names")= chr "Station_ID"


# Just to be sure, display CRS:

st_crs(stations_OS)
Coordinate Reference System:
   User input: EPSG:27700
   wkt:
PROJCRS["OSGB 1936 / British National Grid",
     BASEGEOGCRS["OSGB 1936",
         DATUM["OSGB 1936",
             ELLIPSOID["Airy 1830",6377563.396,299.3249646,
                 LENGTHUNIT["metre",1]]],
         PRIMEM["Greenwich",0,
             ANGLEUNIT["degree",0.0174532925199433]],
         ID["EPSG",4277]],
     CONVERSION["British National Grid",
         METHOD["Transverse Mercator",
             ID["EPSG",9807]],
         PARAMETER["Latitude of natural origin",49,
             ANGLEUNIT["degree",0.0174532925199433],
             ID["EPSG",8801]],
         PARAMETER["Longitude of natural origin",-2,
             ANGLEUNIT["degree",0.0174532925199433],
             ID["EPSG",8802]],
         PARAMETER["Scale factor at natural origin",0.9996012717,
             SCALEUNIT["unity",1],
             ID["EPSG",8805]],
         PARAMETER["False easting",400000,
             LENGTHUNIT["metre",1],
             ID["EPSG",8806]],
         PARAMETER["False northing",-100000,
             LENGTHUNIT["metre",1],
             ID["EPSG",8807]]],
     CS[Cartesian,2],
         AXIS["(E)",east,
             ORDER[1],
             LENGTHUNIT["metre",1]],
         AXIS["(N)",north,
             ORDER[2],
             LENGTHUNIT["metre",1]],
     USAGE[
         SCOPE["Engineering survey, topographic mapping."],
         AREA["United Kingdom (UK) - offshore to boundary of UKCS within
49°45'N to 61°N and 9°W to 2°E; onshore Great Britain (England, Wales
and Scotland). Isle of Man onshore."],
         BBOX[49.75,-9,61.01,2.01]],
     ID["EPSG",27700]]


If you have your station Long/Lat coordinates in a CSV file, then
`st_read()` will read directly into an sf object.


On 23/10/2022 12:08, Nick Wray wrote:
> I am trying to convert longitude and latitude values for UK weather
> stations to UK Ordnance Survey National Grid References.  There are sites
> where one can do them one at a time but I have a large number.  I have
> found some code which does the conversion and include the first three
> points which I want to convert as an example
>
> library(sp)
> xy<-as.data.frame(cbind(c(55.8533,55.8551,55.8473),c(-2.38588,-2.39620,-2.11304)))
> colnames(xy)<-c("lon","lat")
> xy
> coordinates(xy)<-~lon+lat
> ## see site in jounral convert llon lat
> proj4string(xy)<-CRS("+init=epsg:4326")
> ptsOS<-spTransform(xy,CRS("+init=epsg:27700"))
> ptsOS
>
> but it doesn't give the right answers - the output doesn't correspond to
> either what I get whn I put an individual point into a conversion site nor
> with other points in the same region which I have got from other sources
> and are mutually compatible
>
> Can anyone explain where I'm going wrong and how I can fix this?  Thanks
> Nick Wray
>
>       [[alternative HTML version deleted]]
>
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo using r-project.org
> https://eur02.safelinks.protection.outlook.com/?url=https%3A%2F%2Fstat.ethz.ch%2Fmailman%2Flistinfo%2Fr-sig-geo&data=05%7C01%7CRoger.Bivand%40nhh.no%7C68be57c5516345d3e61a08dab4edf0ee%7C33a15b2f849941998d56f20b5aa91af2%7C0%7C0%7C638021231300728219%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C3000%7C%7C%7C&sdata=mqCLcKYV2K%2FuxI3wmkhC6QMbALy5Y7uuVXBHkVoyAxE%3D&reserved=0

--
Micha Silver
Ben Gurion Univ.
Sde Boker, Remote Sensing Lab
cell: +972-523-665918
https://eur02.safelinks.protection.outlook.com/?url=https%3A%2F%2Forcid.org%2F0000-0002-1128-1325&data=05%7C01%7CRoger.Bivand%40nhh.no%7C68be57c5516345d3e61a08dab4edf0ee%7C33a15b2f849941998d56f20b5aa91af2%7C0%7C0%7C638021231300728219%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C3000%7C%7C%7C&sdata=5IJ6yY1pJaYtRglQLNV%2FJY3SfktRn0864FPres66CSI%3D&reserved=0

_______________________________________________
R-sig-Geo mailing list
R-sig-Geo using r-project.org
https://eur02.safelinks.protection.outlook.com/?url=https%3A%2F%2Fstat.ethz.ch%2Fmailman%2Flistinfo%2Fr-sig-geo&data=05%7C01%7CRoger.Bivand%40nhh.no%7C68be57c5516345d3e61a08dab4edf0ee%7C33a15b2f849941998d56f20b5aa91af2%7C0%7C0%7C638021231300728219%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C3000%7C%7C%7C&sdata=mqCLcKYV2K%2FuxI3wmkhC6QMbALy5Y7uuVXBHkVoyAxE%3D&reserved=0

More information about the R-sig-Geo mailing list